FFT time domain average vs frequency bin average

Let me clarify.

Fourier transform does not represent the histogram of the signal.  Fourier transform is a linear transform that takes signal from time domain (complex function) to frequency domain (another complex function).  It takes a complex function to another complex function.
Fourier transform is linear as the poster above pointed out.
Phase in your samples matters as pointed out above.  If trial-by-trial data varies in phase, then you do not want to average before doing a Fourier transform, but you also do not want to average after Fourier transform.  You want to average after Fourier transform and norm.  I will elaborate below as far as exactly what needs to be done.

The main issue here is that the question is posed wrong.  It is not "should I take the Fourier transform before averaging or after averaging".  Because it makes no difference due to linearity of Fourier transform.

The correct question to ask is "should I take the amplitude of the Fourier transform before averaging or after averaging".  For this question, the answer is before.

Here are the details.

Suppose your sampled data is represented by the sequences:

$d_1=d_1[n_1],d_1[n_2],...d_1[n_N]$

$d_2=d_2[n_1],d_2[n_2],...d_2[n_N]$

$d_3=d_3[n_1],d_3[n_2],...d_3[n_N]$

...

$d_M=d_M[n_1],d_M[n_2],...d_M[n_N]$

where $d_1,...d_M$ are data from M trials and $n_1,...n_N$ are sampled time points, then:

$F_1 = \sum_{j=1}^M{|\mathcal{F}\{d_j\}|} \neq  |\mathcal{F}\{\sum_{j=1}^M{d_j}\}| = F_2$

So while the transform $\mathcal{F}$ is linear, $|\mathcal{F}|$ is not.

Furthermore, while $d_j[n_i]$ is real for all $i,j$, $\mathcal{F}\{d_j\}$ is not, but $|\mathcal{F}\{d_j\}|$ is.

As for what you should do, you should take the Fourier transform of individual trials (via FFT), get the amplitude of individual trials, and a average them together.

Finally, what is $1/f$.  $1/f$ is a short term for the frequency spectrum of "natural" signals (usually people think of images).

When people say there is a large $1/f$ component, it means that the amplitude as a function of frequency looks like $1/f$.  It's totally hand-wavy... probably coming from a biologist :p

The inverse Fourier transform of $1/f$ is some sign function, but that is useless.  It is an imaginary sign function!  Real functions generate symmetric Fourier transform.

In fact saying that the spectrum is $1/f$, tells you something about the signal, but it doesn't let you recover the signal.  All that you know is that $|\mathcal{F}\{x(t)\}| = |1/f|$.  This doesn't let you uniquely determine $x(t)$ because all the phase information is gone, and we know that the structure of a signal relies heavily on its phase.

What does $1/f$ tell you?  Simply that it contains a lot of low frequency and a little high frequency.

Just as important a question, what does averaging buy you?  and more important is how to interpret the result?  Tune in tomorrow for a more in depth discussion :p


How Are Images Considered Non Stationary Signal When They Are Invariant to Time?

Stationarity is a statitical concept defined through stochastic processes. However the following is a practical account of why most typical images are nonstationary :

Throughout the image window, consider a small rectangular test block, say 16x16, and observe the frequency spectrum estimation of that block. Now move this test block accross the image and re-compute the resulting spectrum estimation each time. Now, if the spectrum does not change while the test block is moving accross the image, then you can loosely say that the image under concern is stationary.

Since for most natural images, such a moving block will be crossing unrelated distinct objects during its travesral of the image, the frequency spectrum will also be changing, i.e., image being nonstatinary

Basic Fourier type transforms are not very efficient and proper for this kind of local variations in spectra whereas wavelet or packets or short-space Fourier methods will be more efficiently indicating those local features.


Attempting to use the trapezoidal rule to form a difference equation representing a circuit

The Laplace transform of the original differential equation is
$$Cs V_o(s) = V_i(s) \left( \frac{1}{R_1} \right) - V_o(s) \left( \frac{1}{R_1} + \frac{1}{R_2} \right)$$
$$R_1 R_2 C s V_o(s) = R_2 V_i(s) - (R_1 + R_2) V_o(s).$$
In the Laplace domain, the ratio $V_o(s)/V_i(s)$ is the transfer function for the circuit,
$$H(s) = \frac{ V_o(s) }{ V_i(s) } = \frac{ R_2 }{ R_1 + R_2 + R_1 R_1 C s }
= \frac{K}{1 + R_p C s},$$
where $K = R_2 / (R_1 + R_2)$ is the voltage divider ratio for $R_1$ and $R_2$
in a series connection, and $R_p = R_1 R_2 / (R_1 + R_2)$ is the equivalent
resistant for $R_1$ and $R_2$ in a parallel connection.
We can obtain the $z$-domain transfer function for an equivalent digital filter
using the bilinear transformation,
$$s = \frac{2}{T} \left( \frac{z-1}{z+1} \right)$$.
In this transformation, the constant is set to $2/T$, which will produce a digital filter that is an approximation to the analog filter that is equivalent to approximating the convolution integral by the trapezoidal rule.  I assume this is what you are looking for.
By making this substitution, the transfer function for the digital filter is
$$H(z) = \frac{K}{1 + R_p C \left( \dfrac{2}{T} \right)
\left( \dfrac{z-1}{z+1} \right)}
= \dfrac{ KT (z+1) }{ T(z+1) + 2 R_p C (z-1) }$$.
This can be further simplified to
$$H(z) = \frac{V_o(z)}{V_i(z)} =
\frac{ \left( \dfrac{KT}{T + 2 R_p C} \right) \left( 1 - z^{-1} \right) }
{1 + \left( \dfrac{T - 2 R_p C}{T + 2 R_p C} \right) z^{-1}}.$$
Taking the inverse $z$-transform of this equation, we obtain a
difference equation for the digital filter,
$$Y(z) \left[ 1 + \left( \dfrac{T - 2 R_p C}{T + 2 R_p C} \right) z^{-1} \right]
= X(z) \left[ \left( \dfrac{KT}{T + 2 R_p C} \right) \left( 1 - z^{-1} \right) \right].$$
$$y[n] = \left( \dfrac{KT}{T + 2 R_p C} \right) \left( x[n] - x[n-1] \right)
- \left( \dfrac{T - 2 R_p C}{T + 2 R_p C} \right) y[n-1].$$
You should be able to use this equation to recursively compute $y[n]$.


Directly compare subpixel shifts between two spectra — and get believable errors

I think using cross-correlation and interpolating the peak would work fine.  As described in Is up-sampling prior to cross-correlation useless?, interpolating or upsampling before the cross-correlation doesn't actually get you any more information.  The information about the sub-sample peak is contained in the samples around it. You just need to extract it with minimal error.  I gathered some notes here.

The simplest method is quadratic/parabolic interpolation, which I have a Python example of here.  It's supposedly exact if your spectrum is based on a Gaussian window, or if the peak happens to fall exactly on the midpoint between samples, but otherwise has some error.  So in your case you probably want to use something better.

Here is a list of more complicated, but more accurate estimators.  "Of the above methods, Quinn's second estimator has the least RMS error."

I don't know the math, but this paper says that their parabolic interpolation has theoretical accuracy of 5% of the width of an FFT bin.

Using FFT interpolation on the cross-correlation output does not have any bias error, so that's the best if you want really good accuracy.  If you need to balance accuracy and computation speed, it's recommended to do some FFT interpolation and then follow it with one of the other estimators to get a "good enough" result.

This just uses the parabolic fit, but it outputs the right value for the offset if the noise is low:

def parabolic_polyfit(f, x, n):
    a, b, c = polyfit(arange(x-n//2, x+n//2+1), f[x-n//2:x+n//2+1], 2)
    xv = -0.5 * b/a
    yv = a * xv**2 + b * xv + c

    return (xv, yv)

arraysize = 1001
fluxlevel = 100.0
noise = 0.3 # 2.0 is too noisy for sub-sample accuracy
signal_std = 15.0
signal_depth = 40.0
gaussian = lambda x: np.exp(-(mu-x)**2/ (2 * signal_std))
mu = 500.1
a_flux = np.array([ra.normalvariate(fluxlevel, noise) for x in range(arraysize)] - gaussian(np.arange(arraysize)) * signal_depth)
mu = 500.8
b_flux = np.array([ra.normalvariate(fluxlevel, noise) for x in range(arraysize)] - gaussian(np.arange(arraysize)) * signal_depth)

a_flux -= np.mean(a_flux)
b_flux -= np.mean(b_flux)

corr = ss.fftconvolve(b_flux, a_flux[::-1])

peak = np.argmax(corr)
px, py = parabolic_polyfit(corr, peak, 13)

px = px - (len(a_flux) - 1)
print px

The noise in your sample produces results that vary by more than a whole sample, so I reduced it.  Fitting the curve using more of the points of the peak helps tighten the estimate somewhat, but I'm not sure if that's statistically valid, and it actually makes the estimate worse for the lower-noise situation.

With noise = 0.2 and 3-point fit, it gives values like 0.398 and 0.402 for offset = 0.4.

With noise = 2.0 and 13-point fit, it gives values like 0.156 and 0.595 for offset = 0.4.


Frequency bin of the positive and negative frequency

In the tutorial it is stated that

  Because the frequencies are spaced at 1/72 and the first bin corresponds to 0 frequency, the correct bin is 72/12+1. This is the frequency bin of the positive frequency. You must also include the frequency bin corresponding to the negative frequency: -1 cycle/12 months. With MATLAB indexing, the frequency bin of the negative frequency is 72-72/12+1.

So, regarding your questions:

It's not that there is only one bin for positive frequencies. The thing is that you are looking for a linear regression model, so you check what frequency is the most significant and then you create what would be the resulting vector if a DFT was made to a discrete sinusoid of that frequency only. That is your model.
The bin number is not any weird formula or anything. You have data for every month, but in the tutorial they are looking for cycles per year. You have 6 years. So $72/12 =6$ and that's where it comes from. (As you said, the $+1$ comes from the fact that there is a bin corresponding to the mean value.)


How to prove energy is preserved in sampling using Parseval's Relation?

Since your signal is band-limited, it has a representation in the orthogonal system of shifted sinc functions,

$$x_c(t)=\sum_{n\in\Bbb Z} x_c(nT)\operatorname{sinc}\left(\tfrac tT-n\right)$$

where $\operatorname{sinc}(x)=\frac{\sin(\pi x)}{\pi x}=\int_{-1/2}^{1/2}e^{i2\pi fx}\,df$.

By the Plancherel identity of the continuous Fourier transform this leads to
$$
\langle \operatorname{sinc}(x-m),\,\operatorname{sinc}(x-n)\rangle
=
\int_{-1/2}^{1/2}e^{i2\pi f(n-m)}\,df
=
\delta_{m,n},
$$
 i.e., orthogonality of the unscaled functions. In the scaled variant above one gets, first by orthogonality and then by rescaling of the time parameter, the Parseval identity for the orthonormal basis $\tfrac1{\sqrt{T}}\operatorname{sinc}(\tfrac tT-n)$, $n\in\Bbb Z$, of the space of band-limited functions,

$$
\|x_c\|_2^2=\sum_{n\in\Bbb Z} |x_c(nT)|^2 \langle \operatorname{sinc}(\tfrac tT-n),\,\operatorname{sinc}(\tfrac tT-n)\rangle=T\,\sum_{n\in\Bbb Z} |x_c(nT)|^2
$$


what is the maximum instantaneous frequency of Dirac delta?

It's not very sensible to even try to assign an instantaneous frequency to signals like yours.

Your specific example fails for more fundamental reasons than usual however, namely for the fact that the Dirac delta is not a function but something mathematicians call a "distribution".

This for example implies that we cannot assign a value to $\delta(0)$ even though we know $\delta(x)=0$ for all $x\neq0$. The definition of the delta distribution is, that $\int_\mathbb{R} \delta(x) f(x) dx=f(0)$. From here it follows that $\delta(0)$ cannot have a finite value. However, stating that it is "infinite" there is also not enough. There's simply no way to think of any value at $0$ that makes sense.

So the answer to your question is that the instantaneous frequency is $0$ almost everywhere, and undefined at $t=0$. So the maximum (and minimum!) value of $\omega$ is $0$!

You can avoid this situation by replacing the unit step with a smooth approximation of a unit step. With this you can give your maximum instantaneous frequency any value you desire, even if you take a limit of a function sequence that results in the unit step itself. So the answer is absolutely arbitrary.


Root Locus, Transfer Functions and Unit Step Response?

The closed loop poles are the roots of the polynomial
$$D(s)=s^2+2s+2+K\tag{1}$$
and, according to the root locus plot, they are $s_{1,2}=-1\pm 2j$. Consequently, we get
$$2+K=|1+2j|^2=5\quad\Longrightarrow\quad K=3\tag{2}$$
With $K=3$ we obtain $$H(0)=\frac{K}{2+K}=\frac35\tag{3}$$
which leaves step response $C$ as the only option.


Question about Difference of Gaussian (DoG) algorithm

Difference of gaussian is the difference in the output of two Gaussian filters with different blur amounts (sigma).

Sigma is the size of the Gaussian filter. A bigger sigma gives you a bigger amount of blurring. A good way to think about it is a Gaussian filter with variance sigma is very roughly like averaging 3 x sigma samples wide (or 3 x 3 in an image)

e.g. from wikipedia:

Very Important when making a Gaussian filter in MATLAB make sure the size of the filter is at least 6 x sigma. In your above code you have 5 x 5 which is fine for sigma = 0.5, but for sigma = 1 you would want 6 x 6 or bigger.

The k is simply a multiplier for sigma.

e.g.

sigma = 0.5;
gauss1 = fspecial('gaussian', round([10*sigma 10*sigma]), sigma);
sigma = 1;
gauss2 = fspecial('gaussian', round([10*sigma 10*sigma]), sigma);
blur1 = imfilter(img, gauss1, 'replicate', 'same');
blur2 = imfilter(img, gauss2, 'replicate', 'same');
dog2 = blur1 - blur2;

A more complete code example that allows you to set the number of octaves and the steps per octave:

%% Filter using DoG
stepsPerOctave = 5;
octaves = 4;
mult = nthroot(2,stepsPerOctave);

% Create blurry images
sigma = 0.5;
kernelSize = [10*sigma*2^(octaves),10*sigma*2^(octaves)]
for k = 1:octaves*stepsPerOctave+1
    disp(['Sigma is ' num2str(sigma)]);
    gauss = fspecial('gaussian', kernelSize, sigma);
    blur(:,:,k) = imfilter(I, gauss, 'replicate', 'same');
    imagesc(blur(:,:,k)); colorbar; title(['Gaussian ' num2str(k)]); pause;
    sigma = sigma * mult;
end

% Create DoG
for k = 1:octaves*stepsPerOctave
    dog(:,:,k) = blur(:,:,k+1) - blur(:,:,k);
    imagesc(dog(:,:,k)); colorbar; title(['DoG ' num2str(k)]); pause;
end


Assuming a channel whose input at each time is a continuous random variable $X$ and its output is $Y=X+Z$, where $Z\sim\mathcal{N}(0,N)$ and $Z$ is independent of $X$, then $$C_{\text{CI-AWGN}}=\frac{1}{2}\log_2\left(1+\frac{P}{N}\right)$$ is the capacity of the continuous-input channel under the power constraint
$$\mathsf{E}X^2\le P$$
The mutual information $I(X;Y)$ is maximized (and is equal to $C_{\text{CI-AWGN}}$) when $X\sim\mathcal{N}(0,P)$.

This means that if $X$ is a continuous Gaussian random variable with the given variance, then the output has the highest possible mutual information with the input. That's it!

When the input variable $X$ is discretized (quantized), a new formulation is required. Indeed, things can easily become difficult. To see it a little bit, one can consider the simple case of a very coarse discritization of $X$ where it can only have two values. So assume that $X$ is selected from a binary alphabet, for instance let $X\in\{\pm1\}$ (or a scaled version to satisfy a power constraint). In terms of modulation, it is identical to BPSK.

It turns out that the capacity (even in this simple case) has no closed form. I report from "Modern Coding Theory" by Richardson and Urbanke:

$$\begin{align}C&_{\text{BI-AWGN}}=1+\\&\frac{1}{\ln(2)}
\left(\left(\frac{2}{N}-1\right)\mathcal{Q}\left(\frac{1}{\sqrt{N}}\right)-\sqrt{\frac{2}{\pi N}}e^{-\frac{1}{2N}}+\sum_{i=1}^{\infty}\frac{(-1)^i}{i(i+1)}e^{\frac{2i(i+1)}{N}}\mathcal{Q}\left(\frac{1+2i}{\sqrt{N}}\right)\right)\end{align}$$
A comparison between the two cases can be seen in the figure below:


What does orthogonality mean in the context of adjacent channels in a filter bank?

The mathematical definition of orthogonality between two vectors is that their dot product is zero.  It just means that there is no correlation between the two- at least at that "phase".  It is often the case that if you shifted one of the vectors you would get strong correlation.

Infinite vectors of different frequencies are always orthogonal, so in an ideal world the output of adjacent channels in a filter band would always be orthogonal.  There are two ways, though, that the real world is not ideal.

First, time limitations can introduce non-orthogonality.  The non-orthogonality is usually trivial if both of the vectors have "many" sinusoidal cycles, but can be substantial if the length is less than a cycle.

Second, non-ideal filters means that attenuated stop-band frequencies get into the filter output, which means that the adjacent channels do have frequencies in common, just at different amplitudes.


Blind Estimation of Signal Parameter and Noise Variance

Squaring Method

The input signal $x[n]$ can take two values: $+1$ or $-1$. After being multiplied by $h$, the signal becomes $|h|e^{j\angle h}$ (if $x[n]=+1$) or $|h|e^{j(\pi + \angle h)}$ (if $x[n]=-1$). By squaring this we get, $|h|^2e^{j2\angle h}$. Both points get mapped to this. Now we can pick out the magnitude and phase of $h$. In summary, the steps are:

Square $y[n]$ to get $z[n]=y^2[n]$.
Get the magnitude, $\hat{|h|}=\sqrt{\bigg|\frac{1}{N}\sum_n z[n]\bigg|}$. It was pointed out in the comments that this is a biased estimate with bias $\sigma^2$, so it is not favorable for a large $\sigma$. I describe another way to estimate $\sigma^2$ (without using $\hat{h}$) so that we can use this $\hat{\sigma}^2$ to remove the bias from $\hat{|h|}$.
Get the phase, $\angle h=\frac{1}{2N} \sum_n \angle z[n]$
Use $\hat{h}$ to get the noise variance, $\hat{\sigma}^2=\bigg(\frac{1}{N}\sum_n z[n]\bigg)-\hat{h}^2$. This comes from the fact that $E\big[z[n]\big]=h^2+\sigma^2$ (also see the clustering method is a better way to go about this as it does not depend on any other estimates).

Caveat: This method is valid only for phase offsets within $\pm \frac{\pi}{2}$. Anything beyond that, the received signal becomes indistinguishable and there are multiple phase shifts which could give the same signal.

Clustering Method

If $h$, can be correctly determined, then you can take hard decisions and calculate the noise variance using the distance of the received symbols from the decision you made.

If $h$ can't be determined right away, then this clustering based estimate of $|h|$ and $\sigma^2$ can be used: We know there should be two clusters ($\pm 1$) so we can use a simple method like k-means (https://en.wikipedia.org/wiki/K-means_clustering) with $k=2$. The algorithm is:

Run k-means on the received samples and get a list of labels back labeling each sample to either cluster $1$ or $2$. Let $S_1$ be the set of sample indices in cluster $1$ and $S_2$ be the set of sample indices in cluster $2$.
Now we are going to take all the points in each cluster, shift them be zero mean and take the variance (this will give us a noise variance estimate for each cluster). Then we will average over all of the cluster noise variance estimates to form the final estimate. $\hat{\sigma}^2=\frac{\bigg( \text{var}\big(y[n]-E[y[n]] \big)\big|_{n \in S_1} + \text{var}\big(y[n]-E[y[n]] \big)\big|_{n \in S_2} \bigg)}{2}$
Now we can take the estimate from step $2$ and subtract off the bias term which we estimated here as $\hat{\sigma}^2$ to get the refined estimate $\hat{|h|}_0=\hat{|h|}-\hat{\sigma}^2$. You could also get $|h|$ using the magnitude of the cluster means, $|\hat{h}|=\frac{\text{E}\big(y[n]\big|_{n \in S_1}\big)+\text{E}\big(y[n]\big|_{n \in S_2}\big)}{2}$.

There is also a way get $\angle h$ from the clusters but it too has the criteria that $\angle h$ must be less than $\pm \frac{\pi}{2}$ so I will not detail it here.

Long story short: $|h|$ and $\sigma^2$ are doable, but $\angle h$ is only doable for some cases.

Code: https://github.com/B-William/DSPSE/blob/master/blindEstimationOfSignalParameterAndNoiseVariance.m


How can faux correlated data be generated for testing or training?

Remember that the power spectral density (PSD) of anything is the Fourier transform of the autocorrelation function of that process.

Hence, if you want a defined autocorrelation function (ACF), do its inverse Fourier transform, and you'll get a PSD that would represent a signal of given ACF.

Now, you can use all the beautiful methods of filter design to come up with a filter that, when applied to a white process (e.g. random numbers from a random number generator) will shape that process's PSD as desired!

So, that's kind of the general look at things: If you just want something like

  "Each sample should be correlated much with the previous one, but only slightly with the one before"

then maybe understanding it as roughly

$$y[n] = 0.9x[n] - 0.1x[n-1]$$

would get you there quicker: $(0.9, 0.1)$ is actually the coefficient vector of a FIR filter that you can apply to your uncorrelated input vector $x$ and be done with it.


Why Wavelet based Transform Is More Suitable for Image Compression Compared to DCT?

Both JPEG and JPEG 2000 use the change of basis compression type.

Namely, we transform the data into a different representation assuming in this representation the number of parameters needed to describe to data is lower.
Or to the least, most of the information is gathered within few parameters.

Now, if you look at the energy level of the DCT coefficients of real world images you'd see most of the energy is limited to the very few coefficients.
This is what JPG does, keeps only the few dominant coefficients and throws the rest.

How efficient is DCT?
It depends on how less coefficients are needed to describe the image in an acceptable quality.

It turns out that for Wavelets the situation is better.
Namely on real world images less coefficients are needed to describe the image with the same perceived quality.

This property is called the ability to decorrelate the data / being energy densed, etc...

At the end, it is all up to how much less parameters are needed to describe the same data.

P.S.
Using SVD / KLT (They are equivalent in this sense) would be even more efficient in the sense of how much less coefficients are needed for the same quality (In the $ {\ell}_{2} $ sense).
The problem is those bases are adaptive to the data and hence if you use them you have to include them to be able to decompress the data.
What's so nice, you can ask which, in general, approximate the SVD better?
It turns out that Wavelets are better doing so.


A signal either has finite energy, finite power or even infinite power. If it has finite energy, it will have zero average power, according to your definition
$$P_x=\lim_{T_0\rightarrow\infty}\frac{1}{T_0}\int_{-\frac{T_0}{2}}^{\frac{T_0}{2}}|x(t)|^2dt.$$

Knowing that the sincs are orthogonal to each other, as well as the cos function, the average power is given by

$$P_x=P(2)+P(2\cos(2\pi f_0t))+P(\operatorname{sinc})$$

where $P(...)$ is the power of the components in brackets. The sinc function has finite energy, hence it has zero power on average. Hence, we have

$$P_x=P(2)+P(2\cos(2\pi f_0 t)$$
with
$$P(2)=4$$
and $$P(2\cos(2\pi f_0 t)=2^2\cdot\frac{1}{2}=2.$$

Therefore, the overall power is given by

$$P_x=6.$$

For the calculation of the $\cos$, consider
$$P_{\cos}=\frac{1}{T}\int_{0}^{T}\cos^2(2\pi t/T)dt=\frac{1}{2}$$


Is a Fourier transform a sound way to analyse a transient signal?

If the signal is zero-valued outside the time window being analyzed ("analysis window"), then the discrete Fourier transform (DFT) exactly samples the discrete-time Fourier transform (DTFT) of the signal at harmonic frequencies. DTFT implies no time-domain periodicity. The sampling of DTFT by DFT is evident from definitions of the two transforms, if we apply a finite support signal constraint:
$$\text{DTFT:}\quad X(\omega) = \sum_{n=-\infty}^{\infty} x[n] \,e^{-i \omega n}$$
If $x[n] = 0$ for all $n$ such that $n < 0$ or $n \ge N$ (finite support signal constraint), then:
$$= \sum_{n=0}^{N-1} x[n] \,e^{-i \omega n}$$
The frequency variables of the two transforms are related by $\omega = \frac{2\pi}{N}k$, so:
$$= \sum_{n=0}^{N-1} x[n] \,e^{-\frac {i 2\pi}{N}kn}$$
This is identical to the definition of DFT:
$$\text{DFT:}\quad X[k] = \sum_{n=0}^{N-1} x[n] e^{-\frac {i 2\pi}{N}kn}$$
Choosing a longer window (larger $N$) enables DFT to sample DTFT more densely.
For a signal that has non-zero values outside the analysis window, doing the analysis on just the window is equivalent to first multiplying the signal by a rectangular window, and then sampling the DTFT of this windowed signal using DFT. There are other window functions that may have more favorable characteristics for the type of signal and for the particular analysis task, allowing to isolate a transient of interest without too much distortion in time domain or smearing in frequency domain.
Fast Fourier Transform (FFT) is often used because it is very fast to compute, even though a decomposition into harmonic sinusoids might not always be the best way to describe a signal.


How to enhance Pedestrian Dead Reckoning using Particle filter?

How to determine the initial position pdf distribution? I have a Bluetooth beacon to detect the starting position. Particle filter requires a cluster of point at the start, right?

No, the distribution of a particle filter does not need to be a cluster when it starts.  A classic case is for a robot to figure out where it is.  At first it does not know where it is.  It takes sensor readings wherever it happens to be at, but that is usually not enough to localize it (if all you needed was a single set of readings, you wouldn't need a particle filter).  Instead, it moves around and takes more readings.  All of the many spread out particles represent possible positions.  As it moves and takes more readings more and more of the possible positions become highly unlikely, until the distribution collapses to a single cluster.  At that point the robot "knows" (or at least thinks it knows) where it's at.

  I don't think in my case the distribution is Gaussian. Then what should it be?

Without giving us more information there is no way to answer that question.  Since it sounds like your initial estimate comes solely from your bluetooth beacon, you will probably need to generate an empirical data-based model of its accuracy.

  Should I do the particle filtering on the smartphone or pass the data back and do it on the computer using Matlab? Will the computation very intense.

Either way.  The computational load will largely depend on how many particles you use.  There is no right answer as to how many particles you need.  More particles means more computation but also lower chance of getting a wrong answer.


Chroma-based approaches are used when the tracks to align are two performances of the same material with widely different instrumentations - in which case you need to find the "lowest common denominator" which is harmony, melody, chord progression etc. Chroma are good for that - they totally abstract timbre, recording techniques. But you don't need those in your case.

As you have guessed, cross-correlation of the raw audio would be 1/ not very robust to noise, 2/ not very robust to differences in transducers (the two microphones might have very different responses) 3/ very costly if you have minutes of audio.

I suggest you to:

Extract a sequence of MFCC vectors from both signals. This way you'll have something of lower dimensionality and a bit more robust to noise and differences in transducers.
Optionally normalize the MFCC ("by design" coefficient 0 has more variance than coefficient 1 and so on...)... Ideally you'd want each row of the MFCC matrix to have variance 1.
Compute the cross-correlation of two matrices you get, along the time axis (that is to say, compute the 2D cross-correlation and just keep the data for the time axis).


Channel Vocoder producing output with "click" sounds

I think it could be one of two issues.

What is your frame overlap? if the frame overlap is something less than 40% or so you will begin to get undesired artefacts when you perform overlap add synthesis at the end. Try an overlap of 50-80% (although more computationally expensive, it provides greater frequency resolution).

The second consideration which goes hand in hand with the above suggestion. Are you using any kind of analysis window on your framed signal? Perhaps try using a tapered window such has hamming, Blackmann or Hanning. This will help minimise spectral leakage between frames. Hope this helps.


When is a network called end-to-end training?

From feature extraction to learning the desired result, deep learning algorithms can act as full pipelines for solving tasks at hand. End-to-end learning usually refers to omitting any hand-crafted intermediary algorithms and directly learning the solution of a given problem from the sampled dataset. This could involve concatenation of different networks such as multiple CNNs and LSTMs, which are trained simultaneously.

For the OCR example, instead of trying to classify characters and clustering them into words, it is shown to be a better approach to directly use CNNs to regress the words themselves.

In the optical flow, one could go directly from images to the final flow field, omitting the image difference computations.

In self driving cars, the network can be trained to directly learn how to drive.

In all such examples, the idea is to let the network go from the "raw-est" possible data to the final-most output. This is found to perform better. End-to-end learning reduces the effort of human design and performs better in most applications.


How do I implement a footstep recognition algorithm? One that recognizes the start and end points of each footstep?

This signal looks pretty clean.  I would stick to time domain feature extraction.  Maybe 1st and 2nd derivative of time domain data.  If you're just counting, you can almost just use the magnitude and see that these steps have the most negative values or acceleration.  You could use the negative peaks to center a search for the edges of the "step", and to discriminate against the other negative peaks that are in between the steps.


Sample Frequency - Recommended Practical Nyquist, Lower Used by majorities appications

44.1 kHz and 48 kHz are sampling frequencies that are almost always supported by software and drivers but if you need to go lower than that then stop at 8 kHz, as it is commonly supported.


Estimating the Signal by Deconvolution with a Prior on the Filter Coefficients and the Signal Samples

I would take approach based on Blind Deconvolution.

Since we're dealing with ill posed problem some assumptions should be made.
The intuitive approach would be using the information as a prior for the signal. Another idea is to add LPF assumption of the Filter by setting the sum of its coefficients to be 1 and non negative. Yet since we have Discrete Prior on the signal we're getting into combinatoric problem.
Which means brute force solution where the number of combinations is $ {2}^{n} $ where $ d $ is the number of the signal samples.

For $ n \leq 16 $ I'd say it would work for the given input size.
Yet for solution with higher number of samples this method isn't feasible.

In order to deal with higher dimensions (More samples) I'd use a GMM:

Namely the Prior model is a Gaussian Mixture Model (GMM) with 2 Gaussian's centered at $ \left\{ -1, 1 \right\} $ with very small variance in order to approximate discrete probability function.

So the problem I'm looking to solve is given by:

$$\begin{aligned}
\arg \min_{h, x} \quad & \frac{1}{2} {\left\| h \ast x - y \right\|}_{2}^{2} \\
\text{subject to} \quad & \sum h = 1, \, h \succeq 0
\end{aligned}$$

Step 1 - Solving for the Filter $ h $

Given the signal $ x $ is known, solving for the filter is pretty easy using matrix form of the problem (Which is convex):

$$\begin{aligned}
\arg \min_{h} \quad & \frac{1}{2} {\left\| X h - y \right\|}_{2}^{2} \\
\text{subject to} \quad & \sum h = 1, \, h \succeq 0
\end{aligned}$$

I showed, code included, how to solve such problem in my answer to How to Project onto the Unit Simplex as Intersection of Two Sets (Optimizing a Convex Function)?

Step 2 - Solving for the Signal $ x $

The model is $ y \mid h \sim \mathcal{N} \left( h \ast x, {\sigma}_{n} I \right) $ and the prior $ {x}_{i} \sim 0.5 \mathcal{N} \left( {\mu}_{1} = -1, {\sigma}_{1}^{2} = {0.1}^{2} \right) + 0.5 \mathcal{N} \left( {\mu}_{2} = 1, {\sigma}_{2}^{2} = {0.1}^{2} \right) $.

I'd use the MAP so we have:

$$\begin{aligned}
\arg \max_{x} p \left( x \mid y \right) & = \arg \max_{x} p \left( y \mid x \right) p \left( x \right) \\
& = \arg \max_{x} \log p \left( y \mid x \right) + \log p \left( x \right) \\
& = \arg \min_{x} -\log p \left( y \mid x \right) - \log p \left( x \right) \\
& = \arg \min_{x} \frac{1}{2} {\left\| h \ast x - y \right\|}_{2}^{2} - \lambda \log p \left( x \right) \\
& = \arg \min_{x} \frac{1}{2} {\left\| h \ast x - y \right\|}_{2}^{2} - \lambda \sum \log p \left( {x}_{i} \right)
\end{aligned}$$

Where $ \lambda \propto N {\sigma}_{n}^{2} $ where $ N $ is the number of samples (Dimension of $ y $).

This can be solved by any solver. I'd use MATLAB's fminunc().
Though one could alter (For MAP Estimation) the Expectation Maximization (EM) process for faster and better converging algorithm. Another option would be using Probabilistic Programming with one of the options available today.

Remark: The above is Bayesian Modeling of the problem. One could build optimization problem with some intuition in the form of:

$$ \arg \min_{h, x} \frac{1}{2} {\left\| h \ast x - y \right\|}_{2}^{2} + \lambda \sum_{i = 1}^{m} {\left( {x}_{i}^{2} - 1 \right)}^{2} $$

Which isn't motivated by Bayesian model but still drives the solution to where we want it.

Initialization

One approach to initialization of the estimated $ x $ signal is to use hard threshold. So we set $ {x}_{i} = 1 $ if $ {y}_{i} \geq 0 $ and $ {x}_{i} = -1 $ otherwise.

This approach could also be a greedy method to solve the step for $ x $. Yet it doesn't take under account the delay of the filter.

After iterating enough for a stable solution one could round the result of $ x $ such that $ {x}_{i} \in \left\{ -1, 1 \right\} $.

I haven't tested this approach myself, but I really like its model.
I will publish the results of MATLAB simulation soon.


Fast approximate optical flow / image shift

This might be a slow terrible solution, but you could do a FFT-based cross-correlation of subsequent frames and then find the peak to identify the offset between frames.  Maybe only do it on a small subset of the image to save processor cycles.

It would not work with rotation or drastic scene changes from one frame to the next, and there are probably better methods.  This is kind of a "I have a hammer so everything looks like a nail" solution.  I guess this is just like your naive solution, except there's no need for edge detection and the FFT makes it a lot faster than explicitly shifting one pixel at a time.

This question is similar, and no one's suggesting anything other than cross-correlation, so maybe it's not so bad: Using MATLAB to calculate offset between successive images


If you are optimizing for engineering time and are on a platform that supports large FFTs well (i..e not fixed point), then take hotpaw2's advice and use fast convolution.  It will perform much better than a naive FIR implementation and should be relatively easy to implement.

On the other hand, if you have some time to spend on this to get the best implementation or are on a fixed-point platform, you should use a multirate down-filter-up-subtract structure.  But it's a bit trickier to get everything right.

I've got access to trusted and highly optimized implementations of both fast convolution and multirate filtering tools. The fast convolution takes about 3x longer to get equivalent signal performance compared to the multirate structure.  Furthermore, that is even on a floating point platform. The gap would widen considerably on a fixed point dsp.

In general terms:

Down-conversion:

Use 8 stages of halfband,decimate-by-2 filters to transform your 48kHz signal into a 187.5 Hz signal.  The first stage of this downsampling can have a very wide transition band, allowing energy to alias as long as it doesn't alias back into the sub 60 Hz range. As the stages progress, the number of taps needs to increase, but they will be applied at a progressively lower sampling rates, so the overall cost of each stage remains small.

Filtering:

Perform your tight filtering around the 60 Hz bw to keep the energy you will eventually want to subtract.  There is a double advantage to doing your tight filtering at the low rate:

1Hz of transition bandwidth is 256 times larger in terms of digital frequency at the low rate vs. the original rate. So every tap of your filter is 256 times as powerful.
The signal itself is at a lower rate, so the filter only needs to process 1/256 the data.

Up-conversion:

Essentially, this is the reverse of the decimation stages.  Each of the 8 interpolator stages doubles the rate by estimating the sample that goes between consecutive input samples.  The transition band gets wider as the sample rate gets higher.

Subtract:

Subtract your full-rate low-pass filtered signal from the original signal. If you've adjusted for all the group delays properly, the overall structure will be a highpass filter with a narrow transition bandwidth.


I hope I understand your question correctly, in that you're trying to produce a mirror image of the convolution kernel (filter) and then convolve. Flip you convolution kernel. In MATLAB, you can use the flip command. If you flip it left-to-right, this should do it.

However, if you're simply saying that you want to use the same kernel but start the convolution process from a different point, this would not give you a different result back at all. This is one of the reasons you can do FFT based convolution.


Why use $\chi^2$ test to determine the presence of white noise?

If I and Q are normal, then their magnitude square is $\chi^2$ distributed.

Hence, if you check the distribution of the instantaneous power against that distribution, you get information about how much the complex baseband signal is circularly complex normal.

The autocorrelation function has power of the signal as its value for zero shift. For whiteness, the sample autocorrelation approaches an accordingly weighted Dirac impulse with growing observation length.


How to downsample a matrix in columns ? MATLAB

Maybe Off-Topic here, but anyway !

Try it

y = downsample(x',2)'


Is initializing a digital filter's output with no "momentum" a non-trivial task?

If you want a stable linear time-invariant system to output constant $y$, it must have received input $x$ that is the ratio of the constant output and the zero frequency response of the system $H(1):$

$$y = H(1)x\quad \Leftrightarrow \quad x = \frac{y}{H(1)}$$

Or, you'd want to change the state of the system so that it reflects that situation. The system must satisfy the condition $H(1) \ne 0.$

For a composite filter that consists of a serially connected cascade of second-order sections, you can start with the first section with input $x$ and calculate sequentially using $y = H(1)x$ (recycling the variable names to mean the inputs and outputs of a single section) the constant intermediate output of each section. Or you could start from the other end and use $x = y/H(1).$ Knowing both the input $x$ and output $y$, use the signal flow diagram of each section to calculate its state variables, starting with the summation point with dependency only to known constants. Because all oscillation has settled due to constant input and stability of the filter, all signals are constants rather than functions of time, and the delays in the flow diagram have become identities.

Figure 1. Transposed direct form II biquad signal flow diagram. (CC BY-SA 3.0 by Fcorthay)

Given constant input $x$, the constant output $y$ of a biquad section decipted in Fig. 1 is:

$$y = H(1)x = \frac{b_0+b_11^{-1}+b_21^{-2}}{1+a_11^{-1}+a_21^{-2}}x = \frac{b_0 + b_1 + b_2}{1 + a_1 + a_2}x$$

As for the state variables, the top sum $m_0,$ the middle sum $m_1,$ and the bottom sum $m_2$ are calculated as (start calculation from the bottom):

$$\begin{align}
m_0 &= b_0x + m_1\\
m_1 &= b_1x - a_1y + m_2\\
m_2 &= b_2x - a_2y\\
\end{align}$$

The top sum $m_0$ should equal $y$, which we can verify:

$$\begin{align}m_0 &= b_0x + b_1x - a_1y + m_2\\
&= b_0x + b_1x - a_1y + b_2x - a_2y\\
&= b_0x + b_1x - a_1\frac{b_0 + b_1 + b_2}{1 + a_1 + a_2}x + b_2x - a_2\frac{b_0 + b_1 + b_2}{1 + a_1 + a_2}x\\
&= \frac{b_0 + b_1 + b_2}{1 + a_1 + a_2}x\\
&= y
\end{align}$$


Why does the unit step function not have infinite poles?

No. Every point that's not in ROC doesn't have to be a pole. ROC is a region where the transform attains a finite value, but this ROC also has to "produce back" the time domain signal via the inverse transform, or else it does not converge to that sequence in time domain. So you have to take both the finiteness of the transform as well as the reproducibility of the time domain signal via the inverse transform into account to define a "region that converges and is consistence for both aspects" while defining an ROC for a signal. I am sure you know that the same transfer function with different ROCs produces different time domain signals. So specifying a transfer function does not uniquely define a time domain signal. Only when an ROC is defined along with the transfer function can we associate to a unique time domain sequence. So the left over region that's not part of the ROC doesn't have a pole, since clearly if we shift the ROC to be that side now, we can define a different sequence.


You need to build a time varying delay, where you can modulate the delay amount over time.

The peak delay modulation is a function of your maximum desired frequency shift and the modulation frequency.

This is not trivial since it will require fractional sample delays with some kind of interpolation algorithm. You can't round to the nearest integer delay without "zipper" noise. The best interpolation method depends on the specific requirements of your application (tolerance for noise, input spectrum, MIPS & memory available, etc.)


Cheaply generating coefficients for IIR/FIR given a cutoff frequency

A computationally cheap way to generate quick & dirty FIR filter coefficients is to evaluate a windowed Sinc function.  Any rectangle in the frequency domain (low/high/bandpass) has a Sinc impulse response in the time domain, which can be windowed to a length roughly proportionally inverse in size to the desired width of the transition bands, and then sampled for some FIR coefficients.

About a dozen lines of BASIC pseudo code on my web page to do this should be trivial to convert to plain C : http://www.nicholson.com/rhn/dsp.html#2


How to detect if the given signal is Data or Clock signal type?

Clocks are very regular, i.e. all the up and down transitions are equally spaced. Data tends  to more irregular, i.e. spacing between transitions varies.


The Y-axis units are the same as the Y-axis units you would have if you plotted the time-domain signal.  For instance, if the samples that you input into the FFT are the signal voltage, then the FFT units are also voltage.  To get energy units you have to square the absolute value of the FFT.


Conceptual question on FFT/IFFT (IFFT existence)

The conditions you mention only apply to the continuous-time Fourier transform, and they only have theoretical relevance. All signals occurring in engineering applications will satisfy these conditions. But if you want an example of a function with an infinite number of maxima and minima on a finite interval, here it is:

$$x(t)=\sin(1/t)$$

This function has an infinite number of minima and maxima on the interval $(0,a]$ for any real-valued $a>0$.

And as for your measured signals, they are all discrete-time signals of finite length, and the above conditions do not apply. So you don't need to worry about them when using the DFT/FFT, because the DFT is defined by a finite sum, unlike the continuous-time Fourier transform, which is defined by an improper integral.


Terminology: Camera Tracking, Camera Relocalization

In the context of SLAM papers, camera tracking means "local relative pose estimation of different views of the scene viewed by a given camera". It's also called visual odometry. What is important in this task is to get accurate estimates of from-to-frame (or in a small time window) camera displacements. It's to make a difference with mapping where you have global map optimization applied (think loop closures or global bundle adjustment). I refer you to the PTAM paper for more details as to why these 2 tasks are considered.
Relocalization means your robot has a map of a known environment, then you kidnap and blind it, and when you uncover its sensors you want it to recognize what is its new position with respect to the map it already knew (you did not bring it to an unknown place).


Understanding Soft Thresholding Operator

The soft-thresholding function finds the minimizer of an objective function  that involves data fitting in an $\ell_2$ sense as well as minimization of the $\ell_1$ norm (i.e. absolute value). The Lecture Notes - Penalty and Shrinkage Functions for Sparse Signal Processing gives a good discussion of how the soft threshold function is derived.


Does the audio signal in the time domain change much after a RMS from one person to another?

In short, yes.

One thing that you'll find is that the human ear is actually relatively insensitive to phase, so you can take an audio sample, apply an allpass filter that doesn't affect the signal's magnitude spectrum but distorts its phase, and it won't make it sound much different to you. However, the time-domain waveform could look completely different. I'm not sure what techniques the researchers used in the link you gave, but I would guess that it isn't likely to be based upon the time-domain structure of the audio waveforms.


Bag of Tricks for Denoising Signals While Maintaining Sharp Transitions

L1 norm minimization (compressed sensing) can do a relative better job than conventional Fourier denoising in terms of preserving edges.

The procedure is to minimize an objective function

$$
    |x-y|^2 + b|f(y)|
$$

where $x$ is the noisy signal, $y$ is the denoised signal, $b$ is the regularziation parameter, and $|f(y)|$ is some L1 norm penalty. Denoising is accomplished by finding the solution $y$ to this optimization problem, and $b$ depends on the noise level.

To preserve edges, depending on the signal $y$, you can choose different penalties such that $f(y)$ is sparse (the spirit of compressed sensing):

if $y$ is piece-wise, $f(y)$ can be total variation (TV) penalty;
if $y$ is curve-like (e.g. Sinogram), $f(y)$ can be the expansion coefficients of $y$ with respect to curvelets. (This is for 2D/3D signals, not 1D);
if $y$ has isotropic singularities (edges), $f(y)$ can be the expansion coefficients of $y$ with respect to wavelets.

When $f(y)$ is the the expansion coefficients with respect to some basis functions (like curvelet/wavelet above), solving the optimization problem is equivalent to thresholding the expansion coefficients.

Note this approach can also be applied to deconvolution in which the objective function becomes to $|x-Hy|+ b|f(y)|$, where $H$ is the convolution operator.


$$ y[n] = x[n-2] - 2 x[x-8] $$

that's actually a linear and time-invariant system.  you invert it by computing the Z transform

$$ Y(z) = z^{-2}X(z) - 2 z^{-8}X(z) $$

solve for $X(z)$, then inverse Z transform.  you might find out that, to invert this system, you will have to look at future samples, so the inverse would not be causal.


You are right. Even though the phase jumps at the zeros of the frequency response, such a phase response is usually still called "linear". For a frequency selective filter with frequency response zeros in the stopband, the phase always has discontinuities at the zeros. A purely linear phase response (without jumps) is only possible for filters with no zeros in their frequency response.

The frequency response of your filter can be written as

$$H(e^{j\omega})=|H(e^{j\omega})|e^{j\phi(\omega)}\tag{1}$$

where $\phi(\omega)$ is the phase response with jumps at the zeros, as shown in your figure. For such a linear phase system, the frequency response (1) can equivalently be written as

$$H(e^{j\omega})=A(e^{j\omega})e^{j\hat{\phi}(\omega)}\tag{2}$$

where $A(e^{j\omega})$ is a real-valued but bipolar (i.e. positive and possibly negative) function satisfying $|A(e^{j\omega})|=|H(e^{j\omega})|$. The phase $\hat{\phi}(\omega)$ in (2) is now a purely linear function without any jumps.


Butterworth low pass filter zeros location after bilinear transformation explanation

An $N^{th}$-order analog prototype system results after bilinear transform in an $N^{th}$-order discrete-time system. So the number of zeros and poles remains the same. All zeros at $|s|\rightarrow\infty$ map to $z=-1$. In your case there will be $6$ zeros at $z=-1$ because there are $6$ zeros at $|s|\rightarrow\infty$.

Of course, the magnitude (squared) response is symmetric, but the zeros at $|s|\rightarrow\infty$ cause the zeros in the magnitude at $\Omega\rightarrow\infty$ as well as at $\Omega\rightarrow -\infty$ (you can imagine those points meet at infinity). This becomes more obvious in the $z$-domain where moving towards positive frequencies means moving counterclockwise from $z=1$ (i.e., $\omega=0$), whereas moving towards negative frequencies means moving clockwise from $z=1$. In both cases you'll end up at $z=-1$.


On the Measurement Matrix Used for Compressing Sensing

It is indeed possible to formulate this setting in terms of matrix-vector products. First, let us re-formulate your $x$ (notice throughout that I use bold letters for vectors and matrices):
$$x = \begin{bmatrix}\mathbf{x}_1 & \mathbf{x}_2 & \ldots & \mathbf{x}_8\end{bmatrix}$$
where $\mathbf x_k$ is the $k$ column of $x$.
I define the vertically stacked vector $\mathbf z$:
$$\mathbf z = \begin{bmatrix}\mathbf{x}_1 \newline
\mathbf{x}_2\newline
\vdots\newline
\mathbf{x}_8\end{bmatrix}$$
Now we can perform your column-wise IDFT of $x$ as:
$$\mathbf Z = \begin{bmatrix}\mathbf F^H \newline
& \mathbf F^H \newline
&& \ddots \newline
&&& \mathbf F^H\end{bmatrix} \mathbf z = (\mathbf I_8 \otimes \mathbf F) \mathbf z$$
$\mathbf Z$ is the column-stacked equivalent of your $X$. Reading $X$ row-wise corresponds to reading $\mathbf Z$ through a permutation:
$$\mathbf Z' = \mathbf P \mathbf Z$$
where
$$\mathbf P = \begin{bmatrix}1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \ldots\newline
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & \ldots\newline
\vdots\newline
0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 &\ldots\newline
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & \ldots\newline
\vdots\newline
0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & \ldots\newline
\vdots\newline
&&&&&&&&&& \ldots & 0 & 1\end{bmatrix}_{(64 \times 64)}$$
I hope the pattern is understandable: the first 8 rows of $\mathbf P$ pick out the 1st, 9th, ... 57th entry of $\mathbf Z$, the next 8 rows of $\mathbf P$ pick out the 2nd, 10th, ... 58th entry. The last 8 rows of $\mathbf P$ pick out the 8th, 16th, ... 64th entry of $\mathbf Z$.
Now $\mathbf Z' = X'$ and we can write:
$$y = H \mathbf P (\mathbf I_8 \otimes \mathbf F) \operatorname{stack}(x)$$
Notice how I mix your notation with mine, so sorry it looks a bit messy. The "stack" operator stacks your matrix $x$ column-wise into the vector $\mathbf Z$.

It is important to notice here that this is mostly a theoretical exercise to see that the problem can indeed be brought on a standard form, familiar in compressed sensing, that can be plugged into for example OMP. You do not necessarily have to implement it this way. You may be able to still use the operations the way you describe them, and in particular it is more computationally efficient to compute the column-wise IFFT of your $x$ rather than forming the actual DFT Kronecker product matrix $\mathbf I_8 \otimes \mathbf F$ that I describe. However, it may be easier to implement by explicitly using my formulas here.
To put it another way, if you have a sparse reconstruction algorithm that understands the general compressed sensing problem $\mathbf y = \mathbf A \mathbf x$, then you can form $\mathbf A = H \mathbf P (\mathbf I_8 \otimes \mathbf F)$ and $\mathbf x = \operatorname{stack}(x)$ and plug those into the reconstruction algorithm. Whether you we will be able to "deconstruct" $\mathbf A$ and retain a more computationally efficient structure depends on your reconstruction algorithm and its interface.


BPSK vs $\pi/2$ BPSK, what the benefits?

The signal always changes phase at bit transition times. This makes it easier to achieve bit synchronism at the receiver. In plain BPSK, long runs of 0 or 1, with no phase transitions occurring at the boundary
between two identical bits, the phase-lock loop that is generating the
signal telling the receiver where the bit transitions are, can lose lock; or the initial acquisition is more difficult. Also, in satellite channels with nonlinear amplifiers such as traveling-wave-tube (TWT)
amplifiers, smaller phase transitions are better. In particular, if the RF carrier phase makes
a transition of $\pi$ radians (meaning the carrier envelope makes a
transition through the origin), the TWT amplifier output contains
all the sidebands that were carefully suppressed in the uplink. So,
filtering in the satellite is required to suppress these sidebands.
Hence, $\pi/2$-BPSK and $\pi/4$-QPSK are preferred over their plain cousins. Note 4, not 2, with QPSK: $\pi/2$-QPSK does not avoid the
transitions through $0$, and adds complexity without necessarily
adding benefits.

Note: in plain BPSK, the signals are $\pm\cos(\omega_0 t)
= \{\cos(\omega_0 t), \cos(\omega_0t + \pi)\}$, and so,
at bit transitions, the phase either stays the same, or changes by
$\pi$. in $\pi/2$-BPSK, the signal set is $\pm\cos(\omega_0 t)$
for one bit and $\pm\sin(\omega_0 t)$ for the next. Thus,
regardless of whether we have two successive identical bits
or two different successive bit, the signal phase changes by $\pm \pi/2$
since the signal changes from a cosine to a sine.  Figure out for
yourself why $\pi/2$-QPSK has transitions through $0$ while
$\pi/4$-QPSK does not; then think of offset QPSK a little.


Your call to fir1 looks OK, but you should check the result yourself. You could probably answer several of your other questions by carefully reading the mathworks documentation of the respective functions. E.g., for fir2 you can read that

  f is a vector of frequency points in the range from 0 to 1, where 1 corresponds to the Nyquist frequency. The first point of f must be 0 and the last point must be 1. f must be sorted in increasing order. Duplicate frequency points are allowed and are treated as steps in the frequency response.

and

  m is a vector containing the desired magnitude response at each of the points specified in f.

So you shouldn't define m as a complex vector; the routine knows the phase lag because it's just half of the filter order. The vector f should be defined as follows:

f = [0 .6 .6 .8 .8 1]

Since this is homework, I leave the definition of the magnitude vector up to you.

Your call to butter looks OK, apart from the fact that you used stop as the filter type. You want a band pass filter, not a band stop filter, so just leave out the additional filter type argument, and the routine will by default design a band pass filter (because there are two edge frequencies in the frequency vector).

For the impulse invariance method, you first have to design an analog Butterworth filter (using butter with last input argument s for analog filter), and then use impinvar to create the corresponding digital filter. The mathworks documentation is very clear, and I think you will succeed after having studied it. Otherwise, add to your original question a very specific question where you're stuck.


I can't make sense of the amplitude spectrum of a signal. Help please

A few suggestions and comments:

You calculate the frequency vector with F = [-N/2:N/2-1]/N;. This is wrong. You need to take the sampling frequency into account.
The way you organized your code, you don't have direct control over the sampling frequency. I recommend specifying both the time length and the sampling frequency explicitly in your code.
You need to set up your analysis so that you get high resolution in the frequency domain (so that you can resolve peaks that are very near to each other). The only way to achieve this is to sample for a long time. The more you oversample the signal (that is, the higher the Nyquist frequency is), the longer time you need to observe the signal. This is the reason you want to have explicit control over sampling rate and signal duration.
Do not limit the FFT size by choosing an explicit N, since that may lower your frequency resolution. If you want to set the FFT to a power of 2, then use nextpow2 or similar.


Damped spring mass system - parameter estimation

I think I got the answer: I look at the step response of the system.

$H(s) = \frac{1}{ms^2+ds+c}$

The limit of the step response equals my parameter $1/c$

$y(\infty)=\lim_{s\rightarrow 0} s\cdot \frac{1}{s} \cdot H(s) = \frac{1}{c}$

Using $c$ and $c/m$ from the measurements above I can calculate the mass $m$.

I hope this is right, but the results look reasonable...


Simulating analog filter using convolution or converting with fft

What's not clear to me is what the fundamental difference (if any) is between simulating an analog filter and making a digital filter.

Either way, these functions will produce "ba" transfer function outputs, but the b and a are totally different.

For a 2nd-order filter, for instance, b = [b0, b1, b2] and a = [a0, a1, a2].  These are the coefficients of the transfer function.

For an analog filter, this represents a transfer function like this:

$$H(s) = \frac{b_2 s^2 + b_1 s + b_0}{a_2 s^2 + a_1 s + a_0}$$

For a digital filter, the transfer function is:

$$H(z)=\frac{b_0+b_1z^{-1}+b_2z^{-2}} {a_0+a_1z^{-1}+a_2z^{-2} }$$

  Note that the input signal is an array of values

Then you are doing digital filtering.  You can simulate analog filtering this way, but it's only going to be a simulation.

All the Python IIR filter functions you're talking about, when outputting a digital filter, design it as an approximation of the analog filter using an analog filter transfer function prototype and then converting it to digital using the bilinear transform and frequency warping.  Yes, they can output analog transfer functions, but I think that's only useful for building actually analog electronic filters, or for plotting the frequency response, impulse response, etc.

The Bessel in particular is a bad approximation as you get near fs/4, because the important property of a Bessel filter is group delay, not amplitude, and phase/group delay is not preserved by the bilinear transform (red vertical line is fs/4):

More plots of digital Bessel filter

You can make the simulation closer to what an analog filter would produce by putting a large margin between the highest frequency in your signal and the sampling frequency, called oversampling.  Assuming your "array of values" is correctly bandlimited, then you can upsample them and then apply the digital filter using lfilter and the ba digital filter coefficients you found earlier.  If not, you need to bandlimit the analog signal and then sample it at a significantly higher rate.  How did you get your sampled signal?  If you are generating it digitally, it is probably not bandlimited, unless you specifically thought of bandlimiting while generating it.

  for example a 1st-order Bessel filter

All 1st-order filters are the same, the "1st-order Bessel" is the same thing as a "1st-order Butterworth" or anything else.


Corner detection using Chris Harris & Mike Stephens

Harris Corner detector tries to quantify the local intensity changes at all the directions for each pixel. The figure below illustrates the basic idea clearly:

So $I(x+u,y+v)$ indicates the pixel intensities of all the neighborhood pixels around $(x,y)$.
The window function is applied for feature localization. For most often used Gaussian function, the choice of sigma value corresponds to different feature width in the image.


Why would I want to define a modulation index for each tone (DSB-FC)?

I'm no expert on this, but until the big guys give an answer, here's my take:

For $f(t) = \cos(100t)+2\cos(2\cdot 100t)$, so $\omega_m = 100\ \rm Hz$ for example:

The amplitude's peak value$\ =3$ exceeds $A=2$, so there is overmodulation. This answer only addresses your last point of "I'm pretty sure there's no overmodulation".


Could multiple echo reflections be combined to improve signal to noise ratio?

You need to be careful to look at your bunch of impulses in the fequency domain, to make sure you are really transmitting more energy at all frequencies. Your chain of 3 identical impulses will contain zero energy at certain frequencies. As a result, your echo-separation method (Tk = Sk – (Tk-10) – (Tk-20)) will have infinite gain at those frequencies (this is an IIR filter with 2 poles on the unit-circle), and hence it will be unstable.

As an example, I built 4 different "pulse train" signals (in Matlab), by the following method (an I've plotted the magnitude of the FFT):

>> S1 = [1;0;0;0;0;0;0;0;0;0;0];
>> S2 = [1;0;0;0;0;1;0;0;0;0;0];
>> S3 = [1;0;0;0;0;1;0;0;0;0;1];
>> S4 = [1;0;0;0;0;1;0;0;0;0;-1];
>> plot( abs([fft(S1,256),fft(S2,256),fft(S3,256),fft(S4,256)]) );

Note that S1 has a constant magnitude respnnse (equal to 1), as expected. S2 (made up of 2 impulses) has a magnitude response that goes to 0 at several frequencies (even though, on average, it has twice the power of S1).

S3 has three times the average power as S1, but unfortunately it also has zero energy at some frequencies.

S4 is built with a -1 in the last sample. As a result, it still has 3 times the average energy compared to S1, with the added property that is it never less than unity energy in any one frequency band.

The (simplistic) recovery-filter, for use with S4, would be :  Tk = Sk – (Tk-5) + (Tk-10)), but unfortunately this is also an unstable filter (having some poles outside the unit circle). The trick to making a filter like this stable is split it into two parts.

First, you run a regular IIR filter over the signal as follows:

$$T'_k = S_k + g \times T'_{k-5}$$

where $g=\frac{\sqrt5-1}{2}$. Next, you run this filter:

$$T_k = T'_k - g \times T_{k+5}$$

and the trick, with this second filter, is to compute the samples of $T_k$ in reverse order (if you want to compute 100 samples, start with sample 99, and work your way down to sample 0).

The overarching principles here are:

Spread your energy out over multiple samples, in the transmitted signal, to allow you to send more energy
Try to always ensure there is lots of energy in all frequencies of the transmitted signal (look at the FFT of the signal)
Build an inverse filter

My inverse filter method shown here is a bit perverse, using a time-reversed IIR filter. If you have plenty of CPU cycles, you can compute the correct inverse filter by building it as a long FIR filter (in theory, you need an infinitely long filter, but a long enough FIR will approximate it OK).


Justification for Squared $ {L}_{2} $ Data and Smoothness Term as an Error Bound

One of the motivations to use the $ {L}_{2} $ norm comes from the Maximum a Posteriori Estimation (MAP) framework.

If you model $ \psi \left( u \right) \sim \mathcal{N} \left( 0, \alpha \right) $ then if you derive the MAP Estimator in case the added noise is Gaussian you'd get the exact model you posted above.

An example of the derivation of MAP model to the above can be seen in my answer to Estimating the Signal by Deconvolution with a Prior on the Filter Coefficients and the Signal Samples.


If you want to do it yourself, programming a linear interpolation is really simple. You interpolate between sample 1 and sample 2 by calculating the slope to generate a line between the samples (y are your measurements and x the times)

\begin{align}
m     &= \frac{y_1 - y_2}{x_2 - x_1} \\
y_1(x) &= m (x - x_1) + y_1
\end{align}

this is the linear interpolation between sample 1 and 2, it allows you to calculate the interpolated values for any point in time between sample 1 and 2. After that you just have to repeat this for sample 2 and 3 and so on. It's always the same, that means it comes down to writing one function and a loop. Just google linear interpolation and you can probably find some nice animations and youtube videos that visualize the process.

If you don't want to write the code yourself Scipy offers you multiple ways of interpolating data you can find a nice documentation here:

https://docs.scipy.org/doc/scipy/reference/tutorial/interpolate.html

the first way is probably the simplest (i just stole the examples from the page there):

from scipy.interpolate import interp1d

x = % array with the timestamps of your measurements
y = % your measurements

f1 = interp1d(x, y) % default is linear
f2 = interp1d(x, y, kind='cubic')

% xnew contains all points at which you want to sample the interpolated function
% linspace should go from your start time to your end time
xnew = np.linspace(0, 10, num=1001, endpoint=True)
ynew1 = f1(xnew) % generate the y values for all x values in xnew
ynew2 = f2(xnew)

now you can just plot(xnew, ynew) and it should look interpolated.

spline interpolation is another popular choice you can find it a bit further down with another example you could try that too if you'd like to. I personally like spline the most but none of them will magically generate data that means there really isn't a magical solution that is best, it depends on what you want to do with the result.


A system is memoryless if its output ($y(t)$) for each value of the independent variable ($t$ in this case) at a given time is dependent only on the input at that same time ($x(t)$).

Every system that consists of a delay (like the one in your example) or an accumulator, for example, are systems with memory.

This can be seen just by replacing $t$ by some number. Take, for instance, $t=0$. In that case, the output for that instant ($y(0)$) depends on the input for another value of the independent variable ($x(C),C>0$); that means that the system has memory.


If you are in a position to acquire $N$ snapshots of data, $s_n(t)$, for $n=0\ldots N-1$ then you can just do:

$$
S_n(k) = |\mathrm{fft}(s_n(t))|^2\\
\hat{S}(k) =  \frac{1}{N} \sum_{n=0}^{N-1} S_n(k)
$$

and a certain amount of averaging will reduce the apparent noise on the spectrum.


Fast & accurate convolution algorithm (like FFT) for high dynamic range?

Disclaimer:
I know this topic is older, but if one is looking for "fast accurate convolution high dynamic range" or similar this is one of the first of only a few decent results. I wanna share my insights I got on this topic so it might help somebody in the future.
I apologize if I might use the wrong terms in my answer, but everything I found on this topic is rather vague and lead to confusing even in this thread. I hope the reader will understand anyway.

Direct convolution is mostly accurate to machine precision for each point, i.e. the relative error is usually roughly or close to 1.e-16 for double precision for each point of the result. Each point has 16 correct digits. Rounding errors can significant for untypically large convolutions though, and strictly speaking one should be careful with cancellation and use something like Kahan summation and high enough precision data types, but in practice the error as almost always optimal.

The error of a FFT convolution apart from rounding errors is "global relative" error, meaning the error in each point depends on the machine precision and the peak value of the result. For example if the peak value of the result is 2.e9, then the absolute error in each point is $2\cdot10^9\cdot10^{-16} = 2\cdot10^{-7}$. So if a value in the result is supposed to be very small, let's say $10^{-9}$, the relative error in that point can be huge. FFT convolution is basically useless if you need small relative errors in the tail of your result, e.g. you have a somewhat exponential decay of your data and need accurate values in the tail. Interestingly if FFT convolution is not limited by that error, it has much smaller rounding errors compared to direct convolution, since you obviously do less additions/multiplications. This is actually why people often claim FFT convolution is more accurate, and they are almost right in some sense, so they can be quite adamant.

Unforunately there is no easy universal fix to get fast and accurate convolutions, but depending on your problem there might be one... I have found two:

If you have smooth kernels which can be approximated well by a polynomial in the tail, then the black-box Fast Multipole Method with Chebyshev interpolation might be interesting for you. If your kernel is "nice" this works actually perfectly: you get both linear (!) computational complexity and machine precision accuracy. If this fits your problem you should use it. It's not easy to implement however.

For some specific kernels (convex functions I think, usually from probability densities) you can use an "exponential shift" to get optimal error in some part of the tail of the result. There is a PHD thesis and a github with a python implementation using that systematically, and the author cals it accurate FFT convolution. In most cases this is not super useful however, since either it regresses back to direct convolution or you can use FFT convolution anyway. Although the code does it automatically, which is nice of course.

--------------------EDIT:--------------------

I looked a little bit at the Karatsuba algorithm (I actually made a small implementation), and to me it looks like it has usually similar error behavior like the FFT convolution, i.e. you get an error relative to the peak value of the result. Due to the divide and conquer nature of the algorithm some values in the tail of the result actually have better error, but I don't see an easy systematic way to tell which ones or in any case how to use this observation. Too bad, at first I thought Karatsuba might be something useful in-between direct and FFT convolution. But I don't see common use cases where Karatsuba should be prefered over the common two convolution algorithms.

And to add to the exponential shift I mentioned above: There are many cases where you can use it to improve the result of a convolution, but again it's not a universal fix. I actually use this together with FFT convolution to get pretty good results (in the general case for all inputs: at the worst same error as normal FFT convolution, at best relative error in each point to machine precision). But again, this only really works nicely for specific kernels and data, but for me both kernel and data or somewhat exponential in decay.


Fourier transform of Autocorrelation function, what am I missing here? (Power spectrum)

That question was posted quite some time ago! Have you tried the discrete time Fourier transform of $R_{XX}[m]$, instead of the Fourier transform? I believe this is what was meant, there (the text talks about a "sequence").

An "easy" way to solve this, according to me, is to compute the Z-transform, then identify $S_{XX}(\omega) = Z[R_{XX}](z=e^{j\omega})$. The Z-transform is given as:
\begin{align}
Z[R_{XX}](z) &= \sum_{m=-\infty}^{\infty} R_{XX}[m] z^{-m} \\
   &= \sum_m e^{-\alpha |m|} z^{-m} \\
   &= \frac{1 - e^{-2\alpha}}{(1-e^{-\alpha}z^{-1})(1 - e^{-\alpha}z)}
\end{align}

Hope that still helps!


I would assume by faults you mean breaks in a signal line such that a reflection occurs. This is indeed an application of the autocorrelation: You transmit a sequence down a transmission line. If there is any change in impendance in the the line (such a break or kink etc) then a portion of your signal will be reflected back according to the reflection coefficient, given as:

$$\rho = \frac{Z_L-Z_o}{Z_L+Z_o}$$

Where $Z_L$ is the impedance of the load and $Z_o$ is the impedance of the transmission line, and $\rho$ is a complex valued reflection coefficient with magnitude ranging from -1 to +1, and is the amount of signal reflected back toward the source at the load (from 0 to 100% and with any phase angle).

Before explaining the autocorrelation further and how you would use it, I think it would be important to know some basic transmission line theory and in particular the properties that are given by this reflection coefficient.

For example, if we had a transmission line that had an impedance of 50 ohms, and if it was terminated with a 50 ohm resistor at the load, the reflection at the load would be 0 as given by the numerator of the equation, meaning none of the signal was reflected (this is the maximum power transfer condition when we match impedance of load to the line and to the source).  In such a condition, if we looked into the properly terminated transmission line from the source, since no reflection comes back it would look no different to us than if we were looking directly at a 50 ohm resistor at the source (this is one way to define the impedance of the cable; an infinitely long cable will have the same impedance as a resistor of the same value-- a shorter cable will only also look like this if it is terminated with a resistor having the value of the cable's impedance.)

Observe these other interesting conditions about the reflection coefficient and see how much it can tell us about the load.

$$Z_L = Z_o \rightarrow \rho = 0$$

$$Z_L > Z_o \rightarrow 0 > \rho \ge 1$$

$$Z_L < Z_o \rightarrow -1 \le \rho < 0$$

If the load is an open, the reflection is 1 (100%) and in phase.

If the load is a short, the reflection is -1 and out of phase.

As the reflection propagates back down the line, it constructively and destructively adds with the forward signal from the source as the line itself has phase shift proportional to the position on the line and the frequency of the signal. (A constant delay has a negative linear phase versus frequency). The velocity of this reflection is the speed of light divided by the square root of the dielectric constant of the transmission line (that sentence may be confusing, bottom line the reflection goes at the speed of light, or slower if in a transmission line with a higher dielectric constant than air).  So we can measure the reflection at the source, but we can't yet tell if the load (or anywhere along the line where a discontinuity occurs, essentially becoming our new "load") was an open or a short (if it was one of the two). If the length of the line was exactly a quarter wave-length in propagation distance of our source signal's frequency, then the source signal will have shifted 90° by the time it reached the load, and any reflection will shift yet another 90° by the time it reaches our source for a total 180° phase shift. This is exactly how a quarter wave line translates the impedance of a source to an open. Similarly a quarter wave 50 ohm line will translate the impedance of a 100 ohm load to 25 ohms (these are called impedance transformers and we use these techniques to match different impedances to minimize reflections that otherwise distort our signals).

With all that said, the great thing we can do with the autocorrelation of special sequences is to accurately resolve our distance in propogation time to the mismatched load or fault (a cable RADAR). As far as special sequences, we need to use sequences that have desirable auto-correlation properties, ideally ones that have a strong correlation to themselves only when they are completely aligned and zero at all other time offsets. White noise has this property specifically, so that is what we seek: good approximations to white noise. Such sequences are called pseudo-random sequences given they appear random yet we know exactly what the sequence is (and we need to in order to correlate to it).  If a sequence didn't have this property, meaning it also correlated to delayed versions of itself, then this would reduce our ability to accurately measure range to a reflection given multiple delays will also create a response from the one reflection.

Thus to do this, you transmit a psuedo-random sequence down the transmission line, while measuring the reflection (a "directional coupler" is a useful low cost relatively simple device for doing this). Correlate the reflection to the sequence at all delay offsets. The resulting plot of correlation vs delay offset will indicate faults AND range to fault based on size of peaks and distance from the origin in this correlation function. Opens and shorts will have the strongest peak. If you don't have an open and short but a damaging kink such as to change the line impedance (which would cause signal distortion and other problems so is a worthy fault to find) then you will see a smaller reflection at the range of this fault along with another likely even smaller reflection from the load that isn't perfectly matched. With a sensitive measurement you will often see something from the load due to imperfect matches at the physical transition from the transmission line to the load. If your source also is not well matched to the line, then you can see multiple reflections as a strong reflection bounces back and forth between source and load, but each of these will be at progressively further offsets from the origin with progressively smaller magnitudes so it will be very clear what is going on, especially if we have an estimate of the actual length of the transmission line.

I don't really understand Question 3. Autocorrelation is indeed simply multiplies and accumulation applied to signal samples, and is indeed a statistical approach (not either/or). One very important aspect of this is when you add samples that have noise values given by an independent identically distributed random process (the background noise you will get when you receive your reflected signal from the transmission line), the standard deviation of the sum will go up at $\sqrt{N}$ where $N$ is the number of samples you are adding. Coherent samples (which is what occurs when the sequence you are testing for is in perfect time alignment with the received reflected signal) will go up at $N$, so thus you achieve a $\sqrt{N}$ increase in the signal to noise magnitude ratio (in dB to be clear $10Log10(N)$).

Consider the simple Barker Code sequence that does have this nice autocorrelation property we seek, and assume you transmit a +1 when the code is 1 and an -1 when the code is 0. When you multiply this with itself in time alignment and accumulate the result, it will add to 11. The standard deviation for the IID random process of any noise on each sample will however increase by only $\sqrt{11}$. I leave it as an exercise for you to to see how it only adds to 1 for any other rotational shift. Consider what would occur with practical sequences that can have much longer lengths given by $2^{k}-1$ where k is a positive integer.

$$\begin{bmatrix} 1 & 0 & 1 & 1 & 0 & 1 & 1 & 1 & 0 & 0 & 0 \end{bmatrix}$$

To answer your last question, you can do all this with an FFT since the circular autocorrelation can be computed using the FFT:

$$XCORR = ifft(fft(a)(fft(b^*)))$$

Where $a$ is the perfect sequence as transmitted and $b$ is the reflected complex signal as received. The * represents the complex conjugate. The FFT is used but you would still be doing the process I described of using the autocorrelation property of the sequence. If you had another FFT approach that could give you the range and intensity of the fault please clarify what that is and we can further compare if you don't see it now yourself.


You are seeing the "Digital IF” carrier assuming you are modelling what you want to implement digitally in VHDL. The sine and cosine blocks together with the multipliers and combiner perform a frequency translation from DC to the carrier mentioned.  Thus the constellation at this point is really not observable as it will be "spinning" at this carrier rate. (The constellation as observed in the upper figure is at DC and to note it appears to be sampled at 2 samples per symbol, which is typical, so it is really every other sample (the ones on the diagonals) that represents the QPSK constellation at the correct symbol sampling locations).

To observe the constellation you would need to frequency shift the signal back to DC while being careful not to introduce phase and frequency offsets, or shifts in sample timing (a receiver does this with carrier and timing recovery loops so that it can run completely on an asynchronous clock from the transmitter).

(To note, an optional approach is "Zero-IF" where you implement the frequency translation in the analog domain using an IQ mixer (which is the same functional block with analog components)


Discrete Wavelet Transform (DWT) Filter Bank

In the DWT scheme, whether it is the classical $2$-band or the $M$-band wavelet setting, the very same analysis filter bank (lowpass/highpass + subsampling) is used at each level. Under this condition, one can derive the cascade algorithm that provides  the spectrum of the scaling function:
$$\Phi(\omega)= \prod_{k=1}^\infty \frac {1} {\sqrt 2} H\left( \frac {\omega} {2^k}\right) \Phi^{(\infty)}(0)$$
where iterated half-cut-off frequencies are apparent.

However, this is especially important when one wants to address the properties of the underlying wavelet, or cascade many wavelet levels. In practice, you can easily choose a different set of perfect reconstruction filter bank at each level, for several reasons:

actual filters down the scales tend to become larger, sometimes too much with respect to the signal features, and the signal length.  When a convolutive filter becomes longer than the data itself, you run into troubles, especially in regard of boundary handling (symmetry, padding). The consequence can be high details coefficients caused mostly by wrapping around border artifacts.
filtered/down-sampled data don't have the same characteristics at each scale, so different filters can do a better job
do we really need multicale theory after all on 3 levels? In practice, it can be useful to use longer filters on the first levels, for better frequency separation, and shorter ones on lower levels, for the multiscale effect.

This is not mainstream, as one has to choose filters wisely, sometimes without many hint. However, there were works with DWT or GenLOT filters on the first scales and DCT or Haar on the lower resolutions.  JPEG XR standard somehow behaves in that philosophy. One of the papers I have met on this direction is: On the frame bounds of iterated filter banks, 2009.


Relation between Covariance matrix & Energy of a random signal

I'm not totally sure I understand your question: could you give a concrete example? For instance, what's the randomness? Anyway, here are some relations between the quantities you're interested in.

The energy of the signal is $$\mid \mid Y \mid \mid_2^2 = \sum_{i=0}^n y_i^2$$

The Auto-correllation matrix $R$ is formed as $YY^T$ (outer product between the signal and it's transpose).

Since this matrix is symmetric, it's matrix 2-norm is the spectral radius: $$||YY^T||_2^2 = \mathrm{max}|\lambda_j|$$

Another way to write the 2-norm would be as $$||YY^T||_F^2 = Tr(YY^T)$$ It's fairly stragtforward to prove that $$||* ||_F^2 \geq ||*||_2^2$$ for any matrix, and they co-incide when $YY^T$ has rank 1. It's also true that $$|| * ||_F^2 \leq r|| *|| _2^2 $$


Calculating incoherency in compressive sensing

As a starting-point, you might try something like this:

Step 1: create random (complex-valued) matrices

rows = 32;   %number of output samples from compression matrix
cols = 2048;  %number of input samples supplied to compression matrix
matrixA = randn(rows, cols) + 1i*randn(rows, cols);
matrixB = randn(rows, cols)+ 1i*randn(rows, cols);

Step 2: Compute Maximum Dot-Product Magnitude

%compute all possible cross-correlations
xc = abs(matrixA * matrixB');
%find the maximum
result = max(max(abs(matrixA * matrixB')));
%todo:  scale by sqrt(n)
%...

Step 3: Beware of Amplitude Scaling Impacts

Since the coherence equation that you provided depends upon dot-products, the amplitude scaling of your matrix will effect the result.  Most likely the equation you provided assumes the matrices are scaled in amplitude.  For example, each column in matrixA & matrixB might need to be scaled by it's "Euclidean Norm" or "L2-Norm".  You might find it helpful to read about correlation coefficients which are scaled accordingly, or revisit your equation reference just to be safe.


I don't really think there's one:
You seem to be in the process of figuring out what about a piece of audio signifies stress.
As long as you don't know that, you can't erase any identifying information from it without standing a high probability of also erasing your "stress" information.


Deriving 2-D discrete Fourier transforms

Here is the first one:

By definition,

$\mathfrak{F}\{ f(x,-2y)\} = \iint_{- \infty}^{\infty} f(x,-2y)e^{-j2\pi (ux+vy)}dxdy$

let $\tau = -2y$ and conversely $y=\frac{-\tau}{2}$

$\mathfrak{F}\{ f(x,-2y)\} = \iint_{- \infty}^{\infty} f(x,\tau)e^{-j2\pi (ux-\frac{v\tau}{2})}dxd(-\frac{\tau}{2})$

$\mathfrak{F}\{ f(x,-2y)\} = -\frac{1}{2}\iint_{- \infty}^{\infty} f(x,\tau)e^{-j2\pi (ux-\frac{v\tau}{2})}dxd\tau$

$\mathfrak{F}\{ f(x,-2y)\} = -\frac{1}{2} F(u,-\frac{v}{2})$

The other one will be harder, but I will leave it to you.


Modulation of data bits into complex symbols in MATLAB

The method step(H,X) used here has nothing to do with the step function. It modulates the input data X with the modulator object H, as described here. See also here.


If you are looking to get more mfcc coefficients for the given signal, you can do that by using librosa.feature.mfcc(y=y, sr=sr, n_mfcc=40), here you will get 40 coefficients for each time component.


Reducing a signal to Acos(omega * t+ phi) with negative omega*t

You flipped a sign on the phase, the answer is -.529....
It's not
$$X(t) = Re(2e^{j(\pi/4)}e^{j\omega_0 t})+Re(e^{j0}e^{j\omega_0 t})$$

It's
$$X(t) = Re(2e^{-j(\pi/4)}e^{j\omega_0 t})+Re(e^{j0}e^{j\omega_0 t})$$


When I interview people, my questions depend a lot on the CV they've submitted.

If they're fresh from university, I am much more interested in whether they have ever heard of some of the basic image processing concepts: What is their definition of noise and background? What does a Gaussian filter do? And of course, I want to know if they know how to program as opposed to hacking together a last-minute solution for homework.

If they claim to be more experienced, I want to know how well they know what they're talking about. What's their opinion on supervised vs. unsupervised segmentation? What are good tracking algorithms, and why? What informs their choice of a filtering algorithm? How do they deal with stochastic data? And of course, do they know how to program?

Finally, I'm always happy to find out that they did some research into the topic ahead of time, even if they just read my papers and totally misunderstood them. I can always teach them knowledge, but it is hard to teach effort.


How to understand the filter coefficients in FIR filter

2-tap FIR filters can be understood with a simple intuition (I use Python for plotting, but freqz is the same there).

Let's start with two taps with the same value: b = [1, 1]. This lets through constant signals (e.g. [1, 1, 1, ...]) but will more or less block rapidly changing ones (e.g. [1, -1, 1, -1, ...]), i.e. it is a lowpass:

w, h = scipy.signal.freqz([1, 1], 1)

plt.plot(w, abs(h))

On the other hand, b = [-1, 1] will do the opposite: let changing signals pass but block constant ones, i.e., we have a highpass.

w, h = scipy.signal.freqz([-1, 1], 1)

plt.plot(w, abs(h))

I'm not sure it is straightforward to expand this intuition to more complex signals like yours. But looking at the shape, I would expect some band-pass behavior (because we have both high- and lowpass components) with a peak frequency around one third of the sampling frequency (because a signal like [1, 0, 0, -1, 0, 0, 1, ...] will be able to pass). This apparently isn't so far from the truth:


A square window seems to be another name for a rectangular window.  A rectangular window is the default window used to apply an FFT to a subset of a longer stream of data and without any additional scaling or weighting of any of the subset's samples.


How does constant frequency response relate to equiripple?

What your teacher tried to show is that a transfer function which is the ratio of two polynomials can never be constant (excluding the trivial case where both polynomials are of degree zero). This means that if your desired frequency response is a piecewise constant function (as is the case for an ideal frequency selective filter), you can only approximate that piecewise constant function by a realizable transfer function, and one popular approximation is an equi-ripple approximation, where the maximum deviation from the desired constant value is minimized.


Do high frequency oscillations use more than one bandwidth

A modulated radio signal has non-zero bandwidth.  The modulation signal that carries the information has a much lower max frequency compared with the carrier frequency, keeping the total bandwidth (carrier plus both modulation sidebands) of an AM signal within FCC legal limits.

Other modulation schemes are used for GSM/3G, such as OFDM, that create more widely spaced sidebands, but the carrier frequency is so much higher (GHz range) than the modulation sideband width, that the resulting signal bandwidth is only a tiny fractional slice of the total RF spectrum.


FFT and DFT are a faster way to implement convolution.
Since convolution in time domain corresponds to multiplication in the Fourier domain, you replace the sliding-window based implementation of convolution (in time domain) by a pointwise multiplication (in Fourier domain). This avoids you an extra inner loop over the input data and saves time.

With respect to boundary conditions, there is a minor difference with finite-length signals.
Using the convolution theorem to compute a convolution product in the Fourier domain implicitely assumes that the input data is periodic, i.e., when you reach the right end of the signal you re-enter by its leftmost part. This is also known as circular convolution.
If you need to have a convolution of explicitely non-periodic signals, you can add zeros at the beginning and end of your signal to avoid these border effects.


MATLAB: Implementing Least Squares Estimator for a Given Model

The equation you're trying to solve is
$$
\mathbf{y}=\mathbf{X}\mathbf{h},
$$
where $\mathbf{h}$ is your unknown. The matrix $\mathbf{X}$ is going to have a time-shifted structure that reflects the convolution operator. If we assume that the $\mathbf{y}$ vector starts with y(3) i.e. ignores the first two zeroed out elements of y,  then the corresponding $\mathbf{X}$ matrix is given by:

$$
\mathbf{X}=\left[\begin{matrix}
x(3) & x(2) & x(1)\\
x(4) & x(3) & x(2) \\
x(5) & x(4) & x(3) \\
... & ... & ...
\end{matrix}  \right]
$$

You can add as many rows as you are have observations to support them. You can then solve for $\mathbf{h}$ by
$$
\mathbf{h}_{est}=\text{pinv}(\mathbf{X})\mathbf{y}
$$


References and resources on hyperspectral (multispectral) image processing

What kind of information are you after.  Are you looking for information on how to process the  images from a data processing point of view?

Or from a theory point of view?

Are you looking at satellite images?  if so land or water?

Proprietary software look at ENVI

http://www.exelisvis.com/language/en-us/productsservices/envi.aspx

for free software look at SeaDAS

http://oceancolor.gsfc.nasa.gov/seadas/


LPF in the stage of IQ demodulator is it a analgor filter or digital filter?

The A/D can be placed as a single real A/D before the multipliers, OR as shown in the diagram as two A/Ds one after each multiplier to sample the I and Q channels. In either case, an analog filter is required before any A/D conversion as an anti-alias filter. This can be a bandpass filter or a low-pass filter, depending on which image in the analog domain's frequency spectrum is desired as the signal of interest. Even if there is no actual additional signal in the image bands, with no filter (and gain in the receiver) the noise figure will be significantly degraded due to the accumulated noise floor from each of the images. An analog filter prior to any sampling is important.
So in this case given the placement of the A/D converter, the filters as shown would be the analog anti-alias filter that I describe.
If the reason for image rejection is confusing, please refer to this posts that details how aliasing occurs in the digital sampling process: Higher order harmonics during sampling
So in summary, for the implementation of an IQ demodulator, there could be two A/D converters if I and Q are to be sampled separately: one for the I channel and one for the Q channel resulting in a complex representation of the baseband signal. Alternatively there can be a single A/D converter prior to the multipliers and then the multiplier and filter structures as shown can be implemented all digitally (referred to as a Digital Down-Converter or DDC). In this case, a single analog anti-alias filter (band-pass or low pass as desired for the implementation) is still needed in front of the A/D converter, in addition to the digital low pass filters after the multipliers to select the low frequency component after the multipliers.
The block diagram as shown in the question does not make that clear that $v(t)$, $y(t)$ and $\hat d[k]$ are complex signals (although clearer if you read through the linked reference).


How are phase values able to capture motion from video?

If you move an impulse in the spatial domain (perhaps representing a pixel on some object in a scene), the unwrapped phase of all the spectra associated with that point will be rotated in the frequency domain by an amount proportional to the distance of the move and the frequency.


Basic question about trigonometric series and transforms thereof

There is no general relationship between the Fourier transform of $f$ and that of $g(f)$ where $g$ is an arbitrary function. The Fourier transform does have the linearity property, so if $g$ is something simple like an affine transform, then the same linear relationship applies to their transforms $F$ and $G$.

With respect to your second question, where $h = \sum_k \left(A_k \sin(\omega_k x)\right)^N$, the presence of more than $k$ peaks in the Fourier transform of $h$ is easily explained using the power-reduction trigonometric identity:

$$
\sin^n(\theta) =
\begin{cases}
\frac{2}{2^n} \sum_{k=0}^{\frac{n-1}{2}} (-1)^{(\frac{n-1}{2}-k)} \binom{n}{k} \sin{((n-2k)\theta)}, & n \text{ is odd} \\
\frac{1}{2^n} \binom{n}{\frac{n}{2}} + \frac{2}{2^n} \sum_{k=0}^{\frac{n}{2}-1} (-1)^{(\frac{n}{2}-k)} \binom{n}{k} \cos{((n-2k)\theta)}, & n \text{ is even}
\end{cases}
$$

(the above is shamelessly borrowed from Wikipedia)

So, when you raise a sinusoid to a power, the result can be expressed as a weighted sum of sinusoids at different frequencies, where the number of individual terms is related to the power. That's why you see additional peaks in the spectrum of $h$.

You can come up with a more general relationship for some cases by taking advantage of the multiplication property of the Fourier transform. That is, if $g = f \cdot e$, then its Fourier transform is $G = F * E$ (where $*$ indicates convolution). You could apply this relationship repeatedly to the sinusoid raised to a power to derive the same result as above.


Why is second-order Gaussian called Laplacian Gaussian?

Applying a Gaussian filter is a linear operation, also taking any type of derivative (finite difference) like Laplacian are also linear operations. Now considering the order of applying these operations doesn't matter (a property of convolution), you could combine them in many different ways, e.g. take Laplacian at first then take Gaussian or Gaussian first and Laplacian next, also you could apply Gausian on Laplacian or Laplacian on Gaussian and find a LoG Kernel, then apply the obtained Kernel on image once, which is more efficient in number of computations.


Confusion in proof of Inverse Laplace Transform

To summarize the discussion:

The usual substitution is $s = \sigma + j \omega$ where $\sigma$ is the real part of the $s$ variable and $\omega$ is the imaginary part.
The equation in the image is for the Fourier transform, not the Laplace transform. The Fourier transform can be thought of as the Laplace transform evaluated on the imaginary axis ($\sigma = 0$).
The differential $ds$, when looking at real and imaginary parts distinctly, becomes $d\sigma + j d\omega$.
Any differential is an infinitesimal (very small) change in that variable. $dx$ is a small change in $x$.


The downsampling by a factor of 4 can introduce aliasing (because you upsampled only by a factor 2 in the first stage). If this is the case, it cannot be undone by the following upsampling operation. So in general this system will not produce at its output the original input signal $x[n]$.