Efficiently calculating autocorrelation using FFTs

pichenettes is right, of course. The FFT implements a circular convolution while the xcorr() is based on a linear convolution. In addition you need to square the absolute value in the frequency domain as well. Here is a code snippet that handles all the zero padding, shifting & truncating.

%% Cross correlation through a FFT
n = 1024;
x = randn(n,1);
% cross correlation reference
xref = xcorr(x,x);

%FFT method based on zero padding
fx = fft([x; zeros(n,1)]); % zero pad and FFT
x2 = ifft(fx.*conj(fx)); % abs()^2 and IFFT
% circulate to get the peak in the middle and drop one
% excess zero to get to 2*n-1 samples
x2 = [x2(n+2:end); x2(1:n)];
% calculate the error
d = x2-xref; % difference, this is actually zero
fprintf('Max error = %6.2f\n',max(abs(d)));


Detecting outliers/noise from sensor data

Regarding your next question, I will try to give you some long description explainig a general way of transforming filters defined by their transfer function into discrete domain. If you are not interested in the details, you can skip to the bottom, where I wrote a C code (tested).

The transfer function of a first-order low pass filter (with unit gain) is defined as follows:
$$K(s) = \frac{1}{1+sT}$$
And its angular cutoff frequency defined as $$\omega_c = \frac{1}{T}$$
Therefore the cutoff frequency expressed in Hertz becomes
$$f_c=\frac{1}{2 \pi T}$$
$$T=\frac{1}{2 \pi f_c}$$
Whatever goes above this frequency, will be cut off with gain equal to -3dB/decade for the first order low pass filter.

Either using bilinear transform or by expressing the gain in form of differential equations and then discretizing them with Euler scheme you can obtain a discrete domain equation that expresses the gain. I prefer bilinear transform as it is more straightforward.
In order to apply bilinear transform we perform the following substitution:
$$ s = \frac{1}{T_s} \frac{1-z^{-1}}{1+z^{-1}}$$
Where Ts (not to be mistaken with T) is sampling period, so in your case it is the number of seconds between each sample. I.e. for 44100Hz sampling frequency, Ts becomes 1/44100 seconds.

Now, after performing this substitution and simplifying terms, the discrete time transfer function (Z transform) of your filter becomes:
$$\frac{1+z^{-1}}{(1+2T_s/T) + (1-2T_s/T)z^{-1}}$$
Now lets recall the definition of Z transform. In general form it is expressed as the following:
$$ H(z) = \frac{b_0 + b_1z^{-1}+ b_2z^{-2} + ... + b_Mz^{-M}}{1 + a_1z^{-1}+ a_2z^{-2} + ... + a_Nz^{-N}} = \frac{Y(z)}{X(z)}$$
Where bare numerator coefficients from 1 to M (highest power) and a are denominator coefficients.

Now, a difference equation can be derived in the following way:
$$y(n) = b_0x(n) + b_1x(n-1) + b_2x(n-2) + ... + b_Mx(n-M)  - a_1y(n-1) - a_2y(n-2) - ... - a_Ny(n-N) $$
Where n means n-th sample, y is your output and x is your input. After applying the above, the simplified form of difference equation for your filter becomes:
$$y(n) = (1-\alpha)y(n-1) + \alpha x(n)$$
Where
$$ \alpha = \frac{T_s}{T + T_s}$$

This means that n-th sample (the current one) in filtered signal, called y(n) here, equals one minus alpha times the previous sample in filtered signal plus alpha times n-th sample in your source (noised) signal. In C language:

double lowpass (double x, double y_prev, double alpha)
{
    return (1.0-alpha)*y_prev + alpha*x;
}

int main()
{
    double x[10] = {1,2,3,4,5,6,7,8,9,10}; //signal with noise
    double y[10]; //filtered signal (will be generated)
    double fc = 500.0; //500Hz cutoff frequency
    double Ts = 1.0/44100.0; //for example
    double T  = 1.0/(2.0*3.1415*fc);
    double alpha = Ts/(T+Ts);

    y[0] = x[0]; //initialize first y sample as equal to first x sample

    printf("Input\tOutput\n");
    int i;
    for(i=1; i<10; ++i)
    {
        y[i] = lowpass(x[i], y[i-1], alpha);
        printf("%.4f\t%.4f\n", x[i], y[i]);
    }
}

Without complicating this code, for higher order filtering, you can swap y with x once it is filtered and repeat what's in the loop. It's slow and naive, but otherwise, you'd have to take a power of the transfer function (i.e. K(s)^2 for second order filter) and derive the discrete equation again. Thankfully with some practice, you will see a pattern in this and (even better!) you can find ready derivations of higher order low pass filters in the Internet :)

Sources:
https://en.wikipedia.org/wiki/Low-pass_filter
https://en.wikipedia.org/wiki/Bilinear_transform
https://ccrma.stanford.edu/~jos/fp/Z_Transform_Difference_Equations.html


Is subtraction a common method for efficient filtering?

You are right about the importance of the phase.

Subtraction in the time domain is equivalent to subtraction in the frequency domain. If the frequency domain phase at some frequency is not equal between the complex number being subtracted and the complex number subtracting, even if their magnitudes are the same, the result will be non-zero:

Figure 1. Subtraction of complex numbers of equal magnitude does not always result in a zero. The calculation is black = red - blue.

If the black dot represents the combined frequency response for a system output = input - filtered input at a frequency that you intend to be part of a stop band, you must ensure that the filter is zero-phase. Or if the filter is linear-phase and has in its impulse response a center of symmetry at some time, you must delay the copy of the input by the same amount of time before subtraction. This will not work with efficient infinite impulse response (IIR) filters because of their non-linear phase.


How to calculate the steady state response $y_{ss}(t)$ of a LTI system given the Laplace transform $Y(s)$?

For these calculations, it is better to give the Wolfram Alpha answers:

inverse Laplace transform 1/(s*(s^2 + 4*s + 5))

Which gives the correct expression, consistent with the Final Value expression of 1/5:

$$\frac{1}{5} - \frac{1}{5} e^{-2 t} (2 sin(t) + cos(t))$$

This is very different from the complex (phasor) expression, which could inexplicably have a 1/10, but when evaluated it is still 1/5:

$$\frac{1}{10} i (e^{(-2 - i) t} ((2 + i) e^{2 i t} + (-2 + i)) - 2 i)$$


What is the antialiasing algorithm employed by video cards?

There are various algorithms to do anti-aliasing. I think that the one that you mean is called super-sampling. The idea is that the scene is rendered in x8 resolution of what you eventually see. Then it is low-pass filtered, to eliminate aliasing of higher frequencies in the next step, which is down-sampling x8.

This is described in detail here.

There is a more general discussion of AA algorithms here, and another type of AA algorithm called Multisample AA here.


Steady State of a Difference Equation using MATLAB

Given the trig identities you can rewrite your input as
$$x(t) = 3\cdot \sin(0.25\cdot \pi \cdot t+ \pi/3) - \sin(0.75\cdot \pi \cdot t+ \pi) $$

So it's just two sine wave going in, so it will be two sine waves coming out. You could potentially calculate the value of the transfer function $\omega = 0.25 \pi$ and $\omega = 0.75 \pi$ and manually calculate the output this way.


Is resizing of frequency spectrum a valid method of resampling?

Interpolation by zero padding and computing the IDFT is used in practice but is not optimum and will have greater error over other approaches with the same or less resources, for the same reason that filtering by zeroing frequency bins in a DFT is not recommended. This point is detailed at this post Why is it a bad idea to filter by zeroing out FFT bins?.
Methods that perform better in terms of efficiency and performance include polyphase resampling structures using filter coefficients based on least squares algorithm, and specifically for multiband fiters that concentrate rejection at the image locations associated with the spectral distortion.
For further details on such approaches see:
How to implement Polyphase filter?
and
Computational Complexity of Polyphase Resampling
When a Sinx/x droop can be tolerated (or easily compensated for, see how to make CIC compensation filter) cascade-integrator-comb (CIC) interpolators are highly attractive for their simplicity and efficiency (there are many posts here that detail the beautiful CIC (otherwise called the Hogenhauer filter after the original author) further, see:
CIC Cascaded Integrator-Comb spectrum
Puzzled by CIC filters
Interpolating/Decimating CIC Filter Group Delay


PyWavelets CWT: normalization? Vs Scipy?

Let's follow the math from incubation to delivery. It begins with psi, a rescaled morlet2 (as shown previously) at a scale $a=64$, and $\sigma=5$:
$$
\psi = \psi_{\sigma}(t/a) = e^{j\sigma (t/a)} e^{-(t/a)^2/2} \tag{2}
$$
gets integrated and L1-normalized: -- (see caveat2 below)
$$
\psi = \psi_{\text{int}}(t) = \frac{1}{a} \int \psi_{\sigma}(t/a)\ dt \tag{3}
$$
The norm is off by a scalar factor, but what matters is it's same for different $a$; you can verify it's so (for large $a$, increase integration limits). Note that it's the wavelet that's L1-normalized, not the integrated wavelet $\psi_\text{int}(t)$; the distinction is essential, as we'll see.
$\psi$ is then  convolved with the signal $s$, and finite-differenced and negated; from here, this is:
$$
 -\text{DIFF}(\text{conv}) = [s \cdot (\psi_{\text{int}}(t - 0) - \psi_{\text{int}}(t - 1)),\ ...] \tag{4}
$$
Inspecting closely, this is equivalent to undoing the former integration, with certain numeric differences (addressed later), so we're left with
$$
 \psi = \frac{1}{a} \psi_{\sigma}(t/a) \tag{5}
$$
Lastly, * sqrt(scale):
$$
 \boxed{\psi = \frac{1}{\sqrt{a}} \psi_{\sigma}(t/a)} \tag{6}
$$
which is thus ultimately L2-normalized. We can verify this on a test signal; L2 is biased in amplitude/power representation toward lower frequencies, suppressing high frequencies. If above is right, not only can we see this, but apply a fix by removing * sqrt(scale) from code:

Implementation-specific caveats discussed next.

How reliable is the normalization across scales? -- Not great (-- code2):

On right-hand side we see the best case using limited number of samples, on left we see the resampling implementation. Clearly due to inherently poor resolution of small number of samples, our norms don't do well for low scales - but the recomputed variant at least maintains a mean norm approximately constant across them.
More interestingly, the resampled one suddenly appears quadratic after $\log_2(a)=6$; what's the deal? Recall "It begins ... at scale $a=64$"; this is when we resample with more samples than in original (int_psi):

whose diff is

In 1LR we see int_psi at scale=256 is resampling the same values from int_psi multiple times, which manifests as zeros in diff (2LR).
Next we square the values to get an idea of the L2 norm, and compare against the L2-normalized values in 3LR; we find normed values blowing unnormed values out of water. Is this excessive, or appropriate?
Compare against recomputed at same scale; in 4LR, we at least see recomputed's shape - it's no longer negligible. This suggests the resampling to be overestimating wavelet amplitude, and since in norm computation this is significantly offset by half the values being zero due to diff, then the overestimation must be great.

How great? We can compare peak wavelet amplitudes across scales of the unnormed resampled diff(int_psi_256) and recomputed w (morlet2/scale); norming changes both by the same factor (*sqrt(scale)), won't make a difference.

Surely enough, again at $a=64$ we see critical behavior; peak amplitude of resampled saturates at a minimum, whereas the recomputed continues to decrease. The difference at $a=256$ is more than great enough (~12db=4x, if squared -> 16x) to offset the zeros in diff, hence why the L2 norm grows.
But notice the L1 norm remains well-behaved at higher scales. To find out why, we must precisely compare the proportion of zeros with the overestimation factor. We already know latter to be 4x; for former, I took a non-flat region (otherwise it's zeros anyway) and found ratio of zeros to number of samples: 0.75. Which means, 1 in 4 values are non-zero, and each are 4x greater ... the exact proportion needed to retain the L1 norm.
That zeros proportion isn't very clear from above plots - so instead, dots:

If diff undoes cumsum, why do either at all?
Good question. In fact, by ridding of both, I found the quality degradation (zipper) effect disappear:

But something else also disappears - the high frequency scales (left), and the sign flips (I undid this by removing minus at coef = -). I partly offset this by replacing * sqrt(scale) by 1 / sqrt(scale) (right), but that gets us to the L2 normed case (1 / scale is overkill), and is a baseless modification in light of everything found so far. So can it be done? Sure, one can implement without the int-diff dance, but I won't be exploring it here.

Caveat: effective norm
Recall, conv gets trimmed to len(x), so the entire wavelet may not participate in computation, whereas above computations assume otherwise. The assumption is truer the longer x is - but here's what happens if len(x) == 512, with same wavelet precision=10 (2^10=1024 int_psi), compared against recomputed wavelet:

A mess! Both distort significantly, possibly implying it's inherent to great enough scales pinned against short signals. A simple workaround is zero-padding the signal, but I wouldn't vouch for it without understanding the caveats (which I don't for wavelets, but do for DFT).
Furthermore, note that the norm isn't only per-scale, but also per-shift ($b$); less of the wavelet is in the signal when it's just entering as opposed to when its center is at signal's center.
But keeping it per-scale, I confess I cheated in the earlier example with $f=4$; initially I picked $f=1$, but the L1 norm plot didn't look nice, and I didn't understand why. Since now I do, here are both L1 and L2 per scale for $f=1,64$:

Compare against earlier plot; note how the left modes (low scales) of both norms are unaffected, only the high scale (low frequency) ones are now much greater (L2) and lower (L1), just as expected.

Caveat 2: * step doesn't do the L1-normalizing
Instead, it results from diffing the resampled integrated wavelet; the * step only sets the correct initial scale, and can be omitted without losing the norm. This does not happen with recomputing the wavelet, or resampling the unintegrated wavelet.
(Note that recomputed scipy.signal.morlet2 is L2-normalized, but only because it's explicitly normalized; the diff-resample-int is "self-normalized").


You are on the right track but need perhaps to quantify your goal in mathematical terms.

The sorted coefficients is incidently also the inverse of the cumulative function (see for instance the eq 6 of this paper).

Then to quantify the sparsity of your signal you could

compute some measure using the cumulative distribution function, such as the ratio of coefficients above a threshold (the 5% max quantile),
and compare with a control condition, like doing the same process
with a random signal


Why a DFT of two sinusoids is very noisy even with frequency sampling 5 times higher?

As far as I understood your question, you are complaining about the windowing effect and its consequence of smearing of the spectrum.

So you have a mathematically defined pure sine wave that extends from minus infinity to plus infinity in time, whose theoretical spectrum, you expect to be, a single line in frequency. But you eventually represent that infinitely long ideal sine wave with a finite length record (the windowed version, aka truncated) in practice and take the DFT of the finite length version instead. Then the result will be a broadened single line into a sinc waveform mostly having nonzero amplitudes at all frequencies in general which is what you called noise?

That's the most fundamental manifestation of the effect of time domain windowing (aka aperture effect) on the frequency spectrum; the smearing of the single frequency line and also leakage when multiple tones (including DC) co-exist.

So those are not noise, but an irrecoverable side effect of the practical utilization of finite length DFT. Note however that you might be adding your computational error on top of that. Hence for best verification, you better test your algorithm with an established library from Matlab / Octave kind of platforms.

The following plots are generated by Matlab from your uploaded data, to test your problem and it confirms with the explanation; i.e., those nonzero values you see in the spectrum are due to the windowing effect...

You cannot get rid of it but the effect will minimize by using longer data, however is not the most practical approach. Instead by using special window types such as Hamming, Hanning, Blackman, Kaiser etc, you can tradeoff certain features of it. Incidentally you are using a rectangular window, whose peaks are sharpest but those tails (side lobes) are also largest. All other windows will supress the tails (side lobes) but widen the peaks...

And this is the result of Hamming window applied to your data :

Note that your signal $x[n]$ seems 1 sample longer than it would be...


A very basic question about processing high frequency signals

In digital signal processing, one almost never deals directly with high-frequency signals. The reason is that the frequency band that a signal occupies is completely irrelevant to the information that the signal conveys; the information exists in the envelope of the signal, which usually can be handled either by a generic CPU, a DSP, FPGA, or ASIC.

So, let's say that your GPS signal is $s(t)=x(t)\cos(2\pi f_c t)$, where $f_c\approx1.5\,\text{GHz}$. The actual information is in the envelope $x(t)$; the purpose of the cosine is simply to shift ("upconvert") $x(t)$ to some appropriate band. In the receiver, the signal $s(t)$ is "downconverted" to some lower (or even zero) frequency $f_I$, but this is done in the analog domain (i.e. by RF electronics), to obtain $s_1(t)=x(t)\cos(2\pi f_I t)$. The signal $s_1(t)$ has a small bandwidth, so it requires a small sampling rate, that can be handled in the digital domain.

Note that for quadrature signals the math is very slightly more complicated; see this question and the complex envlope article in wikipedia.

Note also that a similar process exists in the transmitter: the upconversion to RF is done in the analog domain, after the discrete envelope has been converted to an anlog signal.


MATLAB optimization for a subroutine function

To be clear, you have a Matlab "black box" function which takes a $10\times 1$ input and generates a function value and you want to find the minimizer. Such optimization problems are hard to solve in the sense that we cannot apply any heuristics or tricks because we don't know anything about the function you are trying to optimize. Moreover, you don't have any gradient information either.

I would start by trying a derivative free method like Nelder-Mead (fminsearch in Matlab). Do you have a box constraint on the 10 variables, i.e. do you know lower and upper limits for each coordinate? If so you can do a grid search. Brute force search in 10 dimensions may be hopeless and take a very long time. Instead you can use coordinate descent i.e. change one coordinate at a time to the next grid value and step in a direction that reduces the function value.

You can also try other heuristic global optimization routines such as genetic algorithm, simulated annealing and particle swarm optimization. If you have access to Matlab's global optimization toolbox, you can try their "Optimization app" graphical interface to quickly test some of these out.


Loudness detection for natural conversation

I would look into the EBU R-128 method of loudness measuring. Based on an ITU algorithm, it is used to determine loudness for TV which consist mostly of dialogue.

ITU Docs

EBU R-128 Docs


I think you need to change your code as follows

n = 3;
N = 2^n;
signal = rand(1,N);

fft_sig = fft(signal);
fft_sig1 = fft_sig(1:2:end);  % even/odd frequencies
fft_sig2 = fft_sig(2:2:end);

sig1 = ifft(fft_sig1);        % the actual time-domain signals
sig2 = ifft(fft_sig2);

stem(0:N-1,  abs(fft_sig),'k');
hold on
stem(0:2:N-1, abs(fft_sig1), 'b*');
stem(1:2:N-1, abs(fft_sig2),'rs');


Doubts and some confusion on variance for complex rv

If the imaginary and real components each has a variance of 0.5, then what would be the total variance?

The variance of the complex RV in that case would be $0.5 + 0.5 = 1$.

  If on the other hand, if $w∼CN(0,2σ^2_w)$ then does it mean that the real and imaginary components each has variance 1, so the total variance is 2. Is my understanding correct?

Each has variance $2σ^2_w/2 = σ^2_w$.

In other words, when you add two RV, the variance of the result is the sum of the variances of each RV.

  Is there a rule of thumb whether the variance should be 0.5 for each component?

When talking about Rayleigh fading coefficients, you want the variance per dimension to be 0.5 because that results in the average transmitted power being equal to the average received power (that is, the channel neither creates nor absorbs energy).


You haven't specified what part exactly of the solution you didn't understand.

The purpose of the posed question is to show the spectral effects of signal downsampling, in this case by a factor of 2.

The original discrete signal is $x[n]$ (whose DTFT is $X(e^{j\omega})$), with a length of N samples, and the downsampled signal (by a factor of 2) is $g_7[n]=x[2n]$, so $g_7[n]$ is made from the even samples of $x[n]$, with a total length of N/2 samples.

The first experssion of the solution is simply the DFT (by definition) of $g_7[n]$. The solution then proceeds to some series manipulation, then uses the DTFT-derived definition of DFT to show the final expression.

The final expression of the solution shows that downsampling the signal by a factor of 2 produces a double-frequency image of the original signal in the frequency domain. This alludes to the importance of anti-aliasing filtering prior to signal downsampling.


Choosing points in the horizontal plane with proper spacing seems like the right approach. Aneochic HRTFs are probably the easiest but typically lead to "inside the head" localization. Depending on your application this may or may not be a good thing. If you want externalization it's best to use binaural room impulse responses measured in a space that's similar to where the actual playback happens (there are psychological reasons for that). However these binuaral room impulse responses are a lot more complicated and require significantly more computational effort.

You cannot binauralize a measured monaural impulse response after the fact. A real binaural room impulse response consists of the direct sound plus many reflections that all come from different directions. So each part of the room impulse response contains many different anechoic HRTFs. There are however computer programs that allow you to put in a geometry, acoustic properties and then they will calculate a binaural room impulse response.


When you downsample a signal you must make sure that the Nyquist condition remains satisfied.

For example, say that signal $x$ has frequencies between $0-25$ Hz, and that the sampling frequency of $x$ is $100$ Hz. The Nyquist frequency is $100/2 = 50$ Hz, thus you could downsample $x$ by a factor of two with no problems.

Now, say that you would like to downsample $x$ by a factor of $4$ instead of $2$. This would mean that the new sampling frequency would be $25$ Hz. Since the new sampling frequency is $25$ Hz, it must not contain frequencies above $25/2=12.5$ Hz. However, $x$ had frequencies between $0-25$ so there would be aliasing. In order to avoid that, you must low-pass filter $x$ first between $0-12.5$ Hz.

The process of low-pass filtering + downsampling is called decimation. You can find more information about it on Wikipedia.

Upsampling would be basicly interpolation, so I would avoid that where possible. I did not look at the URLs you provided, please phrase your problem in full here.


Attenuating peaks of a signal in a specific direction

My immediate approach to this: try something that I'd call a running-minimum filter:

Look at a fixed amount of $N$ samples around the current sample. If the current sample is significantly below that, modify the sample to be closer to a low-pass variant of the original signal.


Fourier transform of ${\Pi}_{a}(t)\cos(2{\pi}f_{0}t)$

Your last line is wrong. Note that you have

$$H(f)\star\delta(f-f_0)=H(f-f_0)\tag{1}$$

for any $H(f)$.

So your last line should be

$$\begin{align*}
&a\,\text{sinc}(2af)\star\delta(f-f_{0})+a\text{ sinc}(2af)\star\delta(f+f_{0})=\\&a\,\text{sinc}(2a(f-f_0))+a\,\text{sinc}(2a(f+f_0))\tag{2}
\end{align*}$$

You basically get copies of your original spectrum centered at $f_0$ and $-f_0$.


You could try calculating Shannon entropy of the spectrum. Normalize the Fourier transform $f(x)$ of your signal so that $\int_{-\infty}^\infty |f(x)|^2\, dx = 1$ and calculate Shannon entropy as $- \int_{-\infty}^\infty |f(x)|^2 \log |f(x)|^2\, dx$. You can clamp $|f(x)|^2 \log |f(x)|^2$ to zero for very small values of $f(x)$ if its logarithm blows up and ruins the calculations. Zero-padding before windowing & FFT will give you more frequency resolution and a better approximation, especially if you have narrow peaks in the spectrum.

If you want something more standard, have a look at spectral flatness. But it seems even less stable numerically.


Sinusoid amplitude loss in power spectrum estimation using Pwelch method

I think your difficulty lies in having an overlap size of 0.

I'd recommend having an overlap that gives a Constant Overlap Add result when the windows are laid on top of each other.

For a rectangular window that is 0.
Hann it is 50% blocksize
Hamming seems to have many, but try 75%

Overlap values for different Windows

And don't forget that the PSD 'bin' represents all the content within that delta-f.  For example, if you had a signal at 100 Hz and 100.25 Hz, the 1 Hz delta-f wouldn't resolve them to 2 spikes .

RMS from a bin is sqrt(Area of bin).


How to transform SoundWave Picture to sound?

I don't think "I love you" is doable with a sound wave image, because that's already very roughly about 200 cycles of the fundamental frequency of the voice, and you need also the harmonics so that it sounds right and can be understood. Too much horizontal resolution is needed to pack that into a tattoo.

You could probably pull it off as a spectrogram.

I've seen ads for a commercial service that presumably stores at a server a database of sound clips each associated with a sound wave image. Then their phone app sends the tattoo photo to the server and plays back the sound clip associated with the image that most resembles the tattoo photo. Technically that's not far different from having the sound clip number tattooed as a bar code.


Is mean squared error or absolute error on two spectrograms of two different wav files enough to compare how similar they are?

My suggestion would be to partition the spectrum according to the well tempered scale, one frequency band for each half tone. You generate a pattern of 1s and 0s for each chord you want to be able to test, with a maximum of six 1s.
When a chord is played, just do a FFT of the whole thing, normalize the energy and check each band for a threshold.
Your "fitness" value then could be directly calculated from the correct strings vs. wrong strings in some way.


How to add complex WGN to complex damped exponential and compute SNR?

One suggestion is to leverage the vector processing in Matlab and to not add the noise in a for loop, but given the vector of your signal of length $N$ you can simply create the white Gaussian Noise using:
noise = std * sqrt(2)/2 * randn(N,1) + std * sqrt(2)/2* 1j*randn(N,1)
(or randn(1,N) if you prefer row vectors instead of column vectors)
Where std represents the total standard deviation of the noise.
Then the noise can be added to your signal simply using
y = x + noise
Next be aware that the total power as a white noise and the square of the standard deviation std in whatever units are chosen is distributed evenly over the entire frequency range given by -$f_s/2$ to +$f_s/2$ where $f_s$ is the sampling rate. SNR of interest is typically that portion of the noise that shares the same spectrum as the signal. (So if the signal only occupied 1/4 of the total spectrum as sampled, then the noise would be 1/4 or 6 dB lower than that predicted by the total variance of noise).


Fourier Series of Aperiodic convolution of periodic functions

okay, now that the popcorn is eaten and my fingers are clean...

$$ x(t) = \sum\limits_{m=-\infty}^{+\infty} c_m^x e^{j 2 \pi (m/T) t} $$

$$ y(t) = \sum\limits_{m=-\infty}^{+\infty} c_m^y e^{j 2 \pi (m/T) t} $$

then

$$\begin{align}
 z(t) &= \sum\limits_{m=-\infty}^{+\infty} c_m^z e^{j 2 \pi (m/T) t} \\
\\
  &= x(t) \circledast y(t) \\
\\
  &= \int\limits_{-\infty}^{+\infty} x(\tau) y(t-\tau) \, d\tau \\
\\
  &= \int\limits_{-\infty}^{+\infty} \sum\limits_{m=-\infty}^{+\infty} c_m^x e^{j 2 \pi (m/T) \tau} \sum\limits_{n=-\infty}^{+\infty} c_n^y e^{j 2 \pi (n/T) (t-\tau)}  \, d\tau \\
\\
  &= \sum\limits_{n=-\infty}^{+\infty} \sum\limits_{m=-\infty}^{+\infty} c_m^x c_n^y  \int\limits_{-\infty}^{+\infty} e^{j 2 \pi (m/T) \tau} e^{j 2 \pi (n/T) (t-\tau)}  \, d\tau \\
\\
  &= \sum\limits_{n=-\infty}^{+\infty} \sum\limits_{m=-\infty}^{+\infty} c_m^x c_n^y \int\limits_{-\infty}^{+\infty} e^{j 2 \pi ((m-n)/T) \tau} \, d\tau \ e^{j 2 \pi (n/T) t} \\
\end{align} $$

this integral: $\int\limits_{-\infty}^{+\infty} e^{j 2 \pi ((m-n)/T) \tau} \, d\tau$ ain't converging when $n \ne m$ and certainly not when $n=m$.

you will not be getting finite values for $c_m^z$.


How to decompose a signal into exponentally decaying sinusoids?

For a continuous-time signal this decomposition is precisely the Laplace Transform, and for a discrete time signal, the z-transform. Both decompose a time domain waveform into decaying and growing exponentials.
This is clear from observing the formulas and seeing the underlying operation similar to generalized correlation which is a sum of products:
Here is the general form for the Unilateral Laplace Transform:
$$\mathscr{L\{x(t)\}} = X(s) = \int_0^\infty x(t)e^{-st}dt$$
Where $s= \sigma + j\omega$
Note that $e^{-st}$ represents a rotating phasor that is growing or decaying with time!
If it isn't yet clear, $e^{-st} = e^{-(\sigma + j \omega) t}= e^{-\sigma t}e^{-j\omega t}$, and thus we have a rotating phasor ($e^{-j\omega t}$) that is growing or decaying $e^{-\sigma t}$ both as a function of time. If that is confusing, it may further help to know that the general form of the complex exponential $Ke^{j\theta}$ describes a vector (phasor) of magnitude K and angle $\theta$. Rotating phasors are unique frequency tones and fundamental elements (basis functions) of the Fourier Transform. Such a phasor rotating counter-clockwise on a complex plane is a positive frequency  (the phase of the phasor is increasing linearly with time, and frequency is $d\theta/dt$), and one rotating clockwise would be a negative frequency. If we observed the real axis of any such rotating phasor, we would observe a sinusoidal waveform which is how we are typically first introduced to "frequencies". Real sinusoids consist of two such tones, a positive and negative frequency, as given by Euler's equation: $2cos(\omega t) = e^{j\omega t} + e^{-j\omega t}$. Thus complex tones of the form $e^{j\omega t}$ are even simpler subcomponents of real tones and once that is well understood, thinking of "frequencies" in terms of complex tones rather than real sinusoids greatly simplifies signal processing and leads to much greater possibilities! Complex need not be "so complex"!!

Let's call $e^{st}$ our hypothesis of what may be a component of $x(t)$ for any given value of $s$. We vary $s$ over all possible values to establish the result $X(s)$. In the Unilateral Laplace Transform, as given by the formula above, we multiply our signal $x(t)$ in time with $e^{-st}$ and integrate for all time, $t\ge 0$. Note the product is negated in the exponent: We spin exactly opposite to our hypothesis, and decay if it is growing and vice versa at the exact same rate as the mechanism to locate all these fundamental components! When a component is found, the rotation is stopped and the growth or decay is leveled out!

For any particular value (point)  $s_1$ of $e^{s_1 t}$ that may be a fundamental component of our signal $x(t)$, meaning $x(t)$ contains as part of it one such fundamental decaying, rotating complex exponential; then for that component, it's product with $e^{-s_1 t}$ will be $e^{s_1 t}e^{-s_1 t} = 1$ for all time (for t>0)!... which then when integrated to infinity as given by the formula for the Unilateral Laplace Transform would go to infinity-- but notice any other value for $s$ where $\sigma > \sigma _1$  (for $s_1 = \sigma _1 + j\omega_1$) would necessarily decay and therefore converge to a finite value. Thus we see how we get distinct points that go to infinity (poles) and a vertical line as the boundary with everything to the right in the "Region of Convergence".  Anything to the left of this vertical line (as given by the left most poles in the solution to the Laplace Transform) will not converge and therefore is an invalid region where no solution exists (and therefore our signal has no components there).
If plotted the Laplace Transform would be a surface of complex results for every combination on the 2 dimensional surface given by $(\sigma, j\omega)$. We often just depict this as poles (points going to infinity) and zeros (points equal to zero) since every other point on that surface is defined and unique by those pole and zero locations (but we may then miss understanding the full picture).
If we restrict $s$ to be only the $j\omega$ axis, then the Laplace Transform reduces to the Fourier Transform, which is a correlation with every possible rotating phasor that doesn't decay or grow (stays constant for all time). Specifically for the Unilateral Laplace Transform where $t\ge 0$, this would also involve multiplying the signal by a unit step ($u(t)$):
$$\mathscr{F\{u(t)x(t)\}} = \int_0^\infty x(t)e^{-j\omega t}dt$$
Notice how similar that equation is to the Unilateral Laplace Transform first given, and observe the surface plot of the magnitude of an example Laplace Transform below (the results of the Laplace Transform are themselves complex values, so each result on the surface has a magnitude and phase) for a low pass filter, and if we slice that 3 dimensional shape along the $j\omega$ axis we end up with the frequency response (Fourier Transform of the filter's impulse response, where $X(s)$ is the Laplace transform of the filter's impulse response), as given in the graphic below in the upper right hand corner (a low pass filter), and we see from the low pass filter that the impulse response consists of a decaying sinusoid (if the poles were on the $j\omega$ axis, this would be a constant level sinusoid, but given they are shifted to the left, that sinusoid is decaying with time.)

For discrete-time systems, the z-transform does the exact same thing, just with discrete phasors rotating in discrete time that are growing and decaying and instead of an s-plane we have a z-plane where the unit circle corresponds to the $j\omega$ axis. With discrete processing we use algorithms that compute the Discrete Fourier Transform (DFT) such as the ubiquitous FFT. The corollary for computing the Z-Transform is the chirp Z-transform (CZT), and for efficient implementations of this see Bluestein's algorithm.


What measurement errors occur in an inteferogram?

You did an excellent job in presenting the background and context of your question and your listing of error sources is fine. But FTIR also presents some other issues, particularly given the relatively poor SNR performance of IR detectors compared to UV-Vis regime detectors. So two things. First, I recommend this book: P.R. Griffiths, J.A. de Haseth, Fourier Transform Infrared Spectrometry, 2nd Ed., Wiley-Interscience, John Wiley & Sons, Inc., NY, ©2007.

Second, when I simulate FTIR using my optical calculus software, I generally include light source temperature fluctuations, jitter on the moving mirror and, of course, detector and preamplifier noises. The next few figures show the simulation model (with the 3 noise souces I stated), which looks somewhat like a block diagram, but is actually the simulation program, and some results for ethyl acrylate in the 400 to 4000 wavenumber range. Figure 1 follows:

Above is Fig. 1, the FTIR simulation program with ethyl acrylate as the test substance. The simulation conditions are shown in Fig. 1.

Figure 2 is the last of the ten 16k interferograms generated. It takes about 2.5 minutes per interferogram on my 2006 iMac 24. (It would be much faster on a modern Windows-based PC, but speed isn't everything and I am currently holding off buying a new PC.) Figure 2 follows:

The next figure (Fig. 3) is an expanded view of the center of the interferogram:

Now Fig. 4 shows the pseudo-transmittance of the ethyl acrylate test substance:

Then Fig. 5 shows the pseudo-transmittance of the background:

Finally, dividing the two pseudo-transmittances and taking the negative log yields the absorbance spectrum of ethyl acrylate, as shown in Fig. 6:

My software runs on PCs (and my real old Macs) and is free, with full commented source code. However, it is an add-on to a commercial simulation program called ExtendSim® (with which I am not financially associated: I am just a very early user).

If I can find any relevant and detailed papers discussing FTIR considerations in particular, I will edit this answer to include them. And it would benefit me as well, even though the spectroscopy I used to do was not in the IR.


Is convolution distributive over multiplication?

Addition is (by nature) additive. Binary operations like multiply and convolution are both somehow of multiplicative structure. Thus, I suspect there is no generic distributivity (apart for Boole or Boolean algebras, where the restriction on allowed values shrinks the problem, see below).

However, there is a sort of commutativity of multiply and convolution, as follows. It is interesting to first look at such operator  properties from an algebraic point of view. Properties are for instance, for a generic operator binary $\bigcirc$:

commutativity: $a \bigcirc b = b \bigcirc a $
associativity: $a \bigcirc( b \bigcirc c) =  (a \bigcirc b) \bigcirc c  $
alternativity: $a \bigcirc( a \bigcirc b) =  (a \bigcirc a) \bigcirc b  $ or $(a \bigcirc  b) \bigcirc b =   a \bigcirc (b \bigcirc b)$
flexibility: $(a \bigcirc  b) \bigcirc a =   a \bigcirc (b \bigcirc a)$

Distributivity of a binary   operator $\diamond$  over operator $\bigcirc$  is more involved: it can be left-distributed: $$a \diamond (b \bigcirc c) = (a \diamond b) \bigcirc (a \diamond c)$$ or right-distributed: $$(a \bigcirc b) \diamond c = (a \diamond c)  \bigcirc (b\diamond c)$$

Division is right-distributive over add, not left-distributive. Distributivity of multiply over add is at play in rings and fields. But the converse does not apply. There are mathematical structures where
two binary operators distribute over each other, like Boolean or switching algebras. Apart from that, I cannot remember of standard cases of multiply-like $\ast$  or $\cdot$ that would distribute. However, I will details some specific cases, useful to signal processing, ending with the Hilbert transform and Bedrosian theorem.

In the case of multiplication and convolution, or $a$ scalar, this is associative, as :

$$a.[b(t)\cdot c(t)] = [a.b(t)]\cdot c(t)= b(t)\cdot [a.c(t)]$$
and
$$a.[b(t)\ast c(t)] = [a.b(t)]\ast c(t)= b(t)\ast [a.c(t)]$$

You see that you cannot directly"distribute it", as this would become "bilinear".
Convolution and point-wise product are linear. So, if you consider that $c(t)$ (or  $b(t)$ by symmetry), there exists a linear operator $\Lambda$ such that:

$$\Lambda(c(t)) = a(t)\ast[b(t)\cdot c(t)] $$

Whether $\Lambda$ has a simple closed form expression is a difficult matter. As a side note, where $\mathcal{H} $ denotes the Hilbert transform, the Bedrosian theorem (A product theorem for Hilbert transforms states that (under some conditions):

$$\mathcal{H}(a(t)e^{i\theta t}) = a(t)(\mathcal{H}(e^{i\theta t}))\,.$$

where $a(t)$ is an (low-pass) envelope, and $e^{i\theta t}$ a modulation. More generally, if a low-pass signal $x_\flat$ and a high-pass  $x_\sharp$ have non-overlapping spectra, then under Bedrosian theorem:

$$\mathcal{H}(x_\flat(t)\cdot  x_\sharp(t))=x_\flat(t)\cdot \mathcal{H}(x_\sharp(t))$$

Now, let us remind that the Hilbert transform can be seen as the convolution with the distribution ($\operatorname {p.v.}$ is the Cauchy principal value):

$$h(t) =\operatorname {p.v.} {\frac {1}{\pi t}}$$

Thus:

$$h(t) \ast (x_\flat(t)\cdot x_\sharp(t))  =  x_\flat(t) \cdot(h(t)\ast x_\sharp(t))$$

This is the only practical case I can remember of a switching of operators in signal processing, more a form of limited commutativity.


Synchronize device to audio stream - Approach

One problem with an acoustic link is frequency dependent multi-path fading/reinforcement, which can vary widely as the source and mic positions change, or as the room acoustic vary (moving hard objects, etc.)

So I might try a "click" that has a wide and known pattern of several non-harmonically related spectra, and check the Hamming distance between the transmitted and received spectral pattern peaks.

Use a non-harmonically related set of frequencies so that the harmonics of any clipped spectra will be mostly ignored (depending on window duration), and because multi-path fading can act as a comb filter notch.


Compressive Sensing - Incoherence Property

I am trying to answer your question about incoherence here rather than update my previous answer on another question of yours.

Compressive sensing requires low coherent pairs. So the lower $\mu(\Phi,\Psi)$, the better. Actually is $\Phi$ is spike basis (identity matrix) with $\phi_k(t) = \delta(t-k)$, and $\Psi$ is Fourier basis with $\psi_j(t) = 1/\sqrt n e^{-i \cdot 2\pi \cdot jt/n}$, which is just the example I showed you in another question, $\mu(\Phi,\Psi) = 1$, and the maximal incoherence is achieved. It is similar that if you apply the spike basis and the orthogonal basis obtained from SVD.

Other coherence includes but not limited to: coherence between noiselets and Haar wavelets is $\sqrt2$; coherence between noiselets and Baubechies D4 and D8 are $2.2$ and $2.9$, respectively; random matrices are largely incoherent with any fixed basis $\Psi$.

You may refer to Candes-Romberg theorem for more details on incoherent sampling.


Are there any analytic approaches to the design of quadratic polynomial filters?

Not really an answer (hence this is community wiki), but I thought we should capture @Mohammad's code:

%Mohammad Z

%Teager-Kaiser Non-Linear Energy Operator.
function [out] = TKEO(x)
    N = length(x);
    x = x(:).';
    temp = x(2:N-1).^2 - x(3:N).*x(1:N-2);
    out = [temp(1) temp temp(end)];
end


OFDM systems are usually implemented using the FFT and its inverse. At the transmitter, information symbol values are assigned to the various OFDM subcarriers, resulting in a single OFDM symbol in the frequency domain. To obtain a time domain that can actually be transmitted over the channel, an inverse FFT is applied. The receiver on the other end of the link performs the reverse (it computes a forward FFT) on its observed signal, resulting (hopefully) in the original frequency-domain symbols, thus accomplishing OFDM demodulation.

Why add a prefix?

The short answer: multipath through the communication channel. You can think of multipath as a linear filter applied to the transmitted signal via propagation through the communication channel. Thus, the signal observed at the receiver is the original one from the transmitter, convolved with the channel response (and corrupted by noise). This results in intersymbol interference (ISI), which degrades system performance.

Techniques to undo the effects of this channel response are called equalizers, and are widely used in practical communications systems. However, their implementation and analysis can be somewhat complicated. It would be nice if, for your OFDM system, you could come up with a simpler equalization method.

Another way to mitigate multipath would be to insert some sort of guard time between symbols (a prefix) in order to allow the delay spread of the channel from one symbol to die out before the next. This strategy is used in many OFDM systems; the added prefix allows time between symbols for the channel's impulse response to die out, reducing or eliminating ISI. For most other systems, adding this overhead on each symbol would be prohibitively inefficient. However, since OFDM's effective symbol rate is much lower (since many modulated subcarriers are transmitted in each OFDM symbol time in parallel), the overhead of this approach is proportionally lower as well.

Why use cyclic prefixes?

Recall the fact that the OFDM symbols are modulated onto subcarriers in the frequency domain. Also recall the relationship between convolution and the discrete Fourier transform (DFT): convolution in the time domain corresponds to multiplication in the frequency domain. Therefore, the effects of the convolution with the channel response show up as pointwise multiplication of the frequency-domain subcarriers. This makes for a very simple equalizer structure: we just have to multiply each subcarrier by a (per-carrier) constant value at the receiver in order to undo the effects of multipath. This is typically called a "one-tap" equalizer, and is equivalent to assuming that any fading is flat across each subcarrier's bandwidth.

However, recall that multiplication in the DFT domain corresponds to cyclic convolution in the time domain. So, in order for this strategy to work, we would need some way of turning the channel response into a cyclic convolution, not a linear convolution as typically happens in the real world. That's where the cyclic prefix comes in.

The implementation: for each time-domain symbol generated by the transmitter, a cyclic prefix is inserted before it. That cyclic prefix takes the last $L$ samples of the $N$-sample time-domain symbol and concatenates it to the front of the symbol to be transmitted. Therefore, for each $N$-sample OFDM symbol, the transmitter actually emits $N+L$ samples. The length of the prefix $L$ is chosen to be greater than the expected multipath spread.

Why does it work?

With the cyclic prefix inserted, you can see conceptually how the linear convolution induced by the channel turns effectively into circular convolution; the cyclic prefix, which includes the end of the symbol, is smeared onto the beginning of the symbol by the channel's impulse response. This effect is quite literally defined as circular convolution, so the goal of achieving simple equalization is accomplished.


Can I get more information from this signal other than magnitude?

There's a lot of answers -- look at what's different between the metal and no-metal case, and try to make some signal processing that'll do the best job of differentiating between them.
Something that I suspect that will work well would be to start a clock as soon as the pulse goes negative (possibly according to the transmitter timing, if that will be reliable).  Capture the height of the peak, then integrate the area of the whole response over time.  That peak to area ratio should be very indicative of the presence of metal.
My gut feel is that the ratio of the peak to the integral of the overall response should be more indicative than the actual value of the integral -- but if the pulses are repetitive enough, just that integral may be enough.  If you want to go down that road, be sure to try to detect every different thing you may be concerned with, at every different range.
It may be more reliable to capture the peak by integrating for a short time, and comparing that number to the integral of a longer time span.
Continuing in this direction, if you were to integrate the pulse over time and look at the shape of the resulting integrated pulse it should be very indicative of what's going on.  Once again my gut feel is that you want to look at ratios -- I'd normalize the integrated pulse to its final value.  The points at which it reaches 1/4, 1/2, 3/4, etc., of its final value should be very telling.


retrieve phase of an extrapolated periodic 1D signal

if you know, in advance, the frequency of the sinusoid, so only the phase (and perhaps the amplitude) is unknown, the simplest noise-immune method i know of to get amplitude and phase is (i'm gonna do this in continuous time and let you translate it to discrete time and eventually to MATLAB):

$$ x(t) = \sin(\omega t + \phi) + n(t) $$

compute

$$ \begin{align}
a & = \int\limits_0^{\frac{2 \pi}{\omega}} x(t) \sin(\omega t) dt \\
& = \int\limits_0^{\frac{2 \pi}{\omega}} \left(\sin(\omega t + \phi) + n(t) \right) \sin(\omega t) dt \\
& = \int\limits_0^{\frac{2 \pi}{\omega}} \sin(\omega t + \phi) \sin(\omega t) dt + \int\limits_0^{\frac{2 \pi}{\omega}} n(t) \sin(\omega t) dt \\
& = \frac{1}{2} \int\limits_0^{\frac{2 \pi}{\omega}} \cos(\phi) dt - \frac{1}{2} \int\limits_0^{\frac{2 \pi}{\omega}} \cos(2\omega t + \phi) dt + \int\limits_0^{\frac{2 \pi}{\omega}} n(t) \sin(\omega t) dt \\
& \approx \frac{1}{2} \int\limits_0^{\frac{2 \pi}{\omega}} \cos(\phi) dt \\
& = \frac{\pi}{\omega} \cos(\phi)  \\
\end{align} $$

and

$$ \begin{align}
b & = \int\limits_0^{\frac{2 \pi}{\omega}} x(t) \cos(\omega t) dt \\
& = \int\limits_0^{\frac{2 \pi}{\omega}} \left(\sin(\omega t + \phi) + n(t) \right) \cos(\omega t) dt \\
& = \int\limits_0^{\frac{2 \pi}{\omega}} \sin(\omega t + \phi) \cos(\omega t) dt + \int\limits_0^{\frac{2 \pi}{\omega}} n(t) \cos(\omega t) dt \\
& = \frac{1}{2} \int\limits_0^{\frac{2 \pi}{\omega}} \sin(\phi) dt + \frac{1}{2} \int\limits_0^{\frac{2 \pi}{\omega}} \sin(2\omega t + \phi) dt + \int\limits_0^{\frac{2 \pi}{\omega}} n(t) \cos(\omega t) dt \\
& \approx \frac{1}{2} \int\limits_0^{\frac{2 \pi}{\omega}} \sin(\phi) dt  \\
& = \frac{\pi}{\omega} \sin(\phi)  \\
\end{align} $$

in both expressions (on the bottom line) we know the middle integral is zero and expect the integral on the right to tend toward zero because $n(t)$ is bipolar and ain't correlated with anything.

i think, then, that you'll find the phase $\phi$ to be:

$$ \phi = \arg\left\{a + jb \right\} $$

which is, if $a>0$

$$ \phi = \arctan\left(\frac{b}{a} \right) $$

otherwise, you might have to add or subtract $\pi$ to that expression.


FFT of any sinusoidal signal -- how many point fft to take?

Many of your questions are answered elsewhere, so I added additional links below but to summarize in one place relative to your questions I offer the explanations below that together with the links should answer your questions:

Required Sampling Frequency:

If the signal is real, then yes the sampling frequency must be greater than twice the highest frequency for a signal bandwidth that extends down to DC. However if the signal is complex then the sampling frequency should be greater than the highest frequency for same case. This is because a complex signal has a unique spectrum over the Nyquist frequency range of $\pm f_s/2$ where $f_s$ is the sampling rate, while for a real signal the negative frequencies are the complex conjugate of the positive frequencies (identical magnitude, opposite phase).

Note that prior to sampling an analog anti-alias filter is required on the signal, otherwise energy if present in higher Nyquist zones (frequencies above the Nyquist frequency range defined above) will alias into the Nyquist frequency range being observed.  (see this post which explains how aliasing occurs: Aliasing after downsampling). It is not feasible to implement a brick-wall filter, so additional margin to the sampling rate is usually added to allow for the characteristics of the anti-alias filter. For this you need to understand the properties of the signal signal of interest prior to sampling (as once sampled there is no way to distinguish an alias from a signal in the first Nyquist zone!). For example in your case, where you are interested in observing the signal out to 20 Harmonics of 60 Hz: The anti-alias filter should pass the 20th harmonics without distortion and then reject the higher harmonics so that you can minimize the sampling rate needed. You would then choose the sampling rate depending on your filter implementation such that no energy folds in to distort the 20th harmonic.

How Many Point FFT is Needed:

This has been answered several times, so see the links given below for further details and caveats, but the main point is that the frequency resolution of each "point" is 1/T Hz where T is the length of your signal in seconds.  Further the bins in the DFT (the FFT is just an algorithm that computes the DFT) as typically returned extend from DC to one bin less that $f_s$. Due to aliasing the DC is equivalent to $n f_s$ where n is any integer, and likewise for all other frequencies in the DFT, so the DFT that extends frmo DC to $f_s$ equivalently extends from DC to one bin less than $f_s/2$ and then from $-f_s/2$ to one bin less than $f_s$ (see fftshift in Matlab). As detailed in the links given, due to the relationship of the sampling rate and the total length in time T of your signal, each bin will have a frequency resolution of $1/f_s$ when no windowing is used.

Note importantly that each bin of the DFT is the result of a Sinc filter, so if your sampling rate is not an integer multiple of 60 Hz, you will also have the effect of "scalloping loss" that will distort the measurement of each harmonic. (see the "Bank of Filters" explanation at this link: Intuition for sidelobes in FFT). By choosing a sampling rate that is an integer multiple of 60 Hz, each 60 Hz harmonic will be properly centered in each bin.

Suggestion for your specific application:

With relationship to choosing the sampling rate and number of bins for your specific application, if the 60 Hz and its harmonics are the dominant signal and you don't have other signal energy of significance then there is a simple solution that will work quite elegantly in accurately estimating each harmonic with minimum measurement distortion. As described above, choose a sampling rate that is an integer multiple of 60 Hz, high enough to allow for the anti-alias filter used, and then choose the number of points such that each bin passes a harmonic of 60 Hz and rejects all the other harmonics. This means the number of samples in the FFT is $f_s/60$ where $f_s$ is an integer multiple of 60 Hz > 2400 Hz (assuming your signal is real). For example 4800 Hz may be a good choice as you suggest in your options (allowing for a fairly relaxed anti-alias filter!), and in this case $N =4800/60= 80$ points.  So in this case $f_s=4800 Hz$ and with 80 points each point is 60 Hz, so 0, 60, 120, 180 etc centered on each harmonic and perfectly nulling the other Harmonics. The actual magnitude measured for each point is the total area under the Sinc function frequency response for each bin, so if there is other energy present at non-harmonic frequencies, this would distort the magnitude measured in each bin. But if the harmonics are dominant, this would be an excellent and simple approach.  This is somewhat similar to the objective in this link in removing harmonic signals (Different way to separate a particular frequency from a signal) as each bin in the DFT is a rotated moving average filter just as done in this referenced link when only one bin is of interest. Note from that link how every other harmonic is nulled!

Please see these other responses which will further help give insight:

Specific frequency resolution

FFT frequency resolution

Number of FFT points required for a specific frequency resolution for an oversampled signal

FFT and number of samples relations

What happens when N increases in N-point DFT


Effect of Hanning Window on phase of FFT

It't really more of a question of what the rectangular window is doing to your phase relationships, because it doesn't deemphasize the "edge effects" that occur at the beginning and end of your block of samples, which is the whole point of using a non-rectangular window.

Try this experiment: Do a number of FFTs with lengths ranging from 59*48000 samples to 60*48000 samples, using both window types. You'll find that the Hanning window gives much more consistent results than the rectangular window does. The two results will agree best when the length is close to an integer number of beats.


Determining channel frequency response from measured data via IFFT

Without seeing your code directly, my guess is that you're doing something like:

N = length(x);
X = fft(x,N);
Y = fft(y,N);
H = Y./X;
h = ifft(H);

The problem with doing that is this operation performs circular convolution over a length of N.  The trouble is, the "true" length of the convolution of y with 1/X's impulse response will be of length 2N - 1 (assuming x and y are both of length N).

So you are getting time domain aliasing.

The solution is probably to set

N = length(x) + length(y) - 1

or greater.


PSD (Power spectral density) explanation

You are right, PSD has to do with calculating the Fourier Transform of the power of the signal and guess what.....it does. But first let's look at the mathematical relationship between the PSD and the autocorrelation function.

Notations:

Fourier Transform:
$$ \mathcal{F}[  x(t)] = X(\omega) = \int_{-\infty}^{\infty} x(t)e^{-j\omega t}dt $$
(Time) Auto-Correlation Function:
$$ R(\tau) = x(\tau) * x(-\tau) = \int_{-\infty}^{\infty} x(t)x(t + \tau)dt $$

Let's prove that the Fourier Transform of the Auto-Correlation function does indeed equal to the Power Spectral density of our stochastic signal signal $x(t)$.

$$ \mathcal{F}[  R(\tau)]= \int_{-\infty}^{\infty} R(\tau)e^{-j\omega \tau}d\tau $$
$$   = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} x(t)x(t + \tau) e^{-j\omega \tau} dtd\tau $$
$$   = \int_{-\infty}^{\infty}x(t) \underbrace{\int_{-\infty}^{\infty} x(t + \tau) e^{-j\omega \tau} d\tau}_{\mathcal{F}[x(t + \tau) ] =  X(\omega)e^{j\omega t}} dt $$
$$   = X(\omega) {\int_{-\infty}^{\infty}x(t)e^{j\omega t}  dt}    $$

$$   = X(\omega)X^*(\omega)= |X(\omega)|^2 $$

What does it all mean?
Note: This explanation is a bit "hacky". But here it goes

The Fourier transform tells us the spectral components of a signal. In our case, the signal is Stochastic; So, trying to calculate the spectral components of the signal will be pointless because, for every realisation of the random process, you will have different expressions for $ \mathcal{F}[x(t)]$.

What if you take the Expected Value of the Fourier transform then? This wouldn't work. Let's take a zero mean signal for example.

$$ \mathbb{E}\{ \mathcal{F}[x(t)] \} = \mathcal{F}[\mathbb{E}\{ x(t) \}] = 0$$

Instead, what if you take the Fourier transform of the square of the signal.
$$ \mathbb{E}\{ \mathcal{F}[x^2(t)] \} = \mathcal{F}[\underbrace{\mathbb{E}\{ x^2(t) \}}_{\text{Av. Power of the Signal}}] $$

The autocorrelation function is essentially the $P(t) $ which you were alluding to.

References:

[1] Communications 1, P-L. Dragotti, Imperial College London

[2] White Noise and Estimation, F. Tobar [Unpublished Report]


i see aliasing happening when the image is captured with a sensor that cannot capture higher frequencies.

Exactly!

If i capture a image, is it always aliased?

Not when the image content is sufficiently bandlimited.
That's often the case because of the physics behind things, or simply because the motive isn't higher in frequency.
(Physical effects include motion blur, abberation, depth of field/ring of confusion/bokeh, wave effects on apperture edges, lens imperfections, air density variations and many more. I'd recommend you look for tutorials on how to take really sharp macro and tele photos to learn what effects can reduce your image bandwidth.)
But, sure. Go photograph a fine grid, and you see a Moiré pattern. By the way, your eyes have limited resolution, too, although I'd presume the brain largely adapts that away.

Other types of aliasing are as far as i understand not really aliasing. e.g if a screen has not many pixels, the image is aliased.

You call it aliasing, it aliases things, you describe it as aliasing: It's really aliasing.


Correlation of discrete time periodic sequences

Energy signals (such as an exponentially decaying pulse) are bounded in total energy as time goes to infinity in which case you can use the first definition. Power signals (such as a sine wave) have finite energy over finite time intervals but infinite energy for all time (power is energy/time) so the first equation will not converge in that case.


MFSK vs FSK: Why BER decreases increasing the symbols

(I assume maintaining the same symbol duration)

nope! That is usually not the case, because it would be unfair – if you send the same number of symbols per second, you'd get twice as many bits per second using 4-FSK rather than 2-FSK.

Therefore, the comparison is typically done with identical Bit Energy rather than identical Symbol duration (which would be, for constant envelope signals, be equal to identical symbol energy).

Hence, with higher modulation order, you get less symbols per bit, and hence, for constant $E_b$ more $E_{sym}$.

Now, errors happen on a symbol level – and a Symbol Error Rate of $r$ doesn't imply a BER of $r$ – in fact, a wrong symbol has at most its worth in bits in wrong bits, but usually less (see: Gray code), so $\text{BER}\le\text{SER}$ generally holds, and one can in an uncorrelated noise scenario be almost certain that $\text{BER}<\text{SER}$.


Calculation of cosine in frequency domain instead of calculatin in time-domain followed by a FFT

If you have a pure tone signal defined by:

$$ S[n] = M \cos( \alpha n + \phi ) $$

Where $f$ is the frequency in cycles per frame and $\alpha$ is in radians per sample:

$$ \alpha = f \cdot \frac{ 2\pi }{ N } $$

And use a $ 1/N $ normalized DFT:

$$ Z[k] = \frac{ 1 }{ N } \sum_{n=0}^{N-1} { S[n] e^{ -i\beta_k n } } $$

Where:

$$ \beta_k = k \cdot \frac{ 2\pi }{ N } $$

The most concise formula for calculating the DFT bin values is:

$$  Z[k] = \frac{ M }{ 2N } \left[ \frac{ Ue^{ i\beta_k }  - V }{ \cos( \alpha ) - \cos( \beta_k ) } \right] $$

Where:

$$ U = \cos( \alpha N + \phi ) - \cos(\phi) $$

$$ V = \cos( \alpha N -\alpha + \phi ) - \cos( -\alpha + \phi ) $$

Notice that there is a pluggable discontinuity at the bin value of a whole integer frequency.  When the frequency is very close to an integer value there is a more computationally robust form of the formula:

$$ Z[k] = \frac{ M }{ 2 }
\cdot \frac{ \sin \left( \frac{ \delta N }{2} \right) }{ N \sin \left( \frac{ \delta }{2} \right) }
\cdot \left[ \frac{ \sin \left( \frac{ \delta N }{2} + \phi \right) e^{ i\beta_k }  - \sin \left( \frac{ \delta N }{2} - \delta - \beta_k + \phi \right) }{ \sin \left( \frac{ \delta  }{2} + \beta_k \right) } \right] $$

Where:

$$ \delta = \alpha - \beta_k = ( f - k ) \cdot \frac{ 2\pi }{ N } $$

The derivations of these formulas can be found in my blog articles "DFT Bin Value Formulas for Pure Real Tones" and "An Alternative Form of the Pure Real Tone DFT Bin Value Formula".

I realize this is a reply to a very old post, but Community dredged it up and finding these formulas, and the subsequent inverses which give exact frequency formulas is what got me into DSP.

Ced


Conjugate symmetry of real-coefficient filters in Oppenheime's Discrete Time Signal Processing

The Fourier transform of a real function $x[n]$ is conjugate symmetric. That is,
$$X^*(\omega)=X(-\omega)$$
This is easy to observe. From the definition of FT: $$X(\omega)=\sum_{n=-\infty}^{+\infty}x[n]e^{-j\omega n}$$
Conjugating both sides:

$$\begin{align}
X^*(\omega)=&\sum_{n=-\infty}^{+\infty}x^*[n]\left(e^{-j\omega n}\right)^*\\
&=\sum_{n=-\infty}^{+\infty}x[n]e^{j\omega n}\\
&=X(-\omega)
\end{align}$$

A direct consequence of this property (that you can easily verify) is that

  1- The real part and the magnitude of the Fourier transform of a real valued function are even functions of $\omega$. That is, $\Re\{X(\omega)\}=\Re\{X(-\omega)\}$ and $|X(\omega)|=|X(-\omega)|$.

  2- The imaginary part and the phase of the Fourier transform of a real valued function are odd functions of $\omega$. That is, $\Im\{X(\omega)\}=-\Im\{X(-\omega)\}$ and $\angle X(\omega)=-\angle X(-\omega)$.

Based on the above explanation, for a real-coefficient filter $h$, the magnitude of the frequency response is an even function.

$$\begin{align}
|H(e^{j(\omega-\pi)})|&=|H(e^{j(-\omega-\pi)})|\\
&=|H(e^{j(-\omega-\pi+2\pi)})|=|H(e^{j(\pi-\omega)})|
\end{align}$$
and the second line is because DTFT is $2\pi$-periodic. Hence we can add $2\pi$ and the result does not change.


Understanding an alternative representation for the difference equation of a FIR filter?

Given that
$$
u[n-k_1] = \sum_{k=k_1}^\infty \delta[n-k]
$$

you can take the difference between two steps with different shifts ($k_1 < k_2$) and obtain:
$$
u[n-k_1] - u[n-k_2]
  = \sum_{k=k_1}^{k_2-1} \delta[n-k]
$$

Then the difference

$$
3u[n-2] - 3u[n-6]
  = 3\left(\delta[n-2] + \delta[n-3] + \delta[n-4] + \delta[n-5]\right)
$$

which appears in your second equation can be replaced to give:
$$
\begin{align}
  y[n]
    &= 2\delta[n] + 6\delta[n - 1] + 3u[n - 2] - 2\delta[n - 4] + \delta[n - 5] - 3u[n - 6] + 2\delta[n - 7] \\
  &= 2\delta[n] + 6\delta[n - 1] - 2\delta[n - 4] + \delta[n - 5] + 2\delta[n - 7]
     + 3u[n-2]-3u[n-6] \\
  &= 2\delta[n] + 6\delta[n - 1] - 2\delta[n - 4] + \delta[n - 5] + 2\delta[n - 7] + 3\left(\delta[n-2] + \delta[n-3] + \delta[n-4] + \delta[n-5]\right) \\
  &= 2\delta[n] + 6\delta[n - 1] + 3\delta[n-2] +3\delta[n-3]+ \delta[n - 4] + 4\delta[n - 5] + 2\delta[n - 7]
\end{align}
$$

which corresponds to your first expression.


Exact definition of soft bits in digital receiver

Soft bits or more generally speaking, soft metrics, are usually taken at the output of the demodulator (e.g. the output of a matched filter sampled at $T_b$). A soft bit can be therefore a real value, but it is usually quantified (which is required if a DSP comes after).

A hard bit is obtained by performing a hard decision on a soft metric. So for instance, if $x[kT_b]$ is the output of the matched filter when a BPSK signal is applied, then a hard-bit can be obtained by simply comparing $x[kTb]>0$, which is $0$ or $1$.

This figure from wikipedia is self-explanatory.

As you see, the output of the matched filter is actually a continuous signal. Some part of your receiver (which is not included in the figure) must perform also timing estimation in order to determine at which instant the output of the matched filter must be sampled. Those samples are the red points in the figure, and might be considered as soft bits. So the higher the magnitude of those samples, the more information they provide about what the original transmitted symbol was. That is the advantage of having soft-bits: they don't only tell you what symbol has been probably sent, they also give an idea about that probability.


Illustration of time domain aliasing cancellation (TDAC, MDCT, lapped transform)

Edit: I have recently created two Jupyter Notebooks that illustrate this behaviour and let you play around with some actual matrices and actual signals.

I find understanding MDCT easiest if we define our transforms as matrix operations.

DFT

In case of a DFT such a matrix would look like

$$
F_M = \left( \sqrt{\frac{1}{M}} e^{-2 \pi i k n / M} \right)_{k,n=0 \dots M-1}
\in \mathbb{C}^{M \times M}
$$

$F_M \approx$

and our transform would become a matrix multiplication

$$
X = F_M x
\qquad{} x,X \in \mathbb{C}^M
$$

where $X$ and $x$ are column vectors of size $M$. The inverse of the matrix $F$ exists, so synthesis is possible

$$
x = F^{-1}_M X
\qquad{} x,X \in \mathbb{C}^M
$$

which means (note we are using the unitary matrix definition here)

$$
F_M F^{-1}_M = I
$$

Wide DFT

If we now just defined a "wide DFT" matrix to be larger in one dimension, lets say $M \times 2M$.

$$
G_{2M} = \left( \sqrt{\frac{1}{M}} e^{-2 \pi i k n / M} \right)_{k=0 \dots M-1 \quad{} n=0 \dots 2M-1}
\in \mathbb{C}^{M \times 2M}
$$

this is equivalent to having longer basis functions.

$G_{2M} \approx$

This means for $2M$ samples $x$ we get $M$ coefficients $\hat{X}$.

$$
\hat{X} = G_{2M} x
\qquad{} x \in \mathbb{C}^{2M}
\quad{} X \in \mathbb{C}^M
$$

Doing that, you might also notice is that there are now repititions in our matrix as our basis function $e^{-2 \pi i k n / M}$ starts to repeat after $M$ samples. It looks like two DFT matrices stacked next to each other, which we can represent as a folding matrix followed by DFT

$$
G_{2M} = F_M
\begin{bmatrix}
I_M & I_M
\end{bmatrix}
$$

$G_{2M} \approx$  $=$  $\times$

Where $I$ is the identity matrix. Obviously an inverse of $\begin{bmatrix}I_M & I_M\end{bmatrix}$ doesn't exist, so an inverse of $G$ doesn't exist either.

If we now look at the complete transform

$$
\hat{X} = F_M
\begin{bmatrix}
I_M & I_M
\end{bmatrix} x
\qquad{} x \in \mathbb{C}^{2M}
\quad{} X \in \mathbb{C}^M
$$

we see that always the sum of exactly two samples of $x$ are being transformed. This is called time domain aliasing.

We can do this for any size and the folding matrix will just repeat, e.g.

$$
G_{4M} = F_M
\begin{bmatrix}
I_M & I_M & I_M & I_M
\end{bmatrix}
$$

$G_{4M} \approx$  $=$  $\times$

and we will end up with the sum of more and more samples being transformed together, producing more and more aliasing.

DCT-IV

We can define DCT-IV as a matrix, too.

$$
C_M = \left( {\sqrt{\frac{2}{M}} \cos\left[\frac{\pi}{M} \left(k+\frac{1}{2}\right) \left(n+\frac{1}{2}\right)\right]} \right)_{k,n = 0 \dots M-1}
\in \mathbb{R}^{M \times M}
$$

$C_M \approx$

Interestingly, the matrix is orthogonal and its own inverse

$$
C_M = C^{-1}_M = C^T_M
$$

Wide DCT-IV

Again, we can define a wide DCT-IV of size $M \times 2M$ by just extending the basis functions.

$D_{2M} \approx$

However, this time the repititions in the basisfunction are different, and thus the folding matrix is:

$$
D_{2M} = C_M
\begin{bmatrix}
I_M & -J_M
\end{bmatrix}
\in \mathbb{R}^{M \times 2M}
$$

$D_{2M} \approx$  =  $\times$

where $J_M$ is the $M \times M$ antidiagonal matrix. If we repeat this further, we get

$$
D_{4M} = C_M
\begin{bmatrix}
I_M & -J_M & -I_M & J_M
\end{bmatrix}
\in \mathbb{R}^{M \times 4M}
$$

$D_{4M} \approx$  =  $\times$

after which our folding matrix will start to repeat. You will discover this is exactly the repetition pattern of the image on Wikipedia

We will need this pattern, so lets define it

$$
Q =
\begin{bmatrix}
I_M & -J_M & -I_M & J_M
\end{bmatrix}
\in \mathbb{R}^{M \times 4M}
$$

MDCT

Analysis

$$
X_k = {\sqrt{\frac{2}{M}} \sum^{2M-1}_{n=0} h(n) x_n \cos\left[\frac{\pi}{M} \left(k+\frac{1}{2}\right) \left(n+\frac{M}{2}+\frac{1}{2}\right)\right]} \qquad{} k = 0 \dots M-1
$$

If we look at the definition of the MDCT we note that we are transforming $n = 0 \dots 2M$ samples to $k = 0 \dots M$ coefficients, so we are using a "wide matrix" of size $M \times 2M$.

$$
E_{2M} = \left( {\sqrt{\frac{2}{M}} \cos\left[\frac{\pi}{M} \left(k+\frac{1}{2}\right) \left(n+\frac{M}{2}+\frac{1}{2}\right)\right]} \right)_{k = 0 \dots M-1 \quad{} n=0 \dots 2M-1}
\in \mathbb{R}^{M \times 2M}
$$

$$
X_k = E_{2M} x
\qquad{} x \in \mathbb{R}^{2M}
\quad{} X \in \mathbb{R}^M
$$

But also, if we compare the basis function to DCT-IV, we notice they are almost the same, but with an offset of $M/2$ for $n$. This means we are using a "shifted view" of DCT-IV matrix, shifted to the right:

$$
E_{2M} = \left( {\sqrt{\frac{2}{M}} \cos\left[\frac{\pi}{M} \left(k+\frac{1}{2}\right) \left(n+\frac{1}{2}\right)\right]} \right)_{k = 0 \dots M-1 \quad{} n=\frac{M}{2} \dots \frac{M}{2}+2M-1}
\in \mathbb{R}^{M \times 2M}
$$

We can express such a shifted view as a cyclic permutation matrix $P\in \mathbb{R}^{4M \times 4M}$, but only if our basis function was cyclic, too. We know $D_{4M}$ is expanded to the point where it becomes cyclic, so we can express $E_{2M}$ using $D_{4M}$ and using an extraction matrix and a permutation matrix:

$$
E_{2M} =
D_{4M}
P
\begin{bmatrix}
I_{2M} \\
0
\end{bmatrix}
$$

$E_{2M} \approx$  $=$  $\times$  $\times$

but since we know how to represent wide DCT-IV by a regular DCT-IV like

$$
D_{4M} = C_M Q
$$

we can write

$$
E_{2M} = C_M Q P
\begin{bmatrix}
I_{2M} \\
0
\end{bmatrix}
$$

$E_{2M} \approx$  $=$  $\times$  $\times$  $\times$

We don't need all of those shift matrices individually, so let's multiply them together into one matrix $R$:

$$
R = Q P
\begin{bmatrix}
I_{2M} \\
0
\end{bmatrix}
$$

$E_{2M} \approx$  $=$  $\times$

And end up with

$$
X = C_M R x
$$

This means MDCT is a permutation and superposition and then a DCT-IV of our input samples.

Synthesis

For reconstruction we need to form a linear combination of basis vectors, which can be written as

$$
\hat{x} = (C_M R)^T X = R^T C_M^T X
$$

And substituing our forward transform for $X$ we get

$$
\hat{x} = R^T C_M^T C_M R x
$$

we know

$$
C_M^TC_M = I_M
$$

but

$$
R^T R \neq I_M
$$

instead

$R^T R \approx$

But if we now overlap+add two frames with just the right hopsize

we get

perfect reconstruction in the overlapping part.

This means, if we first transformed our signal in windows with half a window overlap, all the aliasing produced in this overlap will exactly cancel out, making the entire transform perfectly reconstructing.

Caveats/Further reading

Windowing

All of this is obviously disregarding windowing. But if you look at the permutation matrix and look at what sample is superpositioned with what other sample it should be at least somewhat intuitively visible where the Princen Bradley condition comes from.

Nicer Permutation

Some authors prefer the MDCT to be defined as

$$
X_k = {\sqrt{\frac{2}{M}} \sum^{2M-1}_{n=0} h(n) x_n \cos\left[\frac{\pi}{M} \left(k+\frac{1}{2}\right) \left(n+\frac{M}{2}-\frac{1}{2}\right)\right]} \qquad{} k = 0 \dots M-1
$$

(note the minus). If you then solve all of the remaining math, you will find a permutation matrix that looks like

$\hat{R} \approx$

as opposed to

$R \approx$ .

This doesn't change any of the remaining explanations (TDAC still works etc.), but you can see that the samples in the center of your frame are now 'nicer': they're positive valued and transformed forwards, as opposed to negative valued and transformed backwards.

Math for overlap add

If you feel that last bit became a bit too "hand wavey", that's because it is. The traditional math of transforms or matrices starts to break down if you start talking about overlapping transforms. A very neat way of dealing with them are Polyphase Matrices in the z-Domain. They're hard to wrap your head around at first though.


This line is wrong:

xr(k)=xr(k) + xb(p)*sinc(fs(i)*(k-t2(p)));

The correct one is:

xr(k)=xr(k) + xb(p)*sinc( ((k-1)*Tr-(p-1)*Ts)/Ts );

Where $Ts$ is the original sampling period and $Tr$ is the reconstruction sampling period. In your code $Ts$ is given by $T_s = 1/500$ for sampling frequency of $500$ Hz. And $T_r = 10^{-4}$ set by your continuous plot period.

Note that you will still have slight errors due to the fact that the input signal is not being perfectly bandlimited.

Here is the resulting plot for $Fs=500$ Hz, $f_0 = 200$ Hz and $T_r=10^{-4}$ seconds:

it's ok now.


Why is level of power spectrum dependent on FFT resolution?

The observed phenomenon is simply due to the fact that when you set window length to 100 it is as if you were doing fft with 100 samples instead of 1000. So Nw=100 similar to NFFT = 100.

And why does FFT size affect results (level of noise VS level of tone)? Good question

Two types of spectral densities:

1) PWELCH-like PSD: "PSDW"

A tone grows in dB as FFT size increases. As in a true PSD the tone is a Dirac's delta for NFFT->inf. Larger NFFT means smaller dF, and larger A, in such a way that A * dF = constant
A unit power noise stays at 0dB.

This is really a [Watt/Hz] spectral density and total power given by:

totalPower = trapz(F,PSDW)  =  trapz(fs * F,PSDW/fs)  = sum(PSDW)/NFFT

2) OSCILLOSCOPE-like PSD: "PSDO", PSDO = PSDW/NFFT

A tone remains the same dB as FFT size increases
A unit power noise becomes less dBs as NFFT increases.

This is a [Watt/Hz/NFFT] measure, with power:

totalPower = NFFT * trapz(F,PSDO) =  NFFT*trapz(fs * F,PSDO/fs)  = sum(PSDO)


Use predictive coding. If your predictor manages to guess current position based on previous recorded positions with reasonable accuracy — you can omit the recording of current position. Same predictor will also help you to reconstruct the original path.

What predictor you should use is based on the tracked object's nature. You may use a Kalman filter for example. Or you may use something really simple such as 'assume object always stays at the same place. If we are moving — record it'. Or 'assume object always travels along the straight line', etc.

If your predictor will do a good job, you can additionally use Huffman coding — small corrections of prediction errors will have shorter bit sequences, significant corrections — longer bit sequences. It will save you a huge amount of space with near-to-perfect reconstruction.


Difference between Digital signal and Discrete signal

Suppose you have a continuous time analog signal. It is continuous in both time and amplitude. Now when you sample it ,you get discrete samples every Ts seconds. Now you have discrete samples(discrete in time) each of which can take continuous value( in amplitude). This is normally referred to as discrete signal(discrete in time but continuous in amplitude).

Now further when you take this discrete signal and quantise it i.e. say assign each sample which takes continuous amplitude to one of the N discrete quantisation levels of a quantiser, Then your total signal now is a digital signal. So a digital signal is discrete in time and discrete in amplitude.


Using MIMO with DS-CDMA and spread spectrum

Assuming one antenna transmits one symbol per time unit, then 16 symbols require 4 time units to be out. Then it is simply that

r_1 = H_1 * x_1(1:4)
r_2 = H_2 * x_1(5:8)
r_3 = H_3 * x_1(9:12)
r_4 = H_4 * x_1(13:16)

If channel H is fixed during these 4 time units,

[r_1 r_2 r_3 r_4] = H * [x_1(1:4) x_1(5:8) x_1(9:12) x_1(13:16)];

or

r = reshape(H*reshape(x_1,4,[]),[],1);

This can be generalized to any MIMO size and to any spreading size.


Cooley-Tukey FFT: You don't have to zeropad to a power of 2?

Yes.  If you can factor length N, Cooley-Tukey will give you an DFT computation that uses less multiplies than an unfactored DFT.  A power of 2 factorization requires less lines of text or shorter equations to explain.  But any small prime factors allows a straight forward extension to the power-of-2 explanation or proof of the algorithm.

Many newer FFT libraries (iOS Accelerate for one example) allow efficient composite FFT lengths with factors of 2,3,5 and maybe more, thus requiring less padding to get to the next useful composite length.


Intuition behind cepstrum analysis terminology

in filtering we try to separate two signals from each other. if two signals are added together this is the normal filtering. but these two signals could be multiplied or convoluted with each other. considering there are huge scientific background on this kind of filtering ( for additive signals) we try to transform other problems to this kind.

if signals multiplied with each other, we take the logarithm of their multiplication to convert it two additive signal. then we use filtering methods for additive signals.

if they were convoluted with each other, we know their spectrum is multiplied, so by taking the logarithm of their spectrum, we reach two additive spectrum in frequency domain. now we treat these two additive spectrum in frequency domain as two additive signals and apply the filtering methods for additive signals. considering we work on frequency domain signals if we use normal filtering terminology there would be some confusion e.g. if we say we use a low-pass filter someone might think we apply a LPF on our original time domain signal instead of their frequency domain counter part, so in order to avoid these confusions some new terminology is used.

now if we take Fourier transform of these frequency domain signals we go to time domain, so a low-pass filter (for frequency domain signals) becomes a short-pass filter which passes the signal at times before a cut-off time and high-pass filter becomes long-pass filter which passes the signals after a cut-off time.


Trying to run an implementation of lsqcurvefit from the Optimization Toolbox from Matlab in python

As said in the comments this question is better aligned with stack overflow.
Answer
The problem is that curve_fit assumes ydata = f(xdata, *params) + eps, and you assumed ydata = f(params, xdata) +  eps. If I understand your problem correctly your function could be implemented as
def myfun(xdata, a, b, c):
  return np.exp(a)./ xdata**np.exp(b) - c


Streaming windowed sinc interpolation/resampling: trying to understand a Rust implementation

This looks partially similar to the BASIC code in my Q&D resampling cookbook here:  http://www.nicholson.com/rhn/dsp.html#3  which computes both the Sinc and the window (von Hann) on the fly (per sample).
But it appears that the .rs code is not determining the lowpass cutoff frequency (phi upper and lower) correctly.  It needs to know whether the resampling is up or down in frequency, so that the Sinc computation can use something below half the lesser of the two frequencies.
Computation of one new resampled sample at a time works just fine, as long as the rate change is fixed (or increases slightly), and the initial and ending transients (related to the Sinc window width (or depth?)) are ignored (or accepted as transients).


Going from output of D/C to input of D/C

You won't get the discrete $y_d[m]$ by re-arranging the equations, but by sampling the $y_c(t)$ at integral multiples of $T$. When you sample the RHS in equation1 exactly at $t=mT$ instances, you get $y_d[m]$. You need to remember that $sinc(m-n) = 1$ for $m=n$ and $0$ otherwise, and (BTW in equation 1, in $sinc$ argument $\pi$ won't be there). Anyway, take RHS and evaluate the expression at $t = mT$ where $m$ are integers.
$$y_c(t) \big|_{t=mT} = \sum^{\infty}_{n=-\infty}y_d[n] \, \mbox{sinc} \left( \frac{t-nT}{T} \right) \Big|_{t=mT} = y_d[n]$$


The impulse response, $h(t)$, of a system

firstly do partial fraction for H(s) and then use transform tables of Laplace Transform. Using inverse laplace and tables you must be able to find impulse response. I think you can use j^2 instead of -1 and find complex roots for denominator. Then, using transform pairs you will find a complex solution.


Frequency Domain Interpolation: Convolution with Sinc Function

The "discrete sinc kernel" is also known as a Dirichlet Kernel.  It looks like a you got a 2D (u and v) going there.  I'm not that familiar with that, but we just did a pretty good number on the 1D case.

This post explains the inverse DFT can be interpreted as "epicycle drawing":

How to get Fourier coefficients to draw any shape using DFT?

This post then follows up with the relationship between the sampling, epicycle drawing, upsampling, and the normalized sinc function:

Absolute convergence of periodic sinc interpolation

The conclusions are:

Suppose a continuous band limited function (with band limit $L$) is sampled, then reconstructed from the samples.

1) If L < N/2 the function can be unambiguously reconstructed.  The results from an infinite series sinc reconstruction (Whittaker-Shannon) and taking a DFT, zero padding at the Nyquist and then the inverse DFT to upsample, will produce the same path.

2) If L = N/2, the Nyquist bin must be split between N/2 and -N/2 for the sinc reconstruction to match the DFT reconstruction.  It may or not may not match the original function, but it will match at all the sample points.

3) If L > N/2, and N is even, the Nyquist bin must still be split for the sinc reconstruction to match the DFT reconstruction, but neither will match the original function, yet both will still match at all the sample points.

The relationship between the sinc function and the Direchlet Kernel is this:

1) The sinc function is the limit of the Dirichlet kernel as the sample count goes to infinity

2) For odd N, the Dirichlet kernel is an infinite sum of sinc functions.  For even N, it is an adjusted one.

See the posts for details and discussion.

What I forgot to overtly mention is that the Dirichlet Kernel formula is directly derived from, and equivalent to, doing the FFT, zero pad at Nyquist (split Nyquist on even N) then inverse DFT for resampling or coefficients for Fourier Series.  Your link in the comment is about the reverse case, but since the forward DFT and inverse DFT are mirror image twins (especially when $1/\sqrt{N}$ normalized), so what gets said in one direction has to be true in the other.  This means if you zero pad a signal, and have an even number of points, one of the outside points needs to be halved and the other half appended at the other end before zero padding outward, leaving your data centered.

I already explained the relationship between the Sinc and the Dirichlet, just above.  The math can be found in the second link in the section of my answer that has $y_m = Y(w)$ in it.  Towards the bottom.


Calculating phase response of maximum phase filter using Hilbert Transform

Note that a Hilbert transform relationship for maximum phase systems can only exist if a maximum phase system is defined to be an anti-causal stable system with an anti-causal stable inverse (cf. these MIT course notes, p.4). Considering only systems with rational transfer functions, this means that not only the system's zeros but also its poles must be in the right half-plane (for continuous-time systems) or outside the unit circle (for discrete-time systems), respectively. The reason for this is that for the Hilbert transform relationships to hold, $\ln(H(s))$ (or $\ln(H(z))$) must be analytic, either for $\text{Re}\{s\}\ge 0$ (or $|z|\ge 1$) (minimum-phase), or for $\text{Re}\{s\}\le 0$ (or $|z|\le 1$) (maximum-phase). Since both poles and zeros of $H(s)$ cause singularities in $\ln(H(s))$ we need to make sure that for maximum phase systems there are neither poles nor zeros in the left half-plane $\text{Re}\{s\}\le 0$ (or inside the unit circle for discrete-time systems).

Note that the above definition of a maximum phase system is different from another common definition which says that a maximum phase system is a causal and stable system with all its zeros in the right half-plane (or outside the unit circle), but - since it's causal and stable - all its poles in the left half-plane (inside the unit circle).

According to the first definition above, a maximum phase system is obtained from a minimum phase system by time-reversal, which leaves the magnitude of the frequency response unchanged but inverts the sign of the phase. Consequently, for a maximum phase system according to above definition, we simply get the following Hilbert transform relationship between phase and magnitude:

$$\theta(\omega)=\mathcal{H}\{\ln(A(\omega))\}\tag{1}$$

Note that the sign inversion in the Hilbert transform relationship between magnitude and phase can also be seen from the fact that a minimum phase system has a causal cepstrum (from which the Hilbert transform relationship follows), and a maximum phase system has an anti-causal cepstrum, resulting in the same Hilbert transform relationship between magnitude and phase but with an inverted sign.


Is this normal ?

Yes.
It follows the mathematical definition of the DFT.


How to design a variable slope highpass filter

The order of a filter determines its slope. Normally orders are positive integers (1, 2, 3, ...). So the slopes are fixed to multiples of approximately 3db/octave. However, taking the case of a Butterworth filter, someone has applied partial calculus to obtain fractional order Butterworth filters which will effectively allow you to chose the slope. If you do a Google search for fractional order Butterworth filters you will get back some papers on the subject. For instance here: https://arxiv.org/ftp/arxiv/papers/1210/1210.8194.pdf

If you need a different filter type you could apply the same thinking, but the maths may be more complicated for, say, an elliptic filter.


Template matching or object recognition?

That's hard. I'm not even sure I can do this myself before I actually try.
It is definitely not something that can be described in a DSP.SE answer.

About canny edge even if it is precise you will just get an image like a line drawing as an output. The point is what's next. There are many things to try the are many theory about pattern matching from the edge but it not that easy.

I will try something easier before going that path for example I will try to locate coordination and radius of those circles. (nuts, holes, etc) And see if I can identify the part using just those info. Then I'm going for the straight lines.

Please note that I know just an overview of theories in this field and not generally do image processing but you can try what I suggested using OpenCV for sure.


Because frequency modulation is a nonlinear phenomenon, there's no good general-purpose answer to your question.  You can take any specific baseband signal and figure out its modulated spectrum, but that's it.
For lack of a better technical term, FM modulation takes the frequencies in the parent signal and munges them together, spreading them out by more than the deviation frequency.  For a sine wave this is that Bessel function.  For a pair of sine waves, it'd be all the harmonics of each sine wave plus all possible combinations of intermodulation ($n f_1 + m f_2\, \forall\, n, m \in \mathbb{Z}$) -- so for the two-sine case there may be some discernible, identifiable lines in the output.  For sound with a continuous spectrum, you'll get a continuous spectrum out, but nothing you could just recognize on a spectrum analyzer.
For practical purposes, you're left with rules of thumb, inspired by the above results but limited by our inability to solve the nonlinear differential equations in any general way.  If you have a signal with bandwidth $B$, and a maximum frequency deviation of $f_D$, then the bandwidth of the transmitted signal can be limited to $2(B + f_D)$.  This isn't satisfying from a theoretical standpoint, but it's backed by over 70 years of practical experience with broadcast and voice communications systems.


Wind Turbine Vibration Analysis: Which approach?

One technique that I've seen used in determining bearing faults is using the kurtosis of the vibration signal.

You can track as a function of time what Wikipedia calls the sample excess kurtosis. This is the kurtosis that is different from the kurtosis you would see if the signal was Gaussian distributed.

The sample excess kurtosis is defined as:

$$
\frac{m_4}{m_2^2} - 3 = \frac{\frac{1}{N} \displaystyle \sum_{n=0}^{N=1} (x_n - \bar{x})^4}{\left(\frac{1}{N} \displaystyle \sum_{n=0}^{N=1} (x_n - \bar{x})^2\right)^2} - 3
$$

where $\bar{x}$ is the sample mean, $m_2$ is the sample second order moment, and $m_4$ is the sample fourth order moment.

Brüel & Kjær have a nice writeup about using kurtosis, though they use the actual kurtosis rather than the excess kurtosis.

Why?

Bearing faults tend to show up when grit or dirt or broken bearing pieces get into the rotating elements.  When that happens, the vibration signals tend to be more "spiky" as the fault-inducing material gets ground against the rolling elements and bearing races.

Other approaches

You can also look at removing what you know about the harmonic content of the signals and then looking at the remainder.  There's some stuff about that here.  However, try something simple like kurtosis first.


What is the physical meaning of coherence bandwidth in mmWave massive MIMO channel?

You are right, coherence bandwidth is the frequency domain counterpart of delay spread.

However, to achieve diversity across antennas the spacing between the antennas is important "relative to the environment".Let me explain that a bit more. In a user mobile phone you would typically find the order of wavelengths seperation between antennas sufficient to get independent fading accross them. However, the same is not true at the base station. Typically at the base station the seperation between the antennas has to be a 10s of order of wavelength. Why so? Because a user equipment near to ground sees a much richer scattering environment, which are enough to cause phase changes within a wavelength. However at the base station (usually placed high up in the air) the scattering environmental is not that rich and hence more separation is required.

So, to study the effect of independent fading with antenna separations one can use the methodology to see if other antennas are seeing the same coherence bandwidths or not. It is all part of channel profiling whenever new frequencies are deployed.


The problem you want to solve is called automatic audio alignment, or auto-synchronization of audio. Common algorithms are phase correlation, cross correlation on the raw audio waveform. People also use more "cooked" audio features like MFCCs, say each 10ms frame, and measure frame-wise similarity between the audio segments to be compared.

There are many threads on Stack overflow that show how to do this. For example: find-short-clip-of-audio-within-longer-clip-of-audio, cross-correlation-vs-fft-for-finding-phase-between-2-signals. Also some interesting repos on github, audiosync, videosync


Must OFDM signal be time-division in a multiple-transmitter and SFN system when the signal is free of spread-spectrum?

In a multi-user scenario where every transmitter uses the same center frequency and bandwidth, but transmits different data, time-division multiplex has to be applied. Otherwise, interference between the multiple transmit signals will make successful reconstruction of a specific transmit signal at the receiver virtually impossible.

However, in typical SFNs, all transmitters are part of a broadcast network and transmit identical data. Examples are Digital Audio Broadcast (DAB) and terrestrial Digital Video Broadcast (DVB-T). In such networks, the guard interval of OFDM symbols can compensate for the propagation delay between different transmitters and there is no need for time-division multiplex.

Also, I don't think CDMA and OFDM can be combined. As a multiple-access scheme, OFDM(A) allows sharing of the same frequency band by frequency multiplex while CDMA applies code multiplex.