Reconstructing input from the output of a known discrete-time IIR filter

In general, is it possible to design such a filter?

If by "in general" you mean "for any arbitrary IIR filter", then the answer is no.  In fact, for any real filter, retrieving the actual $x(k)$ in a useful way is problematic, at best.
The basic problem is that a filter filters information out, and leaves you with a residue.  At some point, the initial information is simply lost, and you cannot retrieve it.
About the only time that an inverse filter will work is if you're dealing with some natural phenomenon that shapes the spectrum of $x(k)$ but does not significantly remove any information that you need.

If not then under what conditions is it impossible to retrieve the input $x(k)$?

Mathematically, if you have a filter $H(\omega)$, and it is zero for one or more values of $\omega$, then you simply cannot reconstruct any components of $x(n)$ at those frequencies.
Moreover, if the filter defined by $y(k) = h\left ( x(k) \right)$ has pure delay or unstable zeros (i.e., if it is not minimum phase), then $x(k)$ cannot be reconstructed.  There are techniques for coming up with a best guess of the signal, however.
Practically, the filtering problem is never $y(k) = h\left ( x(k) \right)$.  It is always $y(k) = h\left ( x(k) \right) + n(k)$, where $n(k)$ is some noise process.  So where the filter's response comes close to zero in the frequency domain, you cannot practically reconstruct $x(k)$.

I have known discrete-time IIR filter. I pass a signal x(n) through the filter and obtain the output y(n). How would I design a filter to obtain x(n) from y(n)?

In increasing order of both difficulty and possibility of success:

Cheat, and just measure x(n).
Calculate the transfer function, $H(\omega)$, for your filter.  If the filter is minimum phase, then $\dfrac{1}{H(\omega)}$ will be realizable and stable.  Implement that.  Your output will be $\hat X(\omega) = H(\omega) \left(Y(\omega) N(\omega) \right) $ -- which means that at any frequency where the original filter really attenuates the signal, your estimate of $x(k)$ will be noisy.
Study up, and implement a Kalman filter to estimate $\hat x(k) \simeq x(k)$ from $y(k)$.  If you know the filter structure exactly, this is probably the best way to get an estimate $\hat x(k)$ that is optimal in the least-squares sense.  Assuming your Kalman filter is stable, and the turn-on transient is acceptable, you can derive a steady-state Kalman filter from it.

A corollary: there's probably some variant of the Wiener filter that would do this, too -- but it's been a very long time since I've needed to understand those critters, so I can't say for sure.  Certainly, a steady-state Kalman and a properly constructed Wiener filter -- if it exists -- would basically be the same thing, since they're solving the same problem.
If mean-squared error is not the thing you need to minimize, or if you don't know the filter structure exactly, then some other optimal state estimation method will be your friend.  I'd start with H-infinity filters, but I'm not sure that's where I'd end up.
Yet another wrinkle in this is that -- especially if the system isn't minimum-phase -- if you can get by with a delayed estimate (i.e., you're estimating $x(k - k_0)$ for some delay $k_0 > 0$) then it would be better than any estimate of $x(k)$, and this, too, would drop out of a steady-state Kalman or Wiener filter design exercise.

If you take that last course, note that especially if you're dealing with a non-minimum-phase filter you will be able to get much better estimates of delayed versions of $x(k)$.  I.e., if there's some delay $\kappa$ that is perfectly acceptable, then you can make a Kalman estimator that finds $\hat x(k - \kappa) \simeq x(k - \kappa)$ with less (possibly much less) error than any filter could find $\hat x(k) \simeq x(k)$


How do you send signals through geological formations

For digital communications, bit error rate performance is all about $\frac{E_b}{N_0}$, the energy per bit divided by the noise power spectral density. The part that you can control is the energy per bit. Energy per bit is equal to $\text{(average power)(bit duration)}$. You can increase $E_b$ (and therefore decrease your bit error rate) in two ways:

Raising your transmitter power. This may or may not be possible based on the hardware you're using, the frequency band in use, and any governmental regulations for your locale.
Decreasing the bit rate (i.e. give more time per bit). In cases of extremely high path loss, like communication into caves, you can get reliable operation at a very low bit rate (although other issues such as long-term oscillator stability can come into play and diminish the gain that you might expect).

Communication into deep geological formations is something of a niche application, so you may have trouble finding a lot of literature on suggested approaches. There are companies who specialize in making this type of equipment, so it may serve you better to consult the marketplace for existing solutions that might be off the shelf instead of trying to engineer one from scratch (if your budget permits, of course; I expect such systems to be expensive due to the rather small market).


Filter Banks and FFT, are they similar conceptually?

The answer is yes, the discrete Fourier transform (DFT), of which the fast Fourier transform (FFT) is an efficient implementation, can be interpreted quite literally as a filter bank. Recall the formula for the DFT:

$$
X[k] = \sum_{n=0}^{N-1} x[n] e^{\frac{-j2\pi k n}{N}}
$$

The "filterbank" interpretation for the above equation goes something like this:

In order to get the $k$-th DFT output $X[k]$, we start with the input signal $x[n]$.
Multiply $x[n]$ by a complex exponential function $e^{\frac{-j2\pi k n}{N}}$. By the frequency-shifting property of the DFT, this has the effect of shifting the content at frequency $\frac{2\pi k}{N}$ down to zero frequency.
Convolve the result with a filter that has the impulse response:
$$
h[n] = \begin{cases}
1, && 0 \le n \le N-1 \\
0, && \text{otherwise}
\end{cases}
$$
This filter (known as a moving average or a "boxcar") has a lowpass characteristic.

Thus, we can view the DFT as placing a bank of uniformly-spaced critically-sampled filters across the signal's spectrum. Each filter has the same frequency response as an $N$-sample moving average filter, centered at frequencies $e^{\frac{j2\pi k n}{N}}, k = 0, \ldots, N-1$. I described this similarly in an answer to a previous question here.


Does the inverse-CTFT exist for a dirac delta?

Yes, this is a complex exponential $e^{2 \pi i f_0 t}$, at a frequency determined by the delta's "position" $f_0$ (your input being $\delta(f - f_0)$). Write the integral for the inverse Fourier transform, use the definition of $\delta$ and you'll see it "selects" at this particular frequency the complex exponential being integrated.


Multiplication in the wavelet domain, what does it look like in real space?

They can be  used for time-varying filtering of non-stationary signals, under names like Gabor multipliers, STFT filters for instance.  You can find a review reference in the chapter Linear Time-Frequency Filters, 11.1.3 Implicit design and generalization in Spreading function representation of operators and Gabor multiplier approximation.

You also have the 2004 paper The convolution theorem for the continuous wavelet transform. My  experience in this domain is limited to localized matched filtering, by  multiplying coefficients by scalar complex numbers to locally modify amplitude and phase (with complex continuous wavelet frames), see Adaptive multiple subtraction with wavelet-based complex unary Wiener filters.

The redundancy of the continuous transform certainly helps a lot. With DWT (a variant of them) I have also used  masks to remove localized high amplitude disturbances in some discrete wavelet subbands. The masks were initially binary ($0$ on the places to remove, $1$ outside), but they were tapered to have a smooth profile between $0$ and $1$. Not very theoretical, but efficient to separate directional data.

Finally, you can find across-scale multiplication for ROI (region-of interest) coding with JPEG 2000.


Image processing research areas using GPU

I've done a bit of GPGPU (for optical flow) during my PhD, so here are a few hints to help you:

choose a computation intensive task. You're very likely to find them in variational inverse problems (especially when expressed images instead of sparse feature sets). They are very common tasks: image deconvolution, compressed sensing, dense optical flow (and more broadly dense image registration...), dense 3D reconstruction, image segmentation... ;
from my point of view, non-local image denoising (using non-local means) is very easy to understand, and it's computation-intensive enough to have been included in NVidia's Cuda examples for quite a while. It really deserves a new, better GPU implementation (such as the GPU-based kd-trees proposed a few years ago). However, it may prove challenging for an undergrad;
a useful site to see what Image Processing tasks can be done on GPUs: the Austrian (T.U. Graz) lab GPU4Vision.


How do I detect distinct objects when their edges touch each other?

Detecting different components:
If you're trying to detect the different components, there probably are other approaches to do them than detecting the contours. Here's an example in Mathematica. An erosion followed by dilation is used to close the gap in the second component before detection (if you don't do this, it won't detect it).
img = Binarize@Import["http://i.stack.imgur.com/yqDyu.png"];
Colorize[MorphologicalComponents[Dilation[Erosion[img,1],1]]]

The figure on the left below, shows imperfect object detection (without closing the gap) and on the right, shows the correct detection (running the above code).

Detecting the different contours:
However, if you indeed want to separate only the contours, here's an example. The erosion and dilation are performed as before to close the gap and the resulting image is run through a Canny edge detector. I've made the default options explicit, so that you can see what's being used.
img2 = EdgeDetect[Dilation[Erosion[img, 1], 1], Method -> "Canny"]

This will give you both the inside and the outside edge (see figure on the left below), since the pixel width is greater than 1 all around. I haven't had much luck trying to get it thinner, as the performance degrades (may be different for your other images). The inner contours are the ones you want, and the outer contour is just the combined contour of all 4 components. Now all we need to do is drop the outermost one with:
SelectComponents[img2, "EnclosingComponentCount", # > 0 &]

which gives you just the inner contours (see bottom right). In other words, it only picks those contours that are enclosed by at least one other contour, which automatically disqualifies the outermost one. I do not know the equivalent of these commands/operations in openCV.

Note that the apparent breaks in the figure are due to saving to jpeg in a smaller size. It does not look that way on my screen.


Mel filter bank is important due to following reasons:

It applies the Mel-frequency scaling, which is perceptual scale that helps to simulate the way human ear works. It corresponds to better resolution at low frequencies and less at high.
Using the triangular filter-bank helps to capture the energy at each critical band and gives a rough approximation of the spectrum shape, as well as smooths the harmonic structure. In theory you could manipulate on raw DFT bins, but then you are not reducing the dimensionality of your features - this is the whole point of doing filter-bank analysis, to capture the spectral envelope.

Obviously warping of frequency scale can be altered. Other mappings are possible, such as Bark, linear, etc. Mel is just one out of many options. To finalize, you can always calculate more filter-banks to have coarser scale (getting closer to actual DFT), but you will notice that at some point it will degrade the performance and your features do not carry essential information.


Fastest Available Algorithm to Blur an Image (Low Pass Filter)

The fastest blur would be Box Blur.
You can implement it using Running Sum.
I think Intel FilterBoxBorder works in that manner.

If you'd like you can do a few passes of it to approximate the Gaussian Blur.

You can also use IIR Filter Coefficients to blur the image quite easily.

You may have a look at my project Fast Gaussian Blur.


understanding bases and frames for Gabor transform

The frames end up giving you an orthogonal representation, just not exactly a basis representation that you're talking about. Simply by definition that's not what a Gabor filter does. It's a time-frequency representation of your signal, while "orthogonal basis vectors" wouldn't have a time-specific property -- they stretch for as long as your signal stretches.

If you want to use orthogonal basis transforms, use Fourier, cosine, empirical mode decomposition and similar transforms.

If you want time-frequency representation, your can use short-time Fourier transform, wavelet transforms (of which Gabor can be considered a special case) and other similar transforms.


Savitzky-Golay smoothing filter for not equally spaced data

Because of the way the Savitzky-Golay filter is derived (i.e as local least-squares polynomial fits), there's a natural generalization to nonuniform sampling--it's just much more computationally expensive.

Savitzky-Golay Filters in General

For the standard filter, the idea is to fit a polynomial to a local set of samples [using least squares], then replace the center sample with the value of the polynomial at the center index (i.e. at 0). That means the standard SG filter coefficients can be generated by inverting a Vandermonde matrix of sample indicies. For example, to generate a local parabolic fit across five samples $y_0\dots y_4$ (with local indicies -2,-1,0,1,2), the system of design equations $Ac = y$ would be as follows:

$$
\begin{bmatrix}-2^0 & -2^1 & -2^2 \\ -1^0 & -1^1 & -1^2 \\ 0^0 & 0^1 & 0^2 \\ 1^0 & 1^1 & 1^2 \\ 2^0 & 2^1 & 2^2 \end{bmatrix} \begin{bmatrix} c_0 \\ c_1 \\ c_2 \end{bmatrix} = \begin{bmatrix} y_0 \\ y_1 \\ y_2 \\ y_3 \\ y_4 \end{bmatrix}.
$$

In the above, the $c_0 \dots  c_2$ are the unknown coefficients of the least squares polynomial $c_0 + c_1x + c_2x^2$. Since the value of the polynomial at $x = 0$ is just $c_0$, computing the pseudoinverse of the design matrix (i.e. $c = (A^TA)^{-1}A^T y\space $) will yield the SG filter coefficients in the top row. In this case, they would be

$$
\begin{bmatrix}c_0 \\ c_1 \\ c_2 \end{bmatrix} =
\begin{bmatrix}
-3 & 12 & 17 & 12 & -3 \\
-7 & -4 & 0 & 4 & 7 \\
 5 & -3 & -5 & -3 & 5 \\
\end{bmatrix}
\begin{bmatrix} y_0 \\ y_1 \\ y_2 \\ y_3 \\ y_4 \end{bmatrix}.
$$

Note that since the derivative of $c_0 + c_1x + c_2x^2$ is $c_1 + 2c_2x$, the second row of the matrix (which evaluates $c_1$) will be a smoothed derivative filter. The same argument applies for successive rows--they give smoothed higher-order derivatives. Note that I scaled the matrix by 35 so the first row would match the smoothing coefficients given on Wikipedia (above). The derivative filters each differ by other scaling factors.

Nonuniform Sampling

When the samples are evenly spaced, the filter coefficients are translation-invariant, so the result is just an FIR filter. For nonuniform samples, the coefficients will differ based on the local sample spacing, so the design matrix will need to be constructed and inverted at each sample. If the nonuniform sample times are $x_n$, and we construct local coordinates $t_n$ with each center sample time fixed at $0$, i.e.

$$
\begin{align}
t_{-2} & = x_{-2} - x_0 \\
t_{-1} & = x_{-1} - x_0 \\
t_0 & = x_0 - x_0 \\
t_1 & = x_1 - x_0 \\
t_2 & = x_2 - x_0
\end{align}
$$

then each design matrix will be of the following form:

$$
A =
\begin{bmatrix}
 t_{-2}^0 &  t_{-2}^1 & t_{-2}^2 \\
 t_{-1}^0 &  t_{-1}^1 & t_{-1}^2 \\
 t_0^0 &  t_0^1 & t_0^2 \\
 t_1^0 &  t_1^1 & t_1^2 \\
 t_2^0 &  t_2^1 & t_2^2
\end{bmatrix} =
\begin{bmatrix}
 1 &  t_{-2} & t_{-2}^2 \\
 1 &  t_{-1} & t_{-1}^2 \\
 1 &  0 & 0 \\
 1 &  t_1 & t_1^2 \\
 1 &  t_2 & t_2^2
\end{bmatrix}.
$$

The first row of the pseudoinverse of $A$ dotted with the local sample values will yield $c_0$, the smoothed value at that sample.


As Dilip pointed out in the comment above, you can get the impulse response using the inverse Fourier transform. However, a slightly easier method might be to use the Laplace domain instead; it's more amenable to easy inverse transforming via transform tables. First, recall that the frequency response is really just the $s$-plane transfer function evaluated along the imaginary axis. That is:

$$
H(j\omega) = H(s)|_{s=j\omega}
$$

Go through your frequency response expression and substitute $s$ back in for $j\omega$:

$$
\begin{align}
H(j\omega) &= \frac{1}{-\omega^2+3j\omega+2} \\
&= \frac{1}{(j\omega)^2+3j\omega+2} \\
\end{align}
$$

$$
H(s) = \frac{1}{s^2+3s+2}
$$

So now we have the transfer function in the $s$-plane. To get from here to an easily-inverse-transformable expression, we use a partial fraction decomposition on the transfer function (an online calculator is available here:

$$
H(s) = \frac{1}{s+1} + \frac{1}{s+2}
$$

Given this form, we now refer to a table of common Laplace transforms to find:

$$
e^{-\alpha t} u(t) \Leftrightarrow \frac{1}{s+\alpha}
$$

where $u(t)$ is the Heaviside step function. The impulse response of the system is the inverse Laplace transform of its transfer function. Using the linearity of the Laplace transform and the above relation, we can therefore conclude that:

$$
h(t) = u(t) \left(e^{-t} + e^{-2t} \right)
$$


First of all, FT and IFT of a discrete sequence are the same up to complex conjugation/scaling - these two operations will thus yield qualitatively the same result (in the sense that subsequence processing/statistical modeling will make abstraction of these differences anyway). The textbook definition of cepstrum requires a Fourier transform. People sometimes use the inverse Fourier transform to emphasizes that the X axis on the cepstrum plot can indeed be interpreted as a time (rather than as a quefrency).

The MFCC extraction process found in speech and audio analysis systems uses the DCT since it has the advantage of being a real -> real transform, and because the resulting coefficients are in practice decorrelated from each other - which is not the case if you take the sequence of real/imaginary part of a FT or IFT. This makes the resulting data very good for subsequent modeling with gaussian models with diagonal covariance matrix.

Another view on this: take thousands of hours of audio, extract log mel frequency coefficients (through a bank of filters with triangular responses, with center frequencies distributed on the Mel scale) for each analysis frame. You'll get a big dataset... Do a PCA on it to extract the dimensionality reduction transform that captures most of the variance of the mel coefficients. The transform matrix the PCA will produce will be uncannily close to the DCT matrix.


I'm looking to do my senior capstone project in the field of medical tomography, what is a good book to get acquainted with the field?

What about for a start:

Buzug, Computed Tomography
Hermann, Fundamentals of Computerized Tomography


Dominant eigenvectors of an unknown matrix

In the absence of other answers, and sorry if you have already considered this, but the power method can be very fast with a sparse matrix. You may like to check out how Google does it (eigenvectors are used for page rank, which is pretty fast).

http://en.m.wikipedia.org/wiki/PageRank

  Because of the large eigengap of the modified adjacency matrix
  above,[22] the values of the PageRank eigenvector can be approximated
  to within a high degree of accuracy within only a few iterations.

Footnote 22 points to: http: //www-cs-students.stanford.edu/~taherh/papers/secondeigenvalue.pdf

Footnote 23 discusses how to make the calculation even faster
http: //projecteuclid.org/euclid.im/1150474883

I have not read the papers, but I hope that you find one of them to be useful. (I had to break the links as new users can not post many links.)


Cross-correlation - Why is there a complex part after FFT/IFFT when deleting the mirrored part of the spectrum?

If the input signal is real-valued, then the DFT returns a sequence of length $N$ that satisfies $$X[N-n]=X^*[n]$$

As one usually plots the magnitude of the spectrum, then the relationship of interest is:

$$|X[N-n]|=|X^*[n]|\implies|X[N-n]|=|X[n]|$$

When you work with real-valued signals, one tends to "ignore" the second half of the FFT spectrum due to this symmetry property.

However, if the signal is complex, then both sides of the spectrum are fundamental, as the symmetry is not accomplished anymore. When you call the function ifft after discarding one side of the spectrum, the function is assuming that the original FFT (the input) is non-symmetrical (as you deleted one half of it), so the time-domain sequence from which it came from must be complex. That's why you get a complex-valued signal when coming back to the time-domain.


Questions on the generation of R.O.C. curves

Question 1) Yes, you must know what the right answer is before-hand in order to judge the receiver's performance.  I am doing this right now on a project.  I generate data, scramble it, modulate it, add noise and carrier offset, etc., then I give it to the receiver.  The bits I started out with are my "yes/no" labels.

Question 2) All such measurements are statistical in nature.  Noise is statistical.  Jitter is statistical.  A priori channel conditions are statistical.  Of course the ROC curve is also statistical.  How could it be otherwise?

Question 3) I don't quite understand this question.  A picture would help.

EDIT: Receivers are not tested in every set of conditions under the sun.  First, the designers rely on everyone obeying the FCC requirements for a given band, so in a cellular band they assume that there won't be any TV stations, or radar chirps, and that all other cellular signals will be at specified frequencies and power levels, with specified maximum carrier offsets.

The designers then model "typical" environments (i.e. channels).  That implies certain transmission losses due to free-space path loss, antenna gains, multi-path, etc.  They then put together various typical and worst-case scenarios and test it under those conditions.  For example, a worst case interference scenario would be having adjacent signals at the maximum allowed power.

If the performance is as expected for those selected conditions, the designers assume, through "interpolation", that it will perform acceptably in the other non-tested environments.  Conditions can always be worse than the worst-case test scenario of course.  That is when the phone drops the call. :-)


Convolution of two large signals in MATLAB

If your shorther pulse signal is really a rectengular waveform then please look for the other answer, but otherwise for a general waveform pulse the following code snipped (excerpted from Maximillian's previous post) shows the actual results (on my laptop with MATLAB R2015) of timings for a frequency domain vs time domain implementation of the convolution operation in simplest terms.

Note that I have reduced the sampling rate Fs from 10 Ghz to 1Ghz due to memory reasons and also I have slightly adjusted the sequence lengths to result in a power of $2$ FFT length so that FFT could run at full speed. (this is slightly unfair for a practical course but here I would like to see the maximum attainable speed gain, hence it would be ill posed to use inefficient FFT lengths for such a purpose)

Fs = 10E8; % reduced for MEMORY reasons
T = 17e-3; % slightly adjusted
t = 0:1/Fs:T;
t = t(1:(2^24-2^14+1)); % its length is adjusted to create a 2^N factoring later.
y = sin(2*pi*1000*t); % The input signal
L1 = length(y);

w_rect = 17e-6;  % width of rect
R = ones(1,Fs*w_rect);
R = pulse(1:2^14);
L2 = length(R);

tic
y2 = real(ifft(fft(y,L1+L2-1).*fft(R,L1+L2-1)));
tdft = toc

tic
y3 = conv(y,pulse);
tconv = toc

eff = tconv/tdft

The result is the following: tdft = 1.4769 s, tconv = 42.14 s and hence speed gain = eff = 28.5x.

On the other hand, for slightly different but otherwise unfactorable FFT lengths, the speed gains decreased down to 10x as understandable.

This is exactly the situation where you should seriously consider implementing your convolution in the frequency domain (using efficient FFTs), rather than in the time domain!

The condition is triggered by the fact that the shorter pulse signal (with about 150k samples) is already too long for a direct time domain convolution.

Lets make a brief comparison between the two approaches for the signals of lengths about $L_1 = 150 \times 10^6$ samples and $L_2 = 150 \times 10^3$ samples:

In the first case a direct time domain implementation would approximately require $L_1 \times L_2$ many real MACS which is about $(150 \times 10^6) \times (150 \times 10^3) = 2.25 \times 10^{13}$ real MACS.

Now compare this to an FFT based frequency domain implementation (according to the convolution in time domain equals multiplication in frequency domain theorem of DFT). Again using rough numbers there will be three steps: 1- Convert both signals to frequency domain of the same (and sufficiently long) lengths (which is actually $L_1+L_2-1$ )using efficient FFT, 2-Multpily the obtained DFTs and get the output DFT, 3-Convert the result back to the time domain using the inverse FFT.

The first step requires about $2 \times (150 \times 10^6) \times \log_2(150*10^6) \times 4$ real MACS which is about $3.25 \times 10^{10}$ real MACs

The second step would require roughly $4 \times (150 \times 10^6) = 6 \times 10^8$ real MACs

And the final, third, step would require about $1.6 \times 10^{10}$ real MACS.

Hence the totality of the FFT based frequency domain convolution implementation would require roughly: $(3.2 \times 10^{10}) + (6 \times 10^8) + (1.6 \times 10^{10}) = 5 \times 10^{10}$ real MACS.

Now compare this to the time domain implementation which required $2.25 \times 10^{13}$ real MACS. This could potentially result in about 500x speed improvement assuming a basic implementation. So for example if the convolution would take 5 hours in direct time domain implementation then it should take about 36 seconds in DFT based case!

So you better follow the second approach. That being said, however, now current available multicore CPU + GPU based extensively vectorized (parallel) implementations can dramatically reduce the direct time domain computation time of convolutions, however again then the efficiency of FFT will also benefit from this parallelization and therefore you should, nevertheless, see an order of magnitude improvement in speed when using a DFT approach.

PLease also consider that such long sequences can be processed in blocks and that also would account for certain imprevements in memory and speed performances.

EDIT: I've tried to implement the DFT processing in an non-blocked based mode, where the whole signal processed at once in memory but that required too much RAM (about at least 10 GBs for MATLAB alone in my PC which has only 8GBs of RAM for all programs currently so it started paging and went into freezing, (the computer simply stcuk without any response even to ctrl+alt+del, so I had to hard reset it) If you have more memory such as 16 GB or above you may still try and it may probably do ok, as I assumed it would at most require about 10 GBs of RAM (plus any other programs using your remainig RAM). But otherwise you should consider dividing your larger signal (which 150 M samples) into smaller pieces or use some other technique.


FIR filter digital differentiator with low cutoff

What you want is what I call a Jekyl and Hyde design. The low pass is parsimoniously represented in terms of cosines and the derivative in sines. So let Jekyl be Jekyl and Hyde be Hyde.

The way to go very low would be a multirate approach.  Low pass and decimate for a few cascades until the high pass part of your differentiator is a reasonable fraction of the band.  If you need the original sample rate, upsample appropriately.  Your filter plot implies that you really aren't interested in most of your original band, so why retain it?


The required bandwidth depends on the symbol rate (or baud rate). If you have a binary system (i.e. two symbols) then the bit rate and the symbol rate are equal. If you use $M>2$ symbols, you can transmit more than 1 bit per symbol and your symbol rate, hence also the required bandwidth, is decreased. The problem with increasing the number of symbols is that you need more energy to keep the average distance between the symbols (i.e., the robustness against noise) the same. As pointed out by @vaz in a comment, the latter is only true for passband PAM systems, such as QAM or PSK.


Representation, Bandwidth and composition of Digital Signals

First

  (i) I know that analog signals are continuous and digital signals are
  discrete. Analog signal are represented by a continuous time-varying
  quantity like voltage. How are digital signals represented in actual
  practice? I mean to ask how are the 0s & 1s actually
  represented/transmitted over a transmission medium? Are they
  represented by a sudden spike then staying stable and then a sudden
  fall of the time-varying quantity like voltage?

Physical representation of discrete-time (or digital) signals is also through analog waveforms, but the information they carry is discrete. They can be transmitted like pulses of $5$ V, $0 $ V in a digital computer, or sine waves of different frequencies like in a digital telephone line.

  Second,

  (i) I know bandwidth is the range of frequencies that make up a
  composite analog signal. So, if there are two analog signal with
  frequency 5 Hz and 10 Hz making up a composite signal, will it's
  bandwidth be (10-5=) 5 hz?

If you include all frequencies from $0$ Hz upto the max in the signal, then the bandwidth of your composite baseband signal will be $10-0 = 10$ Hz., counting only the positive excursion of frequencies. If you want to consider your signal as bandpass then its bandwidth will be $10-5=5$ Hz.

  (ii) When we say bandwidth of a transmission medium is 5 hertz, do we
  mean to say that it allows any composite signal to pass through it, in
  which difference between the upper frequency and lower frequency
  making up the composite signal is 5 Hz?

No. The actual frequency range of that passband is also important; a medium that can transmit a band in [15-20] Hz cannot transmit a signal in the band [85-90] Hz, even though both have the same bandwidth. You should shift the signal spectrum into the channel band for transmission.

  (iii) Can we define bandwidth of a single sine analog signal?

Bandwidth of $x(t)=\cos(2 \pi f_0 t)$ is its frequency $f_0$ if you treat it like a baseband signal, but otherwise $0$ Hz if you treat it like bandpass signal...

  Third

  (i) Are digital signal nothing but just superimposition of
  theoretically an infinite number of analog signals? If yes, then why
  and how?

digital signal is a misnomer here. It's true that an ideal pulse waveform requires infinite spectrum and infinite number of sine waves to construct it, but a pulse is not the only way to carry digital information as the answer to the first question already indicated...

  (ii) How do we define frequency and bandwidth(in Hertz, not bps) of a
  digital signal?

digital signal is a misnomer here, it's the discrete information that you actually should care, that information is quantified by an entropy measure. Then this will be transmitted via line signalling with analog bandwidth in Hz depending on the efficiency of signalling waveforms.


What the OP shows is the difference between truncation and rounding which in terms of error simply introduces an offset. They generate the same quantization noise spectrally; the standard deviation of the quantization noise (and equivalently noise power) is equivalents since DC offset is ignored in those computations. If there was concern with the offset error introduced (1/2 lsb) it can also easily be subtracted.
Further it should be clarified that DAC's and ADC's will have a finite DC offset regardless due to the analog portion of the electronics, this is often consistent with other DC offsets in the analog chain (consider the offset voltage between the negative and positive Op-amp terminal for example, or carrier leakage in zero-IF receivers), and when an actual issue (for a system where we are trying to measure a DC value) would be incorporated in an overall DC-offset removal algorithm with either calibration or cancellation loop) as part of all offsets combined. Since most DAC's and ADC's have a precision out to what there self-noise factors allow, the ENOB (equivalent number of bits) is often 1 to 2 bits worst than the stated precision (for example, a 12 bit converter may have an ENOB of 10.5 bits). For this reason it is often not even a concern if there is a DC offset introduced of 0.5 bits since that is well below the actual achieved precision overall. These details with regards to ENOB are clarified on ADC and DAC datasheets.
The OP asked for how the noise equivalence can be verified. It is obvious by inspection of the plots given that the only difference between the two is a static offset, but the equivalence can be computed simply enough for further verification by creating a floating point full scale  sine wave and then creating two quantization signals by both rounding and truncating the sine wave to the closest integer. Subtract the quantized fixed-point (integer) signal from the original floating point signal to get the quantization error alone. Then compute the standard deviation of each and they will be equivalent.
This chart I have demonstrates this graphically; specifically what is shown is generalized to a complex waveform (two ADC's with I and Q) but the concept is the same for a single ADC (real waveform), in that the rounded or truncated waveform is the sum of quantization noise with the higher precision waveform (such as floating point or more specifically samples of analog waveform prior quantization through rounding/truncation).

Further, if you use an incommensurate sine wave (choosing a frequency for the sine wave that isn't an integer sub-multiple of the sampling clock), you can also see through subsequent FFT and other statistical analysis how the quantization noise is well modeled as a uniform white noise distribution, and you can also confirm for either case (truncation or rounding) the formula used to estimate the SNR (the total power of the quantization noise floor relative to the full-scale sine wave) as given by:
$$SNR = 6.02 \text{dB/bit} + 1.76 \text{dB}$$
Which I detail further at this link:
What are advantages of having higher sampling rate of a signal?


Why Does the Kalman Filter Remove Only Gaussian Noise?

A Gaussian process is completely specified by its mean and variance.  A Kalman Filter updates the process mean which is the state, and its variance.  These are the sufficient statistics.

The measurement noise is reduced but the process noise  is part of the recursive state history and is tracked.

Your heading question and subsequent paragraph aren't fully consistent.  One can use a linear Kalman Filter when the noise isn't Gaussian in many circumstances but it wouldn' be optimal.  The derivation assumes Gaussian (or Normal) noise and all deterministic inputs are known, as well as knowing initial values.


For such a segmentation, you might want to use graphcuts:
http://en.wikipedia.org/wiki/Graph_cuts_in_computer_vision

However lately there have been better algorithms proposed, such as Branch and Mincut:
http://research.microsoft.com/en-us/downloads/994e80f7-007b-4fed-a6e9-db5f4f32ccd1/

For most of the methods, source codes are available and you could try them out of the box. Also check out:
http://vision.csd.uwo.ca/code/


As with any perceptual model, the answer is not exactly cut and dry.  In my research, I’ve found that many of the hearing models I’ve encountered are based on band pass filter banks.  These banks are numerous and tightly spaced.  The more adjacent banks that are stimulated, the ‘noiseyer’ the signal would sound.

There are examples of inharmonic signal analysis that take noise into account, such as: https://pdfs.semanticscholar.org/b5bf/98425051f47b0685e75da986ee166e41d6ec.pdf.  I’ve never tried to implement such a structure, but as I understand it, you analyze the peaks in the output of a frequency transform to create a ‘deterministic’ signal.  You then subtract the deterministic signal from the original to get the ‘residual’ noise.  You would then perform further analysis on the residual.  I assume that this isn’t exactly what you were looking for in terms of a model, but I figured it may be informative.


Estimate Taylor series coefficients from samples of a function

Instead of exact polynomial fitting, you could use a least-squares fit, which will find the polynomial of specified order that minimizes the total squared error between the fit and the measured $(x_i,y_i)$ pairs. This can help to mitigate the effects of noise on the fit.

Given measurements $y_i$ of a function $y = f(x)$ at domain values $x_i$ ($i = 0, 1, \ldots , N$), choose a polynomial order $M \le N$ (if $M = N$, then you're down to exact polynomial fitting, as the $N$ points uniquely determine an $M$th order polynomial). Then, set up a system of equations that are linear in the desired polynomial coefficients $p_k$:

$$
y_i = p_Mx_i^M + p_{M-1}x_i^{M-1} + \ldots + p_1x_i + p_0, \;\;\;\;\;i = 0, 1, \ldots , N
$$

The least-squares problem can be solved by arranging the measurements into matrix-vector form:

$$
\mathbf{A} =  \left [ {\begin{array}{cc}
 x_0^M & x_0^{M-1} & \cdots & x_0 & 1  \\
 x_1^M & x_1^{M-1} & \cdots & x_1 & 1  \\
 \vdots & \vdots & \vdots & \vdots & \vdots \\
 x_N^M & x_N^{M-1} & \cdots & x_N & 1 \end{array} } \right ],
\;\;\mathbf{y}= \left [ {\begin{array}{cc}
 y_0  \\
 y_1 \\
 \vdots  \\
 y_N  \end{array} } \right ]
$$

The least-squares solution generates the vector of polynomial coefficients $\left[p_M, p_{M-1}, \ldots , p_0\right]$ that minimizes the total squared error in the above linear system. The solution can be computed as:

$$
\mathbf{\tilde p} = (\mathbf{A^TA})^{-1} \mathbf{A^Ty}
$$

It is worth noting that the matrix $(\mathbf{A^TA})^{-1} \mathbf{A^T}$ is also known as the pseudoinverse of the matrix $\mathbf{A}$. You can then use the least-squares polynomial coefficient vector $\mathbf{\tilde p}$ to evaluate the polynomial at any other $x$ values that you wish.


This scilab script implements a simple LMS adaptive filter.

M = 50;
LMS = zeros(M,N);
LMS(:,1) = zeros(M,1);
ERR = zeros(1,N);
y = zeros(1,N);
mu = 0.0005;

for t=M:N,
    Uwindowed = u(t - [0:M-1]');
    y(t) = LMS(:,t)'*Uwindowed.';
    ERR(t) = d(t) - y(t);
    LMS(:,t+1) = LMS(:,t) + mu*Uwindowed.'*ERR(t);
end;

The inputs to it are u(t), the noisy signal and d(t) the desired signal.  For the standard LMS algorithm to work, you need to have an idea of what the noiseless signal looks like.

From your equations:

$\alpha$ appears to be similar to mu (or mu/2).
u is your s+n.
I am not sure what your d is.
Your equation for the error is, er, erroneous. :-)


Pattern recognition using 2D fft in matlab

DFT is a great tool to find periodic patterns in signals.  However, you should understand what kind of period you want to find. If you want to find a series of points that lie on a regular grid (with even spacing), or more mathematically speaking, a subset of points $(x_i,y_i)$, such that there exists $(dx,dy)$ for which the following equations hold:

$x_i = x_1 + dx*i $
   $ y_i =y_1 + dy*j $

(The patterns can also be rotated, and DFT will still find them).

You should compute a 2D histogram of your data, with large number of bins, and run DFT on it. Afterwards, you should look for a peak in any location expect the DC. An ideal noise should have no patters, thus the contrast of the frequencies responses should be very low.


Triangle-kernel filtering problems at non-integer ratios

I went with Jason R's idea of adjusting the linear filtering, except I didn't like the idea of doing a division on every output pixel, using more than just a simple kernel or how normalising that way would ruin the smooth borders. So instead I modified the kernel itself, which it turns out is done by flattening the peak of the triangle.

It seems like a pretty good solution, and you can see how the sum of the triangles progresses into something resembling a Gaussian function.

Here's how it's done, first of all there's a 'knee' joining the flat top with the rest of the triangle. It ping-pongs on the horizontal axis between 0 and 0.5 like so. Note: $n$ is the scale multiplier for the width of the kernel and the size of the downscaling, it's always 1 or lower, it's the reciprocal of the number displayed at the top of the animation. The smaller $n$, the smaller the output image, the wider the kernel over the original image.

$knee = 0.5 - |(n^{-1} \mod 1) - 0.5|$

Next we need to calculate the height of the top of the flattened triangle, which is done like this:

$m = \lfloor n^{-1} + 0.5\rfloor$

$top = \frac{1 - knee \cdot n}{2m - nm^2}$

And finally we turn it into an affine function which we can use directly when calculating weights:

$c_1 = \frac{-top}{n^{-1}-knee}$

$c_0 = \frac{-c_1}{n}$

The above is calculated only once and $knee$, $top$, $c_1$ and $c_0$ are stored and used when computing weights like this:

$$f(x) =
  \begin{cases}
      \hfill top   \hfill & \text{ if $|x| \leq knee$ } \\
      \hfill c_1 x + c_0 \hfill & \text{ if $knee < |x| \leq n^{-1}$ } \\
      \hfill 0   \hfill & \text{ otherwise } \\
  \end{cases}$$


What is usual independence assumptions on adaptive filters

I think there is an error in your referenced independence assumption.
$w(k)$ should be the update part $\Delta w(k)$ i.e the

$w(k+1)=w(k)+\mu \Delta w(k) =w(k)+\mu \frac{x(k)*e(k)}{c+x(k)^Hx(k)}$

The error after convergence is uncorelated and zero mean.

The two other assumptions are correct.

If your filter is fast enough you may use it at nonstationary systems as long as you are able to track the variations.


How to calculate biquad coefficients for the partial (1%, 10%) applying of pre-emphasis effect

So, sox defines its biquad operation parameters as follows:

biquad b0 b1 b2 a0 a1 a2

  Apply a biquad IIR filter with the given coefficients. Where b* and a* are the numerator and denominator coefficients respectively.

  See http://en.wikipedia.org/wiki/Digital_biquad_filter (where a0 = 1).

  This effect supports the −−plot global option.

A biquad filter is just a special case of an IIR filter, which is a linear system!

That means that you can just have two signal streams, one filtered and one un-filtered, and add them up, scaled by two coefficients with constant sum.

The only thing you'd have to care about is delay, which is a bit tricky to calculate for IIRs; now, luckily, the RIAA biquad guarantees phase errors of this filter to be below 30% at any given frequency, so you can probably get away with summing up with a single-sample delayed input signal.

A single sample delay is a FIR filter with coefficients

b0 = 0, b1 = 1

and any FIR is also understandable as IIR with all aN==0 but a0=1, since IIRs are linear, you can directly modify the coefficients of your biquad by mutliplying them with the aforementioned coefficient and adding 1-coefficient to a0 and b1


Non Local Means Filtering Technique Compared to Gaussian Blur and Total Variation Denoising

Gaussian smoothing and TV denoising are local filters that only consider the neighborhood around each target pixel; on the other hand, non-local means updates a pixel's value by using a weighted average of the pixels judged to be most similar.

The reason why it is called non-local is because it does not consider the spatial relationship between pixels, instead, the similarity between different patches and the patch around the target pixel is computed to determine the weight mentioned above (inversely proportional to the distance between its patch and the reference patch of the target pixel).

Yet some studies also showed that Non-Local means is a local image denoising algorithm, in which experimental results indicated that the bias of the non local means estimator is an increasing function of the radius of the similarity searching zone (patch).

The paper (A.Buade et al's On image denoising methods) provided more detailed comparisons between non local means and Gaussian, TV from the aspect of (1) the method noise, (2) the mean square error, and (3) the visual quality of the restored images by simulation. You can skip the theoretical part by reading section 6 for a brief conception of the effect of those methods, which may also be helpful in better understanding the meaning of those mathematical formulas in previous sections.


Remove multiplied white noise from periodic signal

I suppose you can use synchronous averaging in this case due to the signal being periodic and the nature of the noise (as long as it is sufficiently uncorrelated). See MATLAB example:

ts = 1e-4; % sampling time
tend = 10; % end time
t = 0:ts:tend-ts; % time vector

fx = 10; % periodic signal frequency
x = 5*sin(2*pi*fx*t); % periodic signal

switch 2
    case 1
        n = 1*rand(1,length(x)); % noise in (0,1)
    case 2
        R = 50; % upsampling factor
        nn = 1*rand(1,length(x)/R+1);
        tn = 0:ts*R:tend;
        nnn = timeseries(nn,tn);
        nnnn = resample(nnn, t, 'zoh');
        n = nnnn.Data(:).'; % noise that looks like OP's plot
end

%y = x + n; % additive
y = x.*n; % multiplicative

nsp = 1/(fx*ts); % no. samples per period
NP = fx*tend; % no. periods

avg = zeros(1,nsp); % initialise
for k = 1:NP % loop over number of periods available

    Is = (k-1)*nsp+1; % index to start of period
    Ie = k*nsp; % index to end of period

    avg = avg + 1/NP*y(Is:Ie); % add period to average
end
avg = 2.*avg; % scale dependent on noise and how it enters the signal

plot(t(1:nsp),y(1:nsp),t(1:nsp),avg,t(1:nsp),x(1:nsp)) % plot signals
legend('signal with noise', 'average', 'signal')


Does a phase rotation introduce a delay?

I would assume a discrete-time context (though your notation suggests continuous-time).

Then you have a misunderstanding of the relation between the "delay" and "phase". Delay is measured in time-domain whereas the phase rotation is defined in frequency-domain.

So if there was a delay $\pm d$ between $x_1[n]$ and $x_2[n]$, then $y_1[n]$ (defined as the inverse DTFT of $Y_1(\omega)$ defined below) will have zero delay wrt $x_2[n]$ in time:

$$ Y_1(\omega) = e^{ \mp j \omega d } X_1(\omega)$$


Theory behind Linear Predictive Coding (LPC)

First, to say linear predictive coding (LPC) is "more tolerant of transmission or encoding errors" isn't entirely true. The form in which the coefficients is transmitted makes a big difference. For example, if the linear prediction coefficients are solved for, they can be very sensitive to quantization, much like high order IIR filter coefficients (this is because the synthesis filter will be IIR, but more on that later). However, if they are transmitted in some other form, this problem can be mitigated easily.

One way is to transfer the reflection coefficients. If you recursively solve for a k-th order linear prediction filter, the highest order coefficient at each stage is called the reflection coefficient. These can be used together to completely characterize the system (which can easily be seen from the Levinson recursion). In fact, you can use all of them together to form a lattice filter. These filters are often used when quantization is a concern, as they're much more robust to low bit counts. In addition, if the magnitude of these reflection coefficients are bounded by unity, you're guaranteed a BIBO stable filter which is critical for LPC where the filter is used to synthesize your signal. There are other methods such as line spectral pairs which are frequently used, but aren't as intuitively defined as the reflection coefficients in relation to AR modeling for LPC.

Now, to address the first question, the theory of LPC revolves around vocal tract modeling. Essentially, we're modeling speech as air vibrating as an input to a tube of some structure. You can look for some resources that go into much more detail to flesh out this model (length of tubes, intensity of air, structure etc). These resources relate these structures directly to IIR filters responding to various stimuli, white noise for example.

So when we solve for the linear prediction coefficients, we're looking for the coefficients such that if we input our signal (voice for example) into an FIR filter created from the coefficients, we get white noise as an output. So think about what that means. We're inputting a highly correlated signal, and outputting a white noise sequence. So in effect, we're removing all linear dependence of that signal. Another way to look at this, is that all of the meaningful information is contained in the coefficients that remove this linear dependence. Therefore, we can transfer these coefficients (or some form of them as above), and the receiving end can recreate the signal. This is done by inverting the linear predictive FIR filter to create an IIR filter, and inputting white noise. So the compression comes in from removing this linear dependence, and transferring the coefficients. This is why the Burg method is also sometimes referred to as the maximum entropy method, as it aims to maximize the "randomness" or whiteness of the output noise in the linear prediction filter. Another way to look at this, outside of physical modeling, is that vocal spectrum (with exception to unvoiced sounds such as hissing) tend to be very impulsive, making them excellent candidates for AR modeling.

To answer your final question, I'm not sure what you're asking entirely. LPC, or linear predictive coding is meant to "compress" the signal assuming that it can be efficiently modeled as previously discussed. You can certainly use linear prediction to do "shortterm prediction" as you've mentioned. This is the implicit basis behind the high resolution AR methods used for power spectral density estimation. The autocorrelation sequence may be recursively extended from it's finite form from the limited data record to infinity as the theoretical autocorrelation sequence of the unwindowed sequence. This is also why AR methods of PSD estimation don't exhibit sidelobe phenomena.


Lower bound of weighted average of sequence

You can recognize the function $$f: x\mapsto \frac{x}{(1+x)^2}\,.$$
The question can be rephrased as: how do compare the average of the function values $\overline{f(a_i)}$ with the function of the average value $f(\overline{a_i})$? Graphically, this can be illustrated by the following graph (spoiler: from Behold! Jensen's inequality). Here, the average is represented by $\mathbb{E}$, the statistical expectation, akin to the classical average.

The intuition here is the following: for a convex function, as the above, the average of the $y=f(x)$ ($\mathbb{E}[f(x)]$) is higher than the function of the averages of $x$  (or $f(\mathbb{E}[x])$).
In the discrete form, Jensen's inequality can be proven by induction Proof 1 (finite form).
This inequality is reversed for a concave function. The reverse inequality is the one you are interested in. Now, function $f(x)$ is not concave on $[0,\infty[$, graph from RechnerOnline:

With $a_1=0$ and $a_1=1$, we get $$\overline{f(a_i)}-f(\overline{a_i}) = 1/8 - 2/9 \le 0 \,.$$
With $a_1=1$ and $a_1=9$, we get $$\overline{f(a_i)}-f(\overline{a_i}) = 17/100 - 5/36 \ge 0 \,.$$
Hence, the inequality has a sign change. It would be valide for instance on some interval $[\alpha,+\infty)$, such that $f$ is concave.


Assume $N=16$. so the counter $c$ chooses values between $0$ to $15$ and there are $s=1,2,\cdots,\log2(N)$ number of stages.
Based on your comment, the input-output relationship for stage $s$ can be given by:

Represent the counter in binary form in $\log_2 N$ bits (here 4 bits).
Circularly shift the counter $s-1$ times to the right.
Bit reverse the result.
Represent in decimal.

The following Matlab code can give the output in vector form as explained in the comment:

N = 16;
stage = 1              %Do it here for stage #1
c = 0:N-1;
c_bin = de2bi(c);
butterfly_idx = bi2de(fliplr(circshift(c_bin',stage-1)'))

Example for N=16 and stage=2`:

c_bin =

     0     0     0     0
     1     0     0     0
     0     1     0     0
     1     1     0     0
     0     0     1     0
     1     0     1     0
     0     1     1     0
     1     1     1     0
     0     0     0     1
     1     0     0     1
     0     1     0     1
     1     1     0     1
     0     0     1     1
     1     0     1     1
     0     1     1     1
     1     1     1     1

fliplr(circshift(c_bin',1)')

ans =

     0     0     0     0
     0     0     1     0
     0     1     0     0
     0     1     1     0
     1     0     0     0
     1     0     1     0
     1     1     0     0
     1     1     1     0
     0     0     0     1
     0     0     1     1
     0     1     0     1
     0     1     1     1
     1     0     0     1
     1     0     1     1
     1     1     0     1
     1     1     1     1

bi2de(fliplr(circshift(c_bin',1)'))

ans =

     0
     4
     2
     6
     1
     5
     3
     7
     8
    12
    10
    14
     9
    13
    11
    15

which is the expected result for the second stage.


What are Local and Global Inpainting Techniques in Image Processing?

Usually one could divide the methods according to:

Global Method
The function optimized is all over the image with no local data.
Namely it finds best solution involving all pixels.
Diffusion Based is usually iterative method to solve such a Global model.
Local Method
The optimization function is local for a group of patches and holes.
Most "Patch" based methods are local in their nature.

In practice, state of the art methods combine both approaches.
Namely global regulator over a local optimization function or the other way around.


piecewise-quadratic vs piecewise-cubic vs higher order polynomial interpolation?

The mathematically correct way of signal reconstruction would be sinc interpolation https://en.wikipedia.org/wiki/Whittaker%E2%80%93Shannon_interpolation_formula . However that's not very practical so most real world implementation use an approximation. But it's still based on the same principal: it is essentially a low pass filter and NOT a piecewise polynomial interpolation.
That's how most of real world signal reconstruction is done.
Other interpolation methods may be useful in certain special cases for specific applications. To differentiate between these methods, you would have to formulate what application specific metrics distinguish between a "good" and a "bad" interpolation. Then you can analyze different piecewise interpolation techniques with "typical" signals using the defined metrics and see how the performance is impacted and what is optimal.


What's the difference between Channel shortening and time-domain equalizer

Yes,I agreed these terms are very common in communication engineering,

I was asking about channel shortening long time ago, and I got an answer telling me that channel shortening is just calculating the coefficients of the DFE. you can read the details in the below link

Channel Shortening for underwater acoustic channel communication

thanks


What motion estimation techniques that could estimate the motion of every pixel?

You can use Farneback to estimate dense flow.  In OpenCV, the function call is calcOpticalFlowFarneback().

Please note that the initial flow is important and you can use one of OPTFLOW_FARNEBACK_GAUSSIAN or OPTFLOW_USE_INITIAL_FLOW.

It depends on the depth of the pre calculation to estimate the segmentation before calculating the frame to frame flow.


How do I measure the amount of information within an image?

My problem with Shannon's entropy is that it is calculated without considering spatial structures.

That's not true!

Entropy is the expectation of Information, which is the negative logarithm of probability.

So, if you have correlated pixels (for example, in pictures of the first kind, say, a pixel in anything but the top row always has the same value as the pixel above), then that pixel's info drops to zero.

So, your statement is only true if you assume all pixels to be independent, which they aren't.


why do directional derivatives work better with a gaussian filter?

Well, look at your original picture: it's constant for all points but the edges, which means your derivative is zero for all points but these edges.

By applying a "rounding, smoothing" filter to it, you "smear" the edges enough to make the derivative be non-zero for multiple pixels, in every direction.


Because the upper left square would correspond to values of $p=0$ and $q=0$, and $cos(\alpha) = 1 $ if $\alpha=0$, so you get a constant term, that is why it is all gray.
Basically as you increment $p$ and $q$ you get basis images that oscillate more either in vertical or horizontal directions.


How to normalize frequency response charts

First of all, I wouldn't worry too much about the speaker response since it is relatively flat and the microphone has a much bigger roll-off.

Since you've captured the frequency response using sweep, why not to skip the whole part of designing the filter that mimics the frequency response and use the original impulse response? I don't know what kind of software you used for doing your measurements, but a tool such as REW can export the IR wav file. Then it's just a matter of convolving it with your input signal. If your device has some significant amount of self-noise, you can record it in a quiet environment and add it post convolution.

Regarding the decibel scaling, an easy thing would be to play a signal at a known level (for example 94 dB SPL) and record it with the microphone. Then check what is the dBFS level of that recording. This will allow you to map the dB SPL to dB FS. Once it's done, you can convolve your signals and analyse the spectra to see what level they simulate. If any gain adjustment is necessary then you can do it afterwards. This, of course, assumes that you know the original level of those recordings, which is generally not true. If you had the same sort of measurement for the device that was used to capture those recordings, then it's easy. You know what is the sensitivity of that device, you know how much it differs from your microphone, so you can correct for that either by adjusting the level of impulse response or applying gain post convolution. Finally, you can simply ignore all of that and set the gain of your impulse response so the frequency response peaks at 0 dBFS. This means no boosting - only attenuation.


What Are the Best Algorithms for Document Image Thresholding (Example Inside)?

Your image doesn't have uniform brightness,so you shouldn't work with a uniform threshold. You need an adaptive threshold. This can be implemented by preprocessing the image to make the brightness more uniform across the image (code written in Mathematica, you'll have to implement the Matlab version for yourself):

A simple way to make the brightness uniform is to remove the actual text from the image using a closing filter:

white = Closing[src, DiskMatrix[5]]

The filter size should be chosen larger than the font stroke width and smaller than the size of the stains you're trying to remove.

EDIT:
I was asked in the comments to explain what a closing operation does. It's a morphological dilation followed by a morphological erosion. The dilation essentially moves the structuring element at every position in the image, and picks the brightest pixel under the mask, thus :

removing dark structures smaller than the structuring element
shrinking larger dark structures by the size of the structuring element
enlarging bright structures

The erosion operation does the opposite (it picks the darkest pixel under inside the structuring element), so if you apply it on the dilated image:

the dark structures that were removed because they're smaller than the structuring element are still gone
the darker structures that were shrunk are enlarged again to their original size (though their shape will be smoother)
the bright structures are reduced to their original size

So the closing operation removes small dark objects with only minor changes to larger dark objects and bright objects.

Here's an example with different structuring element sizes:

As the size of the structuring element increases, more and more of the characters is removed. At radius=5, all of the characters are removed. If the radius is increased further, the smaller stains are removed, too:

Now you just divide the original image by this "white image" to get an image of (nearly) uniform brightness:

whiteAdjusted = Image[ImageData[src]/ImageData[white]*0.85]

This image can now be binarized with a constant threshold:

Binarize[whiteAdjusted , 0.6]


Convolving Room Impulse Response with a Wav File (python)

The way you read the files is a little odd.  The method used here is clearer with respect to the sampling rate and the signal.
You will note that the above example calls pcm2float after reading.  This may be the cause of your problem.  The wave reader returns an array of ints.  Processing those (quite large) values may well cause problems.  Try converting your wave data to floats before doing the convolution.


Nyquist Rate of cosine modulated function

Update: I feel there is no need to bring concept of Bandpass sampling (as mentioned in comments below) because this is specifically regarding a simple misunderstanding that OP has. The text referenced asks about Nyquist Rate (which wrongly mentions it as Nyquist Frequency) but OP asks about highest Frequency in $y(t)$.

For the component $\frac{1}{2}X(\Omega-\Omega_0)$ the highest frequency content is at $\Omega_0+\Omega_N/2$. The lowest (starting frequency) is at $\Omega_0-\Omega_N/2$. Similary for the other component $\frac{1}{2}X(\Omega+\Omega_0)$, the starting frequency is at $-\Omega_0-\Omega_N/2$ and ending frequency is at $-\Omega_0+\Omega_N/2$.

So if you see the overall frequency content of $y(t)$, lowest = $-\Omega_0-\Omega_N/2$ and highest = $\Omega_0+\Omega_N/2$. This is your signal $y(t)$ whose highest frequency is $\Omega_0+\Omega_N/2$.

The Nyquist rate for this is $2\Omega_0+\Omega_N$ if you consider $y(t)$ as baseband signal in the picture below. But the minimum sampling rate will be much lesser than this if you consider it as Passband signal (concept of bandpass sampling as mentioned in another answer)

  If the Highest frequency of $X(\Omega)$ is $\Omega_N$

No. The highest frequency content in this case is $\Omega_N/2$ which is half the Nyquist rate. Otherwise this would cause aliasing.


Q matrix and updating times in a Kalman filter

Kalman filter works in a predictor-corrector or predict-update sequence. When no new measurement is available the only thing that you can do is to predict(NOT update) the state with your assumed dynamic model. The prediction is done not only on the state vector itself but also on its covariance matrix ${\bf{P}}$ (Not the measurement covariance matrix $\bf{R}$). A popular model commonly used for model regularization (AKA model noise matrix, or the matrix $\bf{Q}$) is the Singer model from 1990 (see [1]).

In your case the transition matrix is given by $\Phi_k =  \begin{bmatrix}1 & dt \\0 & 1 \end{bmatrix}$ , so the state covariance matrix is predicted using
${\bf{P}}_{k|k-1}=\Phi_k {\bf{P}}_{k-1|k-1} \Phi_k^T + {\bf{Q}}_k$. That way the state matrix ${\bf{P}}$ "inflates" to account for the uncertainty.

[1] Singer, R. A., “Estimating Optimal Tracking Filter Performance for Manned Maneuvering Targets,” IEEE Transactions on Aerospace and Electronic Systems, vol. AES-6, no. 4, .pp. 473-483, July 1990.


Deriving the Resolution Equation of an FMCW RADAR

Definitions:

$\tau$: round-trip delay.  $\Delta \tau$ = delay resolution.
$r$: range.  $r=c\tau/2$.  Range resolution $\Delta r = c \Delta \tau / 2$.
$B$: chirp bandwidth
$T$: chirp sweep time
$R$: chirp rate.  $R=B/T$

Basic steps:

We estimate the delay by estimating the frequency shift between the received and transmitted signal.  Typically the resolution of a frequency estimate is the inverse of the observation time, which in this case is $1/T$.  Why?  Think of taking the Fourier Transform of a tone with a rectangular window T.  The window will have nulls at $f-1/T$ and $f+1/T$, where $f$ is the tone's frequency.
Because the chirp rate is $R=B/T$, the time resolution corresponding to a frequency resolution of $\frac{1}{T}$ is $\frac{1/T}{B/T} = \frac{1}{B}$.
Plugging this into the range resolution equation, we have $\Delta r = c/(2B)$.

Intuitively, we improve range resolution by (1) improving resolution of the frequency shift estimate (i.e., increasing sweep time), and (2) decreasing the time delay corresponding to a given frequency shift (i.e., increasing chirp rate).  To do one without sacrificing the other, we need to increase the bandwidth.  If we double the chirp rate but cut the sweep time in half, thus keeping the bandwidth the same, the time delay per frequency shift decreases but the frequency resolution gets worse.  Likewise, if we cut the chirp rate in half but double the sweep time, the time delay per frequency shift increases while the frequency resolution improves.  In either case, the range resolution does not change.


Why is this matrix invertible in the Kalman gain?

Note that $\mathbf{P} _{k\mid k-1}$, just like $\mathbf{R}_k$, is also a covariance matrix, and for this reason it is (at least) positve semi-definite, i.e., $\mathbf{y}^T\mathbf{P}_{k\mid k-1}\mathbf{y}\ge 0$ for $\mathbf{y}\neq\mathbf{0}$. Now set $\mathbf{y}=\mathbf{H}_k^T\mathbf{x}$ to see that also $\mathbf {H} _{k}\mathbf {P} _{k\mid k-1}\mathbf {H} _{k}^{\text{T}}$ is at least positive semi-definite (positive definite if $\mathbf{P}_{k\mid k-1}$ is positive definite and $\mathbf{H}_k$ has full rank). Finally, note that the sum of two positive semi-definite matrices is positive semi-definite.

For invertibility, we require that the sum of the two matrices is positive definite. This is the case if at least one of the two matrices is positive definite. In practice we can rather safely assume that both $\mathbf{R_k}$ and $\mathbf{P} _{k\mid k-1}$ are positive definite. However, $\mathbf {H} _{k}\mathbf {P} _{k\mid k-1}\mathbf {H} _{k}^{\text{T}}$ is usually only positive semi-definite due to rank-deficiency of $\mathbf{H}_k$, so for invertibility of the sum of the two matrices we have to rely on the positive definiteness of $\mathbf{R}_k$.


Spectral Entropy describes the complexity of a system. It is defined as follows:

Calculate the spectrum $X(\omega_i)$ of your signal.
Calculate the Power Spectral Density of your signal via squaring its amplitude and normalizing by the number of bins.

$$P(\omega_i)=\dfrac{1}{N}\left|X(\omega_i) \right|^2 $$

Normalize the calculated PSD so that it can be viewed as a Probability Density Function (integral is equal to 1).

$$p_i=\dfrac{P(\omega_i)}{\sum_iP(\omega_i)} $$

The Power Spectral entropy can be now calculated using a standard formula for an entropy calculation.

$$PSE = -\sum_{i=1}^np_i\ln p_i $$

In case of boosting of your noise signal, without performing any other processing, the Entropy will change. I guess there is no other way around that.


First part for an answer. There are several unknowns in you question. The fist is of course:

why do you really want to bring the signal to the same level?

and the second:

what is this "same level"?

Assuming that the baseline is the lower  levels, do those simulated signals approximately match your expectations?


How to apply the CA-CFAR methods in practice?

Constant false alarm rate detection is more of a general technique that underlies a class of approaches; you can't really say that there's "an" implementation of CFAR processing. CFAR techniques are useful when you need to detect the presence of a signal where there is some uncertainty as to what the magnitude of the signal will be, especially relative to any background noise or interference. The general principle of constant false alarm rate detection is that the detector makes some effort to estimate the background level and set a detection threshold adaptively relative to the observed ambient environment. This approach makes the performance of the detector (typically measured quantitatively using the probability of detection $P_d$ and probability of false alarm $P_{fa}$) insensitive to the bulk magnitude of the signal under observation; only relative changes in the statistics of the signal over time (such as a transition from noise only to signal-of-interest plus noise) are important.


What does the phase of a FFT bin of an audio signal represent?

The frequency bin with the highest magnitude is very often not the fundamental pitch frequency, especially for lower (male) voices and larger stringed instruments.  It could instead represent a strong overtone or harmonic.

The FFT phase result is very important in representing any transient waveform shape or timing, for instance whether a note onset or percussive hit occurred toward the beginning or end of the FFT frame.  It's thus not very useful in matching sounds unless the FFT frame timing is somehow very precisely matched to the sound's timing.

For a very narrow-band signals that are stationary across multiple FFT frames with known time offsets, the FFT phase results can also potentially be used for frequency estimation finer than the FFT bin spacing resolution.


You'll have to convert back into linear units before you can evaluate the SINR formula.

${\textrm{SNR}_{linear}=10^{\big(\frac{\textrm{SNR}_{dB}}{10}\big)}}$ and likewise for CIR...

${\textrm{SINR}_{dB}  = 10 \log_{10}(\frac{1}{\frac{1}{\textrm{SNR}}+\frac{1}{CIR}}})$

${BER = \frac{1}{4{\textrm{SINR}}_{linear}}}$


Computing the Power Spectrum Density (PSD) on a CSV File in Python

I have not used the mlab functions in matplotlib, but from what I can discern from your time plot, I see a pattern that repeats roughly every 15,000 samples, at 1.25 Gsps according to your code.

This means we should see a strong spectral component at $1.25e9/15000=83.3$ KHz.

The reason we do not see this, I suspect is the relatively short FFT length (or very high sampling rate compared to the spectral occupancy, or both!).

Let me explain: The FFT length is currently 256 samples. At 1.25GHz this is only $256/1.25e9 = 204.8$ ns.

The frequency resolution prior to using any additional window (meaning using a rectangular window) is $1/T$ where T is the length of the FFT in time (for further detail on that see:
What happens when N increases in N-point DFT )

So in your case the frequency resolution is greater than $1/204.8e-9$ = 4.88 MHz!!  Any window will reduce sidelobes but at the expense of increasing the mainlobe, so the width of your mainlobe given the Hanning window used will cause the resolution to be even wider than 4.88 MHz. (But definitely good to include windowing, otherwise the additional sidelobes would be much higher). You can see this post for additional details on the windowing choices:Why would one use a Hann or Bartlett window?  (From this we see the Hanning window is approximately twice as wide as the rectangular window, so your example would have approximately 10 MHz of resolution bandwidth per bin in the FFT)

Looking at your plot, it looks like we see exactly this, which is the FFT bin of your strong DC component, and this single FFT bin is approximately 10MHz wide in bandwidth! So that is burying any signals of interest that you want to see.

My recommendation is to increase the FFT length within what your processing will allow. And if that starts to get problematic then reduce the sampling rate to get further resolution (however there are challenges in doing that properly that I won't get into except to say that you must be sure higher frequency components are filtered out prior to resampling). Without requiring a sampling rate change, you would want to use an FFT length of at least $1.25e9/10e3$ = 125,000 to get 20 KHz of resolution BW while using the Hanning window. That's a reasonable length and I suspect you can go much higher than that for finer resolution.

Also as a side note I don't see where the impedance is set, so I guess it is assumed to be normalized to 1 ohm. Just mentioning in case you cared about knowing the actual power levels.


It's about units and dimensions. And I'm probably not the guy to cast an answer on that but I would like to state it as follows. In short; length, time, mass etc. are dimensions, whereas meter, seconds, kg are units of them respectively in SI. (for other dimensions and units, a physics book might help.)

Now, if a wheel rotates 7 times in 2 seconds; it's rotation frequency $f$ is found as:
$$ f = \frac{7 ~\text{ cycles}}{2 ~\text{ seconds}} = 3.5 ~\text{cycles-per-second (cps)}.$$

As can be seen, frequency is given by a (cycle) number divided by time. The cycle number has no physical units or dimensions; it's just a repetition count. Whereas the duration time has a unit of seconds (in SI).

Thus the dimesion of frequency is 1/time, whose SI unit is $1/s = s^{-1}$. This unit is also, and more frequently, denoted as Hz where it (implicitly) stands for (cycles) per second indeed. So, $3.5$ Hz stands for $3.5$ cycles per second.

Now, similarly the period is defined as the duration time per one cycle. Again note that since cycle is a number and does not posses a unit, then the dimension of the period is time and its unit is seconds in SI. But you can always denote it also as seconds per cycle, rarely practiced though.

Finally, using both notations you have
$$ T = \frac{1}{f} $$

So if frequency $f$ is defined with a unit of Hz ($s^{-1}$), then the period will have a unit of seconds (s) in SI. Otherwise if you define frequency as cycles per second then the period unit will be seconds per cycle.

So if you define the frequency as $3.5$ cycles-per-second , then the period associated is:
$$T = \frac{1}{f} = \frac{1}{ 3.5 \frac{\text{cyc}}{s} } = \frac{2}{7} \frac{\text{s}}{\text{cyc}}$$

Otherwise, in accepted norm, as was you would define freqquency as $3.5$ Hz. its period in SI will be:
$$ \boxed{ T = \frac{1}{f} = \frac{1}{ 3.5 s^{-1} } = \frac{2}{7} ~\text{ s} } $$


How to implement a FIR filter with correction of group delay on c#?

Any kind of digital filter will cause the the output signal to be delayed by some amount of samples.

From what I gather, you are trying to run a signal through a high pass filter (is it an FIR or IIR?) and correct the group delay by "filtering the first time, inverting the response in time...". I personally have never been taught or have read of such an algorithm to correct for group delay.

One thing that I am confused about is why you want to fix the group delay. You are not ok with the input signal being delayed by a few hundred samples? All you do is have to wait and all will be fine.

The only reason why I would think you'd need to fix the group delay, is to compare the filtered response to the original signal in the time domain. Or do something much more advanced than garden variety filtering. If you are trying to compare time domain responses, you will have to delay the reference signal by the same amount of FIR filter delay your high pass filter has.

If this answers your question, respond and I'll go into more detail on how to delay your reference signal, and how to find the sample delay of an FIR filter. You can also google "sample delay of fir filter".


Question about eigendecomposition, signal subspace and their properties

My first swing at the answer had some very incorrect claims.

I do not have access to the article, so I am inferring some things from the portion posted in the question.

NOTA BENE: My arguments assume that the eigenvectors of $\mathbf{R}$ are arranged so that the first $n$ belong to the signal subspace and that the last $m-n$ belong to the noise subspace. That is how the formulas appear so clean.

Since $\mathbf{R}$ is a covariance matrix, it is a positive-definite Hermitian matrix. This ensures that

all of its eigenvalues are positive,
it is unitarily diagonalizable.

That last property means that there is a unitary matrix $\mathbf{U}$, whose columns are orthonormal eigenvectors of $\mathbf{R}$, such that $\mathbf{U}^{\dagger}\mathbf{R}\mathbf{U}$ is equal to a diagonal matrix that has $\mathbf{R}$'s eigenvalues on its diagonal:
\begin{eqnarray}
\mathbf{U}^{\dagger}\mathbf{R}\mathbf{U} = \textrm{diag}(\lambda_1,\ldots,\lambda_n,\underbrace{\sigma^2,\ldots,\sigma^2}_{\textrm{$m-n$ $\sigma^2$s}}).
\end{eqnarray}

Because the article mentions $\mathbf{E}_s^{\dagger}$ and $\mathbf{E}_n^{\dagger}$ and not $\mathbf{E}_s^{-1}$ and $\mathbf{E}_n^{-1}$ in the eigenvalue-eigenvector decomposition of $\mathbf{R}$, I suspect that the unitary diagonzaling matrix of $\mathbf{R}$ is $\mathbf{E}_s + \mathbf{E}_n$, where $\mathbf{E}_s^{\dagger}\mathbf{E}_n = \mathbf{O}$, an all-zero matrix. I will demonstrate using the decomposition given in the article.
\begin{eqnarray}
&& (\mathbf{E}_s^{\dagger} + \mathbf{E}_n^{\dagger})\color{red}{\mathbf{R}}(\mathbf{E}_s + \mathbf{E}_n)\\
&=& (\mathbf{E}_s^{\dagger} + \mathbf{E}_n^{\dagger})\color{red}{(\mathbf{E}_s\mathbf{\Lambda}_s\mathbf{E}_s^{\dagger} + \sigma^2\mathbf{E}_n\mathbf{E}_n^{\dagger})}(\mathbf{E}_s + \mathbf{E}_n)\\
&=& (\mathbf{E}_s^{\dagger}\mathbf{E}_s\mathbf{\Lambda}_s\mathbf{E}_s^{\dagger} + \sigma^2\mathbf{E}_n^{\dagger}\mathbf{E}_n\mathbf{E}_n^{\dagger})(\mathbf{E}_s + \mathbf{E}_n),
\end{eqnarray}
where I have already dropped $\mathbf{E}_s^{\dagger}\mathbf{E}_n$ and $\mathbf{E}_n^{\dagger}\mathbf{E}_s$, both of which are all-zero matrices. The next step is
\begin{equation}
\mathbf{E}_s^{\dagger}\mathbf{E}_s\mathbf{\Lambda}_s\mathbf{E}_s^{\dagger}\mathbf{E}_s + \sigma^2\mathbf{E}_n^{\dagger}\mathbf{E}_n\mathbf{E}_n^{\dagger}\mathbf{E}_n.
\end{equation}

The OP states that $\mathbf{E}_s$ is $m\times n$ and has rank $n$, where $n < m$. But I suspect (hope) that $\mathbf{E}_s$ is $m\times m$ with rank $n$ and that and $\mathbf{E}_n$ is $m\times m$ with rank $m-n$, with the following structures:
\begin{eqnarray}
\mathbf{E}_s &=& (\mathbf{u}_1\cdots\mathbf{u}_n\underbrace{\mathbf{0}\cdots\mathbf{0}}_{m-n}),\\
\mathbf{E}_n &=& (\underbrace{\mathbf{0}\cdots\mathbf{0}}_{n}\mathbf{u}_{n+1}\cdots\mathbf{u}_m)
\end{eqnarray}
That is, $\mathbf{E}_s$ should have $n$ orthonormal columns followed by $m-n$ all-zero columns, while $\mathbf{E}_n$ starts with $n$ all-zero columns and ends with $m-n$ orthonormal columns.

If these conditions were true, then we would have
\begin{eqnarray}
\mathbf{E}_s^{\dagger}\mathbf{E}_s &=& \textrm{diag}(\underbrace{1,\ldots,1}_{n},\underbrace{0,\ldots,0}_{m-n}),\\
\mathbf{E}_n^{\dagger}\mathbf{E}_n &=& \textrm{diag}(\underbrace{0,\ldots,0}_{n},\underbrace{1,\ldots,1}_{m-n})
\end{eqnarray}
and
\begin{eqnarray}
&& (\mathbf{E}_s^{\dagger} + \mathbf{E}_n^{\dagger})\color{red}{\mathbf{R}}(\mathbf{E}_s + \mathbf{E}_n)\\
&=& \mathbf{E}_s^{\dagger}\mathbf{E}_s\mathbf{\Lambda}_s\mathbf{E}_s^{\dagger}\mathbf{E}_s + \sigma^2\mathbf{E}_n^{\dagger}\mathbf{E}_n\mathbf{E}_n^{\dagger}\mathbf{E}_n\\
&=& \textrm{diag}(\lambda_1,\ldots,\lambda_n,\underbrace{\sigma^2,\ldots,\sigma^2}_{m-n})
\end{eqnarray}

Even if $\mathbf{E}_s$ is $m\times n$ and $\mathbf{E}_n$ is $m\times(m-n)$, we still have many of the same properties as long as
\begin{eqnarray}
\mathbf{E}_s &=& (\mathbf{u}_1\cdots\mathbf{u}_n),\\
\mathbf{E}_n &=& (\mathbf{u}_{n+1}\cdots\mathbf{u}_m).
\end{eqnarray}
In particular, we have
\begin{equation}
\mathbf{E}_s^{\dagger}\mathbf{E}_n = \mathbf{O}_{n\times(m-n)}
\end{equation}
and
\begin{equation}
\mathbf{E}_s\mathbf{\Lambda}_s\mathbf{E}_s^{\dagger} + \sigma^2\mathbf{E}_n\mathbf{E}_n^{\dagger}
=
\sum_{k=1}^{n}\lambda_k\mathbf{u}_k\mathbf{u}_k^{\dagger}
+
\sum_{k=n+1}^{m}\sigma^2\mathbf{u}_k\mathbf{u}_k^{\dagger},
\end{equation}
which is a spectral decomposition of a Hermitian matrix such as $\mathbf{R}$.


pitch detection and tracking.  use autocorrelation or perhaps AMDF (average-magnitude difference function) or ASDF (average-squared difference function).
either the "heterodyne oscillator" or, i would recommend, wavetable synthesis techniques to get the amplitude and phase of each harmonic.


Can a unit delay in discrete time be represented by exponential functions?

In continuous time, we can, though it is a mathematical curiosity of little practical use. Suppose we have a smooth continuous-time signal $x$. Using the Taylor expansion,

$$x (t-1) = x(t) - \dot x (t) + \frac 12 \ddot x (t) - \cdots = \left( 1 - \frac{\mathrm d}{\mathrm d t} + \frac 12 \left(\frac{\mathrm d}{\mathrm d t}\right)^2 - \cdots \right) x(t) = \color{blue}{\exp \left( - \frac{\mathrm d}{\mathrm d t} \right) \, x (t)}$$


Real and Imaginary Frequency Responses of a single complex pole

If I understood you correctly you want a formula for the magnitude of the frequency response of a filter with a single complex pole. As you've already noticed, the system function is

$$H(z)=\frac{1}{1-az^{-1}}\tag{1}$$

which is valid regardless of whether $a$ is real-valued or complex-valued. Let's assume that $a=a_r+ja_i$ is a complex-valued number. If we also assume that $|a|<1$ (i.e. the system is causal and stable), then we get the systems frequency response by replacing $z$ by $e^{j\omega}$ in (1):

$$H(e^{j\omega})=\frac{1}{1-ae^{-j\omega}}\tag{2}$$

Now you just have to compute the magnitude of (2). Using

$$|1-u|^2=1-2\text{Re}\{u\}+|u|^2$$

for any complex $u$, we get with $u=ae^{-j\omega}$ from (2)

$$|H(e^{j\omega})|=\frac{1}{\sqrt{1-2\text{Re}\{ae^{-j\omega}\}+|a|^2}}=
\frac{1}{\sqrt{1-2(a_r\cos\omega+a_i\sin\omega)+|a|^2}}\tag{3}$$

EDIT:

This is how you can derive the transfer functions relating the real-valued input to the real part and the imaginary part of the output, respectively. Let $y(n)=y_r(n)+jy_i(n)$. Then you have the following two difference equations for the two real-valued output signals:

$$y_r(n)=x(n)+a_ry_r(n-1)-a_iy_i(n-1)\\
y_i(n)=a_ry_i(n-1)+a_iy_r(n-1)$$

Transforming these equations to the $\mathcal{Z}$-domain gives

$$Y_r(z)=X(z)+a_rz^{-1}Y_r(z)-a_iz^{-1}Y_i(z)\\
Y_i(z)=a_rz^{-1}Y_i(z)+a_iz^{-1}Y_r(z)$$

Now you have two equations for the two unknowns $Y_r(z)$ and $Y_i(z)$. After some algebra you arrive at

$$Y_r(z)=H_r(z)X(z)=\frac{1-a_rz^{-1}}{1-2a_rz^{-1}+(a_r^2+a_i^2)z^{-2}}X(z)\\
Y_i(z)=H_i(z)X(z)=\frac{a_iz^{-1}}{1-2a_rz^{-1}+(a_r^2+a_i^2)z^{-2}}X(z)\tag{4}$$

Note that even though

$$H(e^{j\omega})=H_r(e^{j\omega})+jH_i(e^{j\omega})$$

$H_r(e^{j\omega})$ and $H_i(e^{j\omega})$ are not the real and imaginary parts of the total frequency response $H(e^{j\omega})$ of the complex-valued system because they are both themselves complex-valued. You can obtain the magnitudes of the frequency responses by computing the magnitudes of $H_r(e^{j\omega})$ and $H_i(e^{j\omega})$ as given in (4), i.e. by replacing $z$ by $e^{j\omega}$ and computing the magnitudes of the two complex functions.

Another option is to express the frequency responses $H_r(e^{j\omega})$ and $H_i(e^{j\omega})$ using the total frequency response $H(e^{j\omega})$ (2):

$$H_r(e^{j\omega})=\frac12[H(e^{j\omega})+H^*(e^{-j\omega})]\\
H_i(e^{j\omega})=\frac{1}{2j}[H(e^{j\omega})-H^*(e^{-j\omega})]\tag{5}$$

From (5), the magnitudes can be written as

$$|H_r(e^{j\omega})|=\frac{1}{2}\sqrt{|H(e^{j\omega})|^2+2\text{Re}\{H(e^{j\omega})H(e^{-j\omega})\}+|H(e^{-j\omega})|^2\}}\\
|H_i(e^{j\omega})|=\frac{1}{2}\sqrt{|H(e^{j\omega})|^2-2\text{Re}\{H(e^{j\omega})H(e^{-j\omega})\}+|H(e^{-j\omega})|^2\}}$$

where $|H(e^{j\omega})|$ is given by (3).

EDIT 2:

Here's some Matlab/Octave code using equation (4) to plot the different transfer functions. $|H_r(e^{j\omega})|$ and $|H_i(e^{j\omega})|$ are of course symmetrical because they correspond to real-valued systems, whereas $|H(e^{j\omega})|$ is not due to the complex pole.

w = linspace(0,2*pi,300);          % frequency vector
ar = .4; ai = -.3; a = ar+j*ai;    % complex pole
u = exp(-j*w);                     % u = z^(-1)
Hr = (1-ar*u)./(1-2*ar*u+(ar^2+ai^2)*u.^2);
Hi = (ai*u)./(1-2*ar*u+(ar^2+ai^2)*u.^2);
H = 1./(1-a*u);
plot(w/pi,abs(Hr),w/pi,abs(Hi),w/pi,abs(H))
legend('Hr', 'Hi', 'H');

EDIT 3:

If you need to split the numerators and denominators of (4) into real and imaginary parts, here's what you get:

$$H_r(e^{j\omega})=\frac{1-a_r\cos\omega+ja_r\sin\omega}{\left[1-2a_r\cos\omega+(a_r^2+a_i^2)\cos(2\omega)\right]+j\left[2a_r\sin\omega-(a_r^2+a_i^2)\sin(2\omega)\right]}\\
H_i(e^{j\omega})=\frac{a_i\cos\omega-ja_i\sin\omega}{\left[1-2a_r\cos\omega+(a_r^2+a_i^2)\cos(2\omega)\right]+j\left[2a_r\sin\omega-(a_r^2+a_i^2)\sin(2\omega)\right]}$$


Computation of the Inverse DCT (IDCT) Using DCT Or IFFT

Have a look at Fast DCT Algorithm (PDF Version).

It has both DCT and Inverse DCT using DFT (FFT).
They show how to do a DCT and Inverse without the reflection trick.

The standard (Less efficient then above) way doing so is:

numElements = 10;

vX = randn(1, numElements);

disp(vX); %<! Diusplay the input signal

vDctRef = dct(vX);

% Forward DCR using FFT
vXX     = [fliplr(vX), vX]; %<! Mirroring
vXDft   = fft(vXX);

vGrid = [0:((2 * numElements) - 1)];

vShiftGrid = exp(-1j * 2 * pi * (numElements - 0.5) * vGrid / (2 * numElements));

vXDft2 = real(vXDft ./ vShiftGrid);

vDct    = vXDft2(1:numElements) / sqrt(2 * numElements);
vDct(1) = vDct(1) / sqrt(2);

disp(vDctRef)
disp(vDct)

% Inverse DCT Using FFT
vDct(1) = vDct(1) * sqrt(2);
vDct    = vDct * sqrt(2 * numElements);

vXDft = [vDct, 0, -fliplr(vDct(2:numElements))] .* vShiftGrid;
vXX = ifft(vXDft);

vX = real(fliplr(vXX(1:numElements)));
disp(vX);

A reference code for the IDCT by the more efficient method (As in reference):

numElements = 10;

% Signal (DCT Coefficients)
vXDct = randn(1, numElements);

% Reference Inverse IDCT
vXRef = idct(vXDct);

% Inverse IDCT Using FFT
vGrid = [0:(numElements - 1)];

vShiftGrid      = exp((1j * pi * vGrid) / (2 * numElements));
vShiftGrid      = vShiftGrid * sqrt(2 * numElements);
vShiftGrid(1)   = vShiftGrid(1) / sqrt(2);

vTmp = vShiftGrid .* vXDct ;
vTmp = real(ifft(vTmp));

vX = zeros(1, numElements);

for ii = 0:((numElements / 2) - 1)
    vX((2 * ii) + 1) = vTmp(ii + 1) ;
    vX((2 * ii) + 2) = vTmp(numElements - ii) ;
end

disp(vXRef);
disp(vX);

Another reference is given by Lecture Notes for the Course MAT-INF2360 - Fourier Theory, Wavelet Analysis and Non Linear Optimization.
Have a look at 4.2 Efficient implementations of the DCT at page 115.

Enjoy...

Remark
When it is written DCT above it refers to DCT Type II.


It looks like that is a color camera with the Bayer pattern, and the color camera data is directly used as grayscale data, as if it was coming from a grayscale camera. If it is a color camera, the color camera data should be de-Bayerized first (and smoothed etc processing) to get color RGB data pixels as usual, and then the RGB data can be converted to grayscale if grayscale is needed.


Cross Domain Equivalent to Nyquist Sampling Theorem?

I'd say that this is not only "similar to a cross-domain equivalent to Nyquist's Sampling Theorem", but it simply is the sampling theorem. The sampling theorem does not specify the domains of the signals involved; it is rather a mathematical condition that a function of a continuous variable needs to satisfy such that it is perfectly represented by equidistant samples. It is irrelevant if the independent variable of that continuous function is time, frequency, space or anything else.


First I would recommend filling in the contour of the toy - in case it looks like the one in the second image. You could do this by analyzing the hierarchy output from  findContours: make  white all regions having a parent or by using an iterative morphological operations (not directly implemented in OpenCV).

Once the toy is nice and fat (not just the edge), you can remove the cord. I would use here some erosion with a flat kernel (horizontal line).

See the code below (this is in Python, but you'll have a clear overview):

import cv2

lemBGR = cv2.imread("lem.png")
lem = cv2.cvtColor(lemBGR,cv2.COLOR_BGR2GRAY)

# Dilate the image in order to close any external contour of the leming
kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(5,5))
lem = cv2.dilate(lem,kernel)

# Identify holes in the leming contour
# This could be done by iterative morphological operations,
# but this is not directly implemented in OpenCV
contour,hier = cv2.findContours(lem.copy(),cv2.RETR_CCOMP,cv2.CHAIN_APPROX_SIMPLE)
# And fill them
for c,h in zip(contour, hier[0]):
    if h[3]!=-1:
        cv2.drawContours(lem,[c],0,255,-1)

# Now bring the leming back to its original size
lem = cv2.erode(lem,kernel)

# Remove the cord by wiping-out all vertical lines
kernel = cv2.getStructuringElement(cv2.MORPH_RECT,(15,1))
#lem = cv2.erode(lem,kernel) # first wipe-out
#lem = cv2.dilate(lem,kernel) # then bring back to original size
# erode and then dilate is the same as opening
lem = cv2.morphologyEx(lem,cv2.MORPH_OPEN,kernel)

# Find the contour of the leming
contour,_ = cv2.findContours(lem.copy(),cv2.RETR_EXTERNAL,cv2.CHAIN_APPROX_SIMPLE)

# And draw it on the original image
for c in contour:
    # enter your filtering here
    x,y,w,h = cv2.boundingRect(c)
    cv2.rectangle(lemBGR,(x,y),(x+w,y+h),(0,255,0),2)

# Display the result
cv2.imshow("lem",lemBGR)
cv2.waitKey()

cv2.imwrite("lem-res.png",lemBGR)

And below are my results.


What are the possible ways to find spots with most concentrated contrast on an image?

Pass a sliding window of eg 100x100 pixels over the image stepped 1 pixel at a time in each dimension. Note the maximum and the minimum pixel value. The local contrast is now the local_max/local_min.
A more robust method might be to use a tapered window (eg Gaussian) and a sliding histogram, using eg the 10-percentile and 90-percentile to avoid salt and pepper noise corrupting results.
If (as stated in your update) you are interested in contrast of features larger than a pixel, downsample your image first.
-k


Overlap-windowing vs preserved information

Depending what do you mean by preserve $100\%$ of data stream information. Usually $50\%$ of overlap is enough and you are not loosing too much information. Amount of data within your observation window is definitely sufficient for further analysis.

Although there are some windows, such as Hanning or modified Gaussian window which that are complementary. This allows you to perfectly reconstruct original signal from windowed segments.


The technique demonstrated in the video appears to be "image morphing" which is not specific to 3D.  They have "morphed" a left-eye view to a right-eye view and back.  The resulting animation is a form of "wiggle" stereoscopic presentation.

See http://en.wikipedia.org/wiki/Morphing

Here's an example I just made from an old stereo pair
of a Civil War cannon, from the Library of Congress.  I used
ImageMagick's "-morph" option to rock back and forth between
the two images:

http://www.simplesystems.org:/users/glennrp/stereo/loc_cannon.gif

See also the "wiggle method" of stereoscopy, with some more examples,
at http://en.wikipedia.org/wiki/Wiggle_stereoscopy
This method typically just alternates between the left and right images, but sometimes "morphing" or filling in with intermediate photos are used.


Force an unstable IIR filter to be stable by forcing a nonlinearity in the digital block diagram

You won't be able to avoid the instability with a linear model. With your current discrete model, the filter will blow up at a certain resonance setting. Using a better filter topology or higher precision will allow you to go a bit farther, but it will blow up.

Preventing that would require a saturating nonlinearity in the feedback path. That path is hidden in your implementation. To access it, you need to go back to the 4-section model without feedback and add the nonlinearity when you add the feedback path. But then you won't be able to simplify the filter to a fourth order linear filter.

Your options are then to delay the feedback path with its nonlinearity by one sample and feed it back into the 4-section model input. That will create significant deviations from the intended filter behaviour for high frequencies and high resonances. You can compensate for that somewhat, but it is tedious to get right. See Antti Huovilainen's paper for an approach like that.

The other option is to solve the implicit feedback equation using an iterative method for every single sample. That is doable, but requires significantly more computational effort and some experience with numerical methods to guarantee convergence and stability.

If you go that latter route, you might as well go with a proper nonlinear filter model, because that will need very similar methods to solve.


system function $H(\omega)$ relationship to odd and even components of h[n]

The DTFT relationships

$$x_{even}[n]=\frac12\left(x[n]+x^*[-n]\right)\Longleftrightarrow\textrm{Re}\left\{X(e^{j\omega})\right\}$$

and

$$x_{odd}[n]=\frac12\left(x[n]-x^*[-n]\right)\Longleftrightarrow j\,\textrm{Im}\left\{X(e^{j\omega})\right\}$$

hold for any sequence $x[n]$ for which the DTFT exists. There is no assumption about $x[n]$ being real-valued or causal (note the complex conjugation $^*$ in the definition of even and odd signals). If $x[n]$ is real-valued you can leave out the conjugation.

Note that the DTFT of the odd part $x_{odd}[n]$ equals $j$ times the imaginary part of the DTFT $X(e^{j\omega})$, so you have

$$X(e^{j\omega})=\textrm{DTFT}\{x_{even}[n]\}+\textrm{DTFT}\{x_{odd}[n]\}$$

(without a $j$ on the right-hand side).


What is plotted against the DFT of a windowed sampled signal?

As you've noted, for a length-$N$ time-domain sequence $x[n]$, the $\mathrm{DTFT}$ (aka Fourier transform) is given by the equation:

$$
\mathrm{DTFT}\left\{x[n]\right\}=X(e^{j\omega})=X(\omega), \quad 0\leq n\leq N - 1.\tag{1}
$$

And the $N$-point $\mathrm{DFT}$ is given by uniformly sampling equation $(1)$ to obtain $N$ equally spaced frequencies $\omega_k = 2\pi k/N$:

$$
\mathrm{DFT}\left\{x[n]\right\}=X(e^{j\omega_k})=X\left(\frac{2\pi k}{N}\right)=X[k], \quad 0\leq k\leq N - 1.$$

If you plot $X[k]$, you're plotting the $N$ equally spaced samples from the $\mathrm{DFT}$ equation not against the continuous frequency $\omega$ but against the index number $k$ which in frequency corresponds to $2\pi k/N$.
The $\mathrm{FFTs}$ are algorithms to compute the $\mathrm{DFT}$. I don't know what's the goal in plotting the $\mathrm{DFT}$ against $k/F_s$. If you check the units you see that $k$ is just a bin number, unit-less, but $F_s$ is in Samples per second. So, that division has units in Seconds per samples...
The spacing between the discrete frequencies is indeed $2\pi/N$ radians. By the continuous-time frequency spacing, I think he meant the corresponding frequency spacing in $\rm Hz$ and this is the frequency resolution:
$$f_k=k\cdot \frac{F_s}{N}$$
From here you can guess what's the frequency spacing since two consecutive bins from the $\mathrm{DFT}$ samples are at $\Delta k=1$.


How to differentiate an empty aerosol from a full one using only sound measurements?

What I would do is :

1) Acquire the audio signal (or in fact, several, like ten) corresponding to a full container. I'd then quantify the energy of the signal, by squaring it and dividing by the number of samples. The obtained quantity is my reference energy. If several signals were measured, I take the average of calculated energies. Let's call it $E_{full}$

2) Do the same with an empty container (which is basically measuring ambient noise if it does not emit any sound like you say). Let's call it $E_{empty}$

From this point, you can define any energy threshold, with which you will compare the energy of the sound emitted by your experimental aerosol.

You can also do it simpler, by detecting for example the maximum peak amplitude of the full and empty aerosol and defining peak-amplitude threshold accordingly.

The threshold definitions can also be fine-tuned by measuring other aerosols with varying fillings, like half-empty, 25 % etc...

And finally, all the measurements I suggest implies that the sound is measured in steady state.


Record Signal after Threshold Reach in GNU Radio

Sounds like you have a very workable approach:

Write a GNU Radio, Embedded Python, out-of-tree Python or C++, doesn't matter, which:

is a general block (not a sync_block)
has a member property triggered or similar, which is initialized to False in the constructor
has a member property threshold or similar, which is initialized to the value passed as constructor argument in the constructor
If triggered, then

copy input to output (specifically, min(len(input_items[0], noutput_items)), because in a general block, there's no guarantee that in- and output buffers have the same free space)
consume_each the number of items copied (tell GNU Radio you've processed this much input), and
return the number of items copied (you produced this much output).

else,

look for a super-threshold value in the input (numpy.any(numpy.abs(input_items[0]) > threshold)
if no values are found:

consume_each(len(input_items[0]))
return 0 (you didn't produce any output)

if a value is found:

triggered = True
find the position of the first value larger than the threshold
consume_each as much items as are before that position
return 0 (you didn't produce any output, but your work method will be called again, with triggered set)

... and connect its output to a "Head" block, set to pass 10000 items; connect that to your file sink.


You can find what you need here where a comparaison of lot of matting algorithms is done. You can see the rank of those algorithms' performances as well as online demonstration on a natural and complicated image.

From that study, you can guess that A Global Sampling Method for Alpha Matting is the best algorithm to resolve your problem. The algorithm is conceived by A Global Sampling Method for Alpha Matting
Kaiming He, Christoph Rhemann, Carsten Rother, Xiaoou Tang and Jian Su.

Hope this helps.

Begueradj.