TL;DR the x-axis is time and the y-axis is coefficients which means they have no units.

Let's go through the process of manufacturing the MFCC features while emphasizing measure units tracking. We start from a sampled audio signal $s(n)$ with amplitude units on the y-axis and time units on the x-axis.

First step is to frame the audio signal. An audio signal is statistically independent, meaning it is constantly changing and unmanageable. To make it approximately constant (statistically speaking) we frame it into $20_{ms} - 40_{ms}$ frames. This does not change the units so we are staying with seconds on the x-axis for each $s_i(n)$, where $i$ denotes the frame.
Next step would be to calculate the Periodogram estimate of the power spectrum for each frame $s_i(k)$. This represents the power distribution of the signal per frequency bin. the x-axis of the output is in units of frequency (Hz or $\frac{\pi \times rad}{sample}$ depending on normalization) y-axis will then respectively be $\frac{power}{frequency}$ ($\frac{db}{Hz}$ or $\frac{db}{\frac{\pi \times rad}{sample}}$). The intuition here is to imitate the human hearing system which is Identifies which frequency regions have more energy in them.
The human ear is more sensitive to frequency changes in the low-frequency region over the high one. In other words, it is easier for us to distinguish the $100_{Hz}$ frequency difference $400_{Hz}-300_{Hz}$ over $1400_{Hz}-1300_{Hz}$. To imitate this behavior, we use a filter bank to cumulate the energy over a number of bins (typically 26). We use triangular filters, to sum energies in a certain region. We use 26 filters and the higher the frequency the wider the filter. This leaves us with power units.

An example of the triangular filters:

On the next step, we approximate the results using Discrete Cosine Transform (DCT). The output of that process leaves us with 26 coefficients. We usually discard the lower and are left with 12-13. This means we have 12-13 coefficients per frame, without any units.
We produce these coefficients for each time frame and stack them together. as presented in the image below:

Note here, the x-axis represents different time frames and is measured in seconds. The y-axis is the discrete values for 12 coefficients. the level of each coefficient for each time frame, is represented by the color scheme.

An additional big advantage of this representation is also being a very compact representation of sound resulting in 12-13 numbers per time frame.


Calculating SNR of a complex-valued signal correctly

Just simply compute the variance of the signal and the variance of the noise. (And adjust the noise for sampling factor if necessary, see below). Variance is not a complex result and the  complex conjugate is used in the computation so should be invariant to real or imaginary content; 5+j0 and 0+j5 have the same variance; if the entire signal rotated in phase, how would that have any effect on the power?
Note computing SNR when you are creating a signal is trivial as you can completely separate the signal and noise components, scale appropriately and add. Measuring SNR on a combined signal is a little more challenging and for this purpose I recommend using correlation to measure SNR given the direct relationship between rho (correlation coefficient) and SNR.
Also pay attention to the spectral density of the noise to set the SNR properly such that it is scaled to be just the portion of the noise that is in the bandwidth of the signal.


Zernike Moments' implementation in OpenCV

I haven't used the Zernike module from PeRL's myself, but check out their implementation here.

It's written as an OpenCV module and looks well optimized.

However! I have no idea what their license on this is, so definitely check before you integrate it into anything you're going to distribute or sell.


Maybe some of the following steps might help:

Desaturate the image using any method your prefer (either a color register, luminance or lightness). It is much easier to perform the next steps if the image is only 1 dimensional
possibly apply some non-linear function over the pixels gray value. This might increase the visual color difference between light pixels and dark pixels. Additionally, it will reduce the color difference between dark pixels and very light pixels. you could consider it as a smooth threshold function.
As the roots have a light color and the ground is overall dark, we should look for 2 aspects: a pixel with a light value and to which the neighboring pixels have a dark value. To enhance this second aspect, you could add a Laplace-edge detect filter to the image (so the resulting value will be the current value and the value of the laplace filter).
At this point, if properly done, the roots should already somewhat pop out the image. One could follow this up with some Dilate - Erode filters (you can stack those) to even out background color differences.

Hope this helped a bit.


Guess the frequency response of a filter with little information

If you know that the system is linear and time invariant, the easiest method (assuming that you have no noise added in the process) is to let the system act on an impulse function. The Fourier transform of the output is the frequency response of the system.


How to match two signals that have similar shape but are scaled and expanded?

You can always match the extrema via a simple linear mapping, e.g. $y[n] = \frac{x[n]-min(x[n])}{range(x[n])}$ where $range(\cdot)$ is the absolute difference between the maximum and minimum values and $x$ is the input signal. This will get $y \in [0 .. 1]$ and then you can shift that interval to whatever (min,max) and baseline you like.

The dilations in time might be a bit harder to deal with.

If you are positive that the dilations are more or less by a constant factor, then prior to assessing the similarity between the two waveforms you can resample the longer one down to the length of the shorter one and then apply correlation. In fact, you could still base your similarity only on the comparable portions of the signals spectra, either by a direct error term or via something like coherence.

But if the dilations are not stable, then Dynamic Time Warping (DTW) (as was suggested) is the better option to go for as it will give you a "speed up / slow down" vector which you can use to either map one signal on top of the other or directly, as a similarity / disimilarity metric.

Hope this helps.


reconstructing time-stretched signals for filtering purposes

Keeping your “correct pitch” after a time stretch can become a more difficult problem than your original problem of filtering short events.  Both require a time-frequency resolution trade-off.


Why the frequency response plots (of causal filters) only show positive frequency?

This has absolutely nothing to do with causality. The frequency response of a real-valued filter (i.e., one with a real-valued impulse response) is (conjugate) symmetric, i.e., the negative frequencies are redundant. That's why it is sufficient to show the frequency response at non-negative frequencies only.

You can easily see that symmetry as follows. The frequency response of a system with impulse response $h(t)$ is given by

$$H(\omega)=\int_{-\infty}^{\infty}h(t)e^{-j\omega t}dt\tag{1}$$

If $h(t)$ is real-valued, taking the complex conjugate of $(1)$ gives

$$H^*(\omega)=\int_{-\infty}^{\infty}h^*(t)e^{j\omega t}dt=\int_{-\infty}^{\infty}h(t)e^{j\omega t}dt=H(-\omega)\tag{2}$$

So, for real-valued $h(t)$, we have

$$H(\omega)=H^*(-\omega)\tag{3}$$

With $H(\omega)=|H(\omega)|e^{j\phi(\omega)}$ we get

$$|H(\omega)|=|H(-\omega)|\tag{4}$$

i.e., the magnitude of the frequency response is an even function of frequency, and

$$\phi(\omega)=-\phi(-\omega)\tag{5}$$

i.e., the phase is an odd function of frequency.

The same is true for discrete-time systems.


Absolute and Square Summability & Integration of sinc function

Neither the integral

$$I_1=\int_{0}^{\infty}\left|\frac{\sin(x)}{x}\right|dx\tag{1}$$

nor the sum

$$S_1=\sum_{n=1}^{\infty}\left|\frac{\sin(n)}{n}\right|\tag{2}$$

converge. So the sinc function (sequence) is not in $L_1$ ($l_1$). A proof of the fact that $I_1$ diverges can be found here. For the sum $S_1$, note that since $|\sin(n)|\ge \sin^2(n)$ (for $n\in\mathbb{R}$) we have

$$S_1\ge \sum_{n=1}^{\infty}\frac{\sin^2(n)}{n}\tag{3}$$

Here it is shown that the series on the right-hand side of $(3)$ diverges, so $S_1$ must also diverge.

The integral

$$I_2=\int_{0}^{\infty}\left(\frac{\sin(x)}{x}\right)^2dx=\frac{\pi}{2}\tag{4}$$

can be computed in many different ways. Two possibilities are shown in the answers to this question.

Also note that

$$\int_{0}^{\infty}\left(\frac{\sin(x)}{x}\right)^2dx=\int_{0}^{\infty}\frac{\sin(x)}{x}dx\tag{5}$$

The computation of the sum

$$S_2=\sum_{n=1}^{\infty}\left(\frac{\sin(n)}{n}\right)^2=\frac{\pi-1}{2}\tag{6}$$

is discussed in this question and its answers. Note that a relation similar to $(5)$ holds:

$$\sum_{n=1}^{\infty}\left(\frac{\sin(n)}{n}\right)^2=\sum_{n=1}^{\infty}\frac{\sin(n)}{n}\tag{7}$$

So clearly, the sinc function (sequence) is in $L_2$ ($l_2$).


May a causal system depend on future values of its own output

By definition, the output of a causal system cannot depend on any future input or output values. However, I think there is a misunderstanding in your question. Suppose you are given the following difference equation:

$$a_0y[n+1]+a_1y[n]+a_2y[n-1]=b_0x[n]\tag{1}$$

The important point is that it is not possible to say if the system described by Eq. $(1)$ is causal or not, simply because Eq. $(1)$ does not uniquely describe one system, but in general there are three systems described by $(1)$. One of them is causal and its output can be computed from

$$a_0y[n+1]=b_0x[n]-a_1y[n]-a_2y[n-1]\tag{2}$$

Note that for any $n$ the output in $(2)$ only depends on past output (and input) values.

The second system described by Eq. $(1)$ is non-causal. The current output sample depends on past and on future output values:

$$a_1y[n]=b_0x[n]-a_0y[n+1]-a_2y[n-1]\tag{3}$$

And, finally, Eq. $(1)$ can also describe an anti-causal system, the output of which only depends on future output and input values:

$$a_2y[n-1]=b_0x[n]-a_0y[n+1]-a_1y[n]\tag{4}$$

The ambiguity of Eq. $(1)$ is also reflected in the $\mathcal{Z}$-transform domain, where $(1)$ corresponds to the following transfer function:

$$H(z)=\frac{b_0}{a_0z+a_1+a_2z^{-1}}\tag{5}$$

This transfer function, just like the difference equation $(1)$, does not uniquely define a system. Only if you specify the region of convergence (ROC) of $(5)$ will the system be uniquely characterized. Each of the three possible interpretations of Eq. $(1)$ discussed above correspond to a different ROC of $H(z)$.


Problem with Generalized Cross Correlation (PHAT filter)

Using the gcc-phat method, you can use the matlab function "gccphat(s2,s1,fs)" (R2016a) with s1 as the reference signal, s2 the other one and fs the sampling frequency.

Basically it will do a similar correlation calculation :

  Ncorr = 2*N-1;  % N size of the signal
  NFFT = 2^nextpow2(Ncorr);
  R12 = bsxfun(@times,fft(s2,NFFT),conj(fft(s1,NFFT)));
  r12_temp = fftshift(ifft(exp(1i*angle(R12))),1);
  r12 = r12_temp(NFFT/2+1-(N-1)/2:NFFT/2+1+(N-1)/2,:);

The only difference with your code is that the matlab function use the nearest power of two sequence length for the fft (as Peter.K suggest with "Nfft = length(x1) + length(x2) - 1") and doing an exponentiel stuff (and this part make the difference with a regular correlation, what you were doing)

But to return the time delay, the function is doing an other trick :

  lags = (-(Ncorr-1)/2:(Ncorr-1)/2).';
  lags = lags/fs;
  [~,idx] = max(abs(r12));
  tau = ts/(2*fs)+lags(idx); % ts as the total sample of the signal

I'm doing similar experimentation and this method is working great for long and short non-moving signal. the reverberation of the room may cause some issue, but if you can see an obvious delay when displaying the signal it should work fine.

Right now, i'm having difficulty to return a good time delay with a moving source (the function retourn 0 or +/- the length of the signal) And i'm still working on it... So if anybody got a lead, it will be much apreciate ;-)

Hope this will help you.

Peace

PS : Also, when doing an fft/ifft you really should use fftshift/ifftshift because the matlab fft/ifft is mixing some stuff during the process.


As far as I know, FSK-type modulations do not use quadrature. My answer below is about quadrature, not about a specific modulation.

There are two ways to look at quadrature (I and Q) modulation: baseband, and bandpass.

In baseband, the signal is complex. You have two independent signals, each modulated using PAM, $s_I(t)$ and $s_Q(t)$. The quadrature signal is $$s_{BB}(t)=s_I(t)+js_Q(t).$$ Naturally, this signal can't physically exist. To transmit it, you need to upconvert $s_{BB}(t)$ to some carrier frequency $f_c(t)$. It becomes a bandpass signal: $$\begin{align} s(t)&=\Re(s_{BB}(t)e^{j2\pi f_ct})\\ &=s_I(t)\cos(2\pi f_ct)-s_Q(t)\sin(2\pi f_ct). \end{align}$$

Both signals carry the same information and are equivalent in many ways. $s_{BB}(t)$ is useful to simulate and analyze the system, because it does not have a carrier frequency. $s(t)$ is the signal that is actually, physically transmitted.


Can you explain how to imagine frequency in 2D? I need to understand how the concept of frequency exists in images

Frequency, by its definition, means how frequently something occurs. You can treat images as planes where you have data points, ie. the pixels of the data. In images a frequency can be thought of as how frequently the sinusoidal components of the structure repeat per unit distance (where your unit would be the distance of a pixel).


For an $N$-periodic sequence $y[n]$, the sum over one period is always the same, no matter how you choose the first index:

$$\sum_{n=0}^{N-1}y[n]=\sum_{n=\Delta}^{N-1+\Delta}y[n]$$

for any integer $\Delta$. Try to verify that the sequence you sum over is indeed $N$-periodic.


Speech recognition using MFCC and DTW(Dynamic Time Warping)?

At the very beginning let me warn you that DTW approach is suitable only for spoken word recognition. Nonetheless it is interesting as a basic exercise.

I assume that you have a database of training files (templates) and you already extracted MFCC's for those. Each training file contains one utterance for a given word class, i.e. "Hello", "Bye". I also imagine that each class has more than one example (plausible).

Here is what you should do during the recognition stage:

Once you have the recording for the whole word, you calculate the MFCC's for each frame. You want the recording of the whole word, not just part of it.
Now you calculate the distance between your recording and each of the templates in your database. In case of DTW you will be calculating the cost between each of 13 dimensional frames (simple distance metric/norm, i.e. Euclidean, Manhattan, etc.). Once the DTW algorithm is finished, you will end up with the distance value (upper right matrix entry as can be seen below) between your test sample and each of the templates.
The last step is to make a decision: to which class your test sample actually belongs to? You can either do it by picking the class of template with the minimum DTW distance. But what is even better, you can use the k-Nearest Neighbours for that.


Using 8 bit samples for speech comparison

Will that affect the quality of speech comparison, and by how much?

That is impossible to tell without knowing what the Node.js thing does internally; I think it's a bit much too ask for us to search for what you meant. As a comment: signal processing in JavaScript sounds like a bad idea, performance-wise and development-wise; it's really not what JS was designed and optimized for, and there's a significant lack of libraries, let alone efficient ones. I always urge people to use the right tools for their job, and JS isn't that here, I think. </comment>

That being said:

From a pure signal point of view, you can model the effect of quantization as noise. The problem here being that this kind of noise is neither uncorrelated to the signal, nor necessarily white. To make things a bit harder  on the compensation side, its amplitude typically isn't Gaussian even. Oh, well, but here goes quantization noise power, a figure that is both very important to understand the maximum SNR you can get out of a digital system and doesn't say much as long as you don't say how well the signal processing can deal with this specific kind of noise.

Luckily, for "small enough" quantization steps, a little stochastic consideration¹ says that the assumption that quantization noise (QN henceforth) is additive is pretty justified.

Now, assuming your signal amplitude is really uniformly distributed, and your ADC being perfectly uniform, the amplitude $\text{SNR}_Q$ (signal-to-QN-ratio) for an M-bit ADC becomes

$$\begin{align*}
\text{SNR}_Q &= 2^M\\
\text{SNR}_Q \text{[dB]}&= 20 \log_{10}(2^M)\\
 &= 20  \log_{10}(2)\,M \\
&\approx 6 M\text,
\end{align*}
$$
which implies that for 8bit, your $\text{SNR}_{Q,8b} \approx 48\,\text{dB}$, and for 16bit $\text{SNR}_{Q,16b} \approx 96\,\text{dB}$.

Now, speech definitely isn't uniform in amplitude; it's a bit hard to justify this model without knowing what your recording looks like, but I'd rather say it's composed of sines; in that case, you get an additional $1.8\text{dB}$ noise for both cases.

Point here is that I doubt that a "real world, non-studio equipment, non-anechoic-chamber silence" speech recording will ever be any close to $48\,\text{dB}$, so probably, no, that's a fine choice, and by the way you ask this question, I kind of doubt it will be trivial to extend the algorithm (which probably uses a lot of elegant numerical math internally) to 16bit, anyway.

¹ the amount of noise power correlated to the signal necessarily being very limited, and an mutually independent set of sufficiently many i.i.d. realizations (Normal due to CLT) added up still being Normal if their moments weren't too different ...


HINT:

Prove that
$$y(t)=\frac{1}{\pi}\int_{-\infty}^{+\infty}\frac{x(\tau)}{t-\tau}d\tau=\frac{1}{\pi}\int_{-\infty}^{+\infty}\frac{x(t-\tau)}{\tau}d\tau\tag{1}$$
Use the right-most expression in $(1)$ to show that the system is time-invariant.

Equivalently, you can directly show that the original input-output relation is a convolution integral, from which it follows that the system must be LTI. Then find the expression for the impulse response $h(t)$ and show that $h(t)\neq 0$ for $t<0$, which means that the system is not causal.

EDIT:

The same variable substitution that can be used to show $(1)$ can also be used to directly show the equivalence of the response to $x(t-t_0)$ and the delayed response to $x(t)$:

$$y_2(t)=\frac{1}{\pi}\int_{-\infty}^{+\infty}\frac{x(\tau-t_o)}{t-\tau}d\tau$$

$$\begin{align}y(t-t_o)&=\frac{1}{\pi}\int_{-\infty}^{+\infty}\frac{x(\tau)}{(t-t_o)-\tau}d\tau{\Huge|}_{\tau=u-t_0}\\&=\frac{1}{\pi}\int_{-\infty}^{+\infty}\frac{x(u-t_0)}{(t-t_o)-(u-t_0)}du\\&=\frac{1}{\pi}\int_{-\infty}^{+\infty}\frac{x(u-t_0)}{t-u}du\\&=y_2(t)\end{align}$$


Biquad coefficients using Magnitude (or Phase) Invariance Mapping Method

2nd order MIM LPF results looks great for selected cutoff frequency (fs = 44.1kHz) as seen in plot (frequency range upto 23873Hz):

Coefficients for plotted 2nd order filter:

a [0.437417389844334 0.1220119196463959 -0.04071587087954469]
b [1 -0.6951293456643256 0.2138427842755107]

EDIT:

It looks like there are some problems in magnitude response when fs is increased and/or when fc is reduced. MIM/PIM paper mentions this as an issue.

So, it looks like MIM can be used for to improve 2nd order LPF magnitude response especially when using low fs and high enough fc. Magnitude issues can be bypassed by using higher order filters. Fortunately MIM method works correctly for my purposes.

EDIT: First and the second D(s) given in paper examples are wrong. The correct ones are (5s + 2) and (s + 10) / (s^2 + 4.1s + 0.4). With those all the MIM calculations become right. The third example (PIM) was correct.


How to perform a Stockwell Transform on a 2D matrix

There is a MATLAB code available for the Stockwell transform here. It also includes a code for the one-dimensional transform but it is less efficient that the Mathworks file exchange code that is mentioned in the question. This is because one uses for loops and the other doesn't.

Using this 2D S-transform function, you can quickly run into memory issues. You can easily modify the code to specify the exact positions of x and y that you need, so that S is not 4 dimensional.


They will work as you show but you can filter and decimate to 2 samples per symbol followed by resampler and TED if your carrier offset is reliably less than $S/4$ where S is the symbol rate (it will work with offset beyond that but with possibly unacceptable degradation) and after than matched filter at 2 samples per symbol with and you can use the pre-decision 1 sample per symbol samples for carrier recovery. The performance of the TED is better prior to matched filtering typically (the zero crossing jitter is less) unless you do additional prefiltering for timing recovery. The TED works well with large carrier offsets hence you need not do carrier recover first. Another option if your carrier loop bandwidth allows the additional loop delay is to do the complex rotator to correct for carrier offset before the interpolator while still doing the phase measurement on the final pre-decision samples so the TED operates with 0 carrier offset, but I would only consider this when dealing with much larger carrier offsets.

(Also I think your block labeled “DDS” is really an “NCO”; a DDS is the cascade of an NCO with and D/A converter)


Linear phase FIR filters should have some sort of symmetry (even or odd)  in their coefficients which will reduce the number of independent multiplications by about two yielding an efficient realization compared to a nonsymmetric FIR filter of the same length.

However, linear phase FIR filters are important not because of their computational efficiency but for their linear phase property.


White noise vs. delta pulse and Ultraviolet catastrofe

As Jason was saying in the comments, the power spectral density of white noise is flat.  This is equivalent to saying that the autocorrelation of white noise is a delta dirac function (i.e. that there is no correlation, positive or negative, between one noise sample and another), not that the noise itself is a delta dirac function.


Why ideal filter cannot realised practically?

Most typical examples of the ideal filters are the classical brickwall frequency selective filters (lowpass, highpass, bandpass etc) which are defined as ideal because of their impossible to realize ideally selective frequency responses.

Their frequency responses include exactly flat passbands, exactly flat stopbands having absolute zero gain, and zero width transition bands. These features make it impossible to realize them using either existing physical devices or finite time computable algorithms...

Assuming an LTI system model for such filters, then their impulse responses tend to be of infinite length and typically be acausal; another reason for denying their physical realization. Yet causality alone is only a concern for real-time processing of time domain signals. It does not prevent realization of such filters in offline procesing where the whole signal data with its past and future are already available. Yet even in such an offline processing mode the fact that they require an infinite length of data is the reason for preventing their realization.

Another point of view into their unrealizability is that input/output relationship of such ideal filters cannot be described by finite order differential / difference equations. In principle any filter whose input/output relation can not be described by a finite order differential / difference equation is said to be unrealizable; i.e., there is no finite time computable algorithm to implement it...


Question related to Discrete signal convolution

Your solution based on Example 2.2 of A.Oppenheim's Signals & Systems is wrong.

In that example, the output
$$y[n] = \sum_{k=-\infty}^{\infty} {x[k]h[n-k]}$$

is computed by the graphical method (a.k.a. flip-and-drag method).

Based on the nonzero range of input $x[n]$, this becomes:
$$y[n] = \sum_{k=-2}^{2} { x[k]h[n-k] }$$

Also based on the nonzero range of $h[n]$, it can be seen that $y[n]$ has nonzero samples in the range from $n=-5$ to $n=2$.

In the flip & drag method, it's best to draw the shifted functions $h[n-k]$, to guide the computation process.

The key is to recognise the nonzero range for the flipped & dragged function $h[n-k]$ ; time reversed $h[-k]$ shifted by "$n$", drawn on an axis of $k$ and thus considered as a function of $k$.

You can see that the overlap between the functions $x[k]$ and $h[n-k]$ begins at the shift of $n = - 5$, from the left side of $x[k]$ (the sample $x[-2]$) and continues until the shift for $n = 2$ at the right side of $x[n]$ (sample $x[2]$). For any other "$n$" value, there is no overlap between $x[k]$ and shifted $h[n-k]$ and output $y[n]$ becomes zero.


How to account for sample rate in a FIR Digital Filter?

The sampling frequency is hidden in your variable $\theta$. The angle $\theta$ determines the frequency of the notch relative to the sampling frequency, i.e.

$$\theta=2\pi\frac{f}{f_s}\tag{1}$$

where $f$ is the notch frequency in Hertz, and $f_s$ is the sampling frequency (in Hertz). So if you want to change the sampling frequency while leaving the notch frequency unchanged, you need to change the angle $\theta$ according to Eq. (1).

The gain of your filter depends on the normalized frequency variable $\theta$. At DC the gain is $2(1-\cos\theta)$, and at Nyquist you get a gain of $2(1+\cos\theta)$. By normalization you can keep one of those two gain factors normalized to $1$, but never both. If you want a gain of $1$ at DC and at Nyquist you need an IIR filter or a higher order FIR filter.


How large FFTs can pull signals out of the noise floor?

If you look at the formula of a single DFT bin

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

you'll notice that his is essentially a correlation of $x$ with the complex sinusoid $e^{-j2\pi k\frac nN}$.

That means the DFT can just be understood as a filter bank of matched filters for single tones that fall in the DFT "raster".

Hence, you simply get FFT length-based processing gain: The length of the sum.

But: you probably don't have perfect knowledge of the exact frequency of the signal you're trying to detect¹! So, you can't put things into that perfect DFT raster.

Now, the larger you choose the FFT length $N$, the finer that raster will get, but also, the longer your observation has to be, and the more compute power you'll need.

At some point, the DFT stops being the best possible tone detector, and superresolution techniques become relevant. In this case (weak tone, you're sure that you've only got exactly one tone in your signal), the ESPRIT algorithm with a long observation period leading to the autocovariance matrix estimate that it takes as input, would probably work very nicely.

¹ There's inevitably frequency error in your receiver, and in your transmitter. Papers that start with We assume perfect synchronization typically skip the hard part of making a system work...


What DSP method are consumer products like the game Rocksmith and guitar polyphonic tuners likely to use for multiple f0 detection?

Check out chapter 1.3 of this IRCAM paper on multi-F0 estimation. It discusses the difficulties in extracting multiple F0s from a recording, including the handling of overlapping partials, transients, and reverberation, as well as the modeling of domain-specific sources with varied spectral properties.


Comparing two time domain signals in a scale and shift-invariant way

First remove the DC offset (C2) from each signal by subtracting the average level from each sample, then removing the arbitrary scaling (C1) of each signal by dividing each sample by the average sample magnitude.


Strange results from my FM synthesis algorithm

It should be more like this:

modulator = A_mod * sin(2 * pi * fm_freq * i / sample_rate)

carrier = sin (2 * pi * (freq * (1 + modulator)) * i / sample_rate)

The two main changes are:

the modulation signal needs a scale factor (A_mod) which determines the amplitude of the FM - typically A_mod << 1
the carrier frequency needs to be multiplied by (1 + modulator) rather than just modulator, since it needs to be centered about freq with a deviation of +/- A_mod * freq (in your version the centre frequency is zero and your instantaneous frequency is varying between +/- freq !).


First, motivations for an engineer solution, which might be the initial reason of such a question:

The compression ratio of a JPEG encoded image result from a sort of average from different image blocks "compressibility". The JPEG principles rely on two combined aspects:

the average luminosity from block to block varies slowly, hence it can be nicely predicted from preceding blocks. This is the differential prediction used on DC coefficients (here your 240);
local stationarized parts of images can be faithfully approximated by a few quantized 2D cosines.

The scheme you provide is a  simplified version, here are a few. First, the DCT is applied on $8\times8$ blocks, not $4\times4$. Hence, it is more difficult to generate long runs of zeros after quantization. Second, the quantization scheme is linear (dividing by 2). Third, with only one block, you cannot gain from inter-block prediction on DC. Fourth, with  only 15 AC coefficients, the (value,run) Huffman coding might be a little inefficient. To add a few, non-adaptive tables, lack of header, may highly impact the estimated compression ratio  on such a small file.

Last, the "DCT coefficient uses 10 bits" is a bit awkward. So, I interpret it as "The DC coefficient uses 10 bits". With that, you are left with 4 values and their occurrences:

$2$ (2 times)
$1$ (3 times)
$-1$ (5 times)
$0$ (5 times)

from with you can build a Huffman tree, and a non-prefix code-word for each of them. Adding their overall length to the above 10 bits can give you a rough estimate of the total number of bytes actually needed. You can use for instance an  Online Huffman coding tool.

Second, you can work out your own JPEG format, on $4\times4$ blocks, and input an sufficiency large image with those repeated $4\times4$  blocks, to compute actual compression ratios.


space frequency representation of images in matlab

Here is pseudo matlab code for doing it.

X = 2D_image;
Y = 3D_array_for_storage;
Z = 4D_array_for_final_result;
for m = 1 : Num_cols_in_X
  Y(m,:,:) = spectrogram(X(:,m));
end
for m = 1 : size(Y,2)
  for n = 1 : size(Y,3)
    Z(:,:,m,n) = spectrogram(Y(:,m,n));
  end
end

The result is a 4-D array that is row-time x row-freq x col-time x col-freq. (It does not have to be time and frequency per se. The units of the axes are arbitrary; they could be second, meters, volts, etc.) As for how you display, I have no idea what's best for your case. A movie could display any three of the dimensions. Perhaps you could tile the fourth dimension? It would all depend of the specifics.


JPEG: In image compression why are higher frequency components divided by a larger number during quantization

Note that the matrix which divides the DCT coefficients, performs quantization and the larger the divisor the coarser will be the step size at that frequency.

And that's the very essence of the perceptive image coding, where less perceived (roughly less important) content is coarsely quantized while the more important and dominant low frequency content is finely quantized and represented better. This will result in more quantized zeros at high frequencies which eventually yield more compression while maintaining less distortion compared to a uniform quantization.

Hence this takes advantage of human perception characteristics. Had you thrown away the low frequency information most people would object immediately to the result, while they would not even recognize the loss of high frequency information until some aggressive threshold.


Sequence expansion by zeros and interpolation - does it insert additional frequencies?

The key to understanding what inserting zeros does is to understand two things: what samples represent in the time domain (because we want to insert zeros in the time domain), and what they represent in the frequency domain (because we want to know what it did to the spectrum).
First, sampling is a type of modulation (PAM—Pulse Amplitude Modulation), equivalent to multiplying an impulse train by our analog signal, and creates images in the frequency domain. When converted to digital sample values, we have PCM (Pulse Code Modulation), a common term for digital audio. Here is an example spectrum represented by our samples; the original analog signal's spectrum is shown in green, and images in red:

These images are the price we pay to represent the analog signal as samples. As such, the usable bandwidth is from 0 Hz up to (but not including) half the sample rate. Above that is a backwards image of our original spectrum, with images repeated around multiples of the sample rate. It's OK, we remove the images when we convert back to analog, using a DAC's lowpass filter.
In the time domain, samples represent impulses. They are instantaneous values taken at a constant interval. Inserting zeros changes nothing except what we consider the sample rate.
For instance, sample a signal one time per second. That represents an impulse train, so consider playing it back raw, as an impulse train.
Now consider placing a zero-valued sample between each of the original samples. Consider playing it back as an impulse train but at twice the original rate, two times per second.

Can you see that nothing has changed in the signal, except the sample rate? Likewise, if we look at the spectrum, nothing has changed in the frequency domain—this is obvious, since the time domain signal has not changed.

However, our usable bandwidth has doubled. The first inverted image now lies in our usable band, now shown in green. It will not be removed by the DAC when played back through a DAC at the new, higher rate. And it will be a problem with any non-linear processing in the digital domain.
That is why we follow the zero-insertion with a lowpass filter (or combine the two steps for efficiency). Here it is again after proper filtering, below half the original sample rate:

So, the answer is that nothing changes in proper integer sample rate conversion by zero-insertion, and the result is as perfect as the lowpass filter used to clean up the exposed images.


Which parameters should I take care of while converting .au to .wav for Genre Classification

From a raw "music" point of view, both .au and .wav file formats may introduce a form of lossy compression, as compared to the basic Pulse Code Modulation (PCM). Partly because some parties to drift away for the standard standard, including mu-law companded versions, there are some inconsistencies in the AU and WAV formats.

As the .au is lossy compressed, what is important is that you do not lose even more information while converting to .wav. From your discussion with @Olli Niemitalo, the information seems preserved. However, I am a little worried that the file sizes are the same (.au file) (I tend to think the .au ought to be smaller, with the basic settings):

  The Au file format [...] Originally [it] was headerless, being simply 8-bit µ-law-encoded
  data at an 8000 Hz sample rate. Hardware from other vendors often used
  sample rates as high as 8192 Hz, often integer factors of video clock
  signals. Newer files have a header that consists of six unsigned
  32-bit words, an optional information chunk and then the data (in big
  endian format).

So you can check the header information to be sure of how the .au was compressed, what channels are here, etc. As long as this information remains in the .wav file, the chances of affecting genre classification are very mild.


Kalman Filter : How measurement noise covariance matrix and process noise helps in working of kalman filter , explain intuitively please?

Considering a linear dynamic system framework, if there were no noises in the process equation and measurement equation, then one could use deterministic means to compute the state exactly, without any need for an estimator.

However, the unavoidable presence of noise in the state transition and the measurements make it impossible to compute the state exactly but only approximately. As the measurements are not true but in error.

Kalman filter (similar to other statistical estimators) improve the state computation by referring to the statistical character of the noise present in the measurements and state transtitions. Kalman filter also benefits from the dynamic equation of the system in predicting the state update.

So those noise covariance matrices are used to specify the Kalman gain K, which tries to correct the error in the deterministic computation of the state update.


Have you seen this MusicDSP Phaser code and discussion by Ross Bencina? It is commented, and says that its implemented to be "clear, not efficient". Cheers.


How to find the DTFT of unit step with $(-1)^n$ scalar multiplied?

First note that $e^{j\pi n}=(-1)^n$. Furthermore, note that $u[n]-u[n-N]$ is a rectangular impulse of width $N$.

Now, you can calculate
$$
\textrm{DTFT}\left\{u[n]-u[n-N])\right\} = \sum_{n=0}^{N-1}e^{-j2\pi nf}=\frac{1-e^{-j2\pi Nf}}{1-e^{-j2\pi f}}
$$
which is the rule for the geometric sum.
Furthermore, you will find

$$
\textrm{DTFT}\left\{\exp{(j\pi n)}\right\}=\delta\left(f-\frac{1}{2}\right)
$$

Finally, the result follows from the convolution theorem, i.e. multiplication in time becomes convolution in frequency:

$$
\textrm{DTFT}\left\{\exp({j\pi n})(u[n]-u[n-N])\right\}=\delta\left(f-\frac{1}{2}\right)\star\frac{1-e^{-j2\pi Nf}}{1-e^{-j2\pi f}}=\frac{1-e^{-j2\pi N(f-\frac{1}{2})}}{1-e^{-j2\pi (f-\frac{1}{2})}}
$$


It would be better to show that property using pencil and paper, because you're fooling yourself with that Matlab code. Of course we have
$$\delta[n-n_0]\star g[n]=g[n-n_0]\tag{1}$$
Now define $h[n]=g[n-n_1]$. According to $(1)$ we have
$$\delta[n-n_0]\star h[n]=h[n-n_0]\tag{2}$$
And, consequently, with $h[n]=g[n-n_1]$,
$$\delta[n-n_0]\star g[n-n_1]=g[n-n_0-n_1]\tag{3}$$
No need to waste computing power to show this.


How to compute dBFS from rms-values if the samples exceed 32bits

There are a few different ways:

Use a fractional multiply or data types. Unfortunately that's not a standard "C" or "C++" instruction but most embedded processors do support it and have some sort of non-standard way to call it.
Use a "long int" 64-bit integer data type, do a double precision multiply.
Convert everything to floating point (float or double) before multiplying.
Scale the numbers by $2^{-16}$ (arithmetic right shift by 16-bits) before squaring.

Option 1 is typically the "best" if your processor supports it as a native instruction. Option 4 is often "good enough" if you have half-way decent signal to noise ratio but that may depend on the specifics or you application and data.
Note on "fractional" data types: This is a fixed point type that is scaled so that the data range is assumed to be  from -1 to +1, hence a multiply can never overflow. That's actually the most "natural" type for fixed point signal processing. Multiplying two 32-bit numbers yields a 64 bit result. The difference between an integer and fractional multiply is that the integer multiply keeps the lower 32 bit and the fractional keeps the higher 32 bit of the result.
Both integer and fractional format are subsets of the "Q" number format. See https://en.wikipedia.org/wiki/Q_(number_format).  Loosely speaking integers are Q0 and fractionals are Q31(signed).


Designing complex magnitude compensating FIR filter

Here is an example I already had in Python demonstrating the Wiener-Hopf equations detailed here Compensating Loudspeaker frequency response in an audio signal to solve for the channel equalizer. This application by the OP is ideal for this given the static channel condition as long as deep spectral nulls don't exist, over a iterative LMS or RLS equalizer solution such as those detailed here (which are great for dynamic channels that don't exhibit deep frequency nulls in band): Recursive Least Square Adaptive Linear Equalizer
For channels with deep spectral nulls in band (frequency selective fading), a decision-feedback equalizer would be the better approach due to degradations from noise enhancement with the linear equalizer approaches.
Eye diagram of undistorted waveform with no channel distortion, as it would be received after the 2nd RRC filter in the receiver:

Tx Spectrum of undistorted waveform prior to channel distortion:

Example complex channel distortion that we will equalize:
cir = [.019*np.exp(1j*.2), .021*np.exp(-1j*.24), .08*np.exp(1j*2.4),  .032*np.exp(-1j*.14), -.032*np.exp(-1j*.14)]

This is the same eye diagram and Rx spectrum above after passing through our channel (yuck, what a mess! How do we possibly demodulate this????). Notice the assymetrical passband, a tell-tale sign that the compensating filter MUST be complex.

Compute Equalizer Solution Using Tx and Rx Waveforms
(Note this is Python 3.8, but nothing jumps out at me that wouldn't be compatible with earlier versions)
import scipy.linalg as linalg
import numpy as np

def convmtx(h,n):
    # creates the convolution (Toeplitz) matrix, which transforms convolution to matrix multiplication
    # h is input array,
    # n is length of array to convolve
    return linalg.toeplitz(np.hstack([h, np.zeros(n-1)]), np.hstack([h[0], np.zeros(n-1)]))

# Compute equalizer given known tx and rx waveforms

omit = 100    #initial samples to exclude
shift = 0  # number of samples to shift dominant equalizer tap to the right

ntaps = 20
depth = ntaps * 30
delay = ntaps//2
A = convmtx(rx[omit+shift: omit+shift+depth], ntaps)
R = np.dot(np.conj(A).T, A)
X = np.concatenate([np.zeros(delay), tx[omit: omit+depth], np.zeros(np.int(np.ceil(ntaps/2)-1))])
ro = np.dot(np.conj(A).T, X)

equal_coeff = np.dot(linalg.inv(R),ro)

Below is the resulting frequency response for the equalizer and the equalized waveform. The details for the parameters used such as equalizer number of taps and depth are further detailed in the earlier links given. Note how the equalizer can only determine the compensation where signal energy exists (otherwise it has no way to measure the channel), therefore it is important to use a test (sounding)  waveform that fully covers the bandwidth of the channel of interest. Pseudo-random sources or any white noise source is ideal given the broad spectral coverage.


You could try low-pass filtering the input signal to get smoother zero-crossings (or even band-pass filtering if you have a good idea of the frequency location of the sine wave). The risk is that if sample-accurate phase information is essential to your application, the additional lag from the filter might be a problem.

Another approach: instead of trying to transform the sine wave into a square wave, what about getting an independent square wave oscillator to align itself in phase/frequency with the sine wave? This can be done with a phase-locked loop.


Discrete-time Fourier transform of $\frac{\cos(\frac{n\pi} 6)}{(n+3)\pi}$


Reconstruction of a signal using 1D discrete wavelet

This sounds like a typical effect due to downsampling after a convolution with a non perfect lowpass, when one of the frequency in the approximation becomes close to the cutoff $f_c$ at some scale, of the form $f_c = f_s/2^j$ for some $j$. Notice that the sum of $120+129.9 \approx 2*125$. Set the second frequency to $124$, you can see a second peak at $126$ Hz, similarly.

So this is aliasing, resulting from a two-fold subsampling after a FIR lowpass, whose transition band is not sharp enough. This can be avoided be sharper wavelets (like $M$-band wavelets), or resorting to shift-invariant or stationary wavelets (swt.m instead of dwt.m). The following real-life example come from seismic data processing: the data is corrupted by a $60$ Hz powerline sine (US standard), and acquired at $4$ ms. The top graph is a seismic data spectrum, after the data has been processed by the discrete wavelet, and shows a second peak at $65$ Hz ($60+65=125$), while the second one is processed by a SI wavelet.

Your question revived my memories from 2002, Efficient Coherent Noise Filtering, An application of shift-invariant wavelet denoising (source of the picture).


A vocoder is a type of synthesizer originally developed for speech (or the human voice, thus the "vo" prefix), even before the use of digital processors.  The Bell Labs Voder was demonstrated in 1939, and VoCoders were reportedly used for encrypted voice communication during WW2 (according to Wikipedia).

DFTs can be used to analyze which audio spectrum peaks are to be used to "code" the speech, and short STFT windows are used to get a fine enough time resolution to isolate typical phonic atoms over reasonably stationary time periods.

But the frequency bin centers of a short FFT are not spaced with sufficient frequency resolution to capture pitch, overtones, and modulations thereof, accurately enough to sound human.

But, by taking into account the FFT phase information, the difference in the phase of FFT result bins between adjacent offset windows can be used to adjust (or interpolate) the apparent or actual frequency of the voice codec synthesizers to more accurately model narrow-band frequency peaks than by just using the magnitude alone of the STFT FFT results.  Thus the "phase" prefix to the analysis-resynthesis vocoder process.

Note that if you slide two Windows apart, the amount of phase change of stationary spectrum should be linear with frequency.  If analysis shows a little more phase change, the frequency is higher, a little less phase change, lower.  And this phase difference will be measurable even for differences in frequency much less than the Fs/N bin spacing of each FFT.  Thus there is information gain in using succesive frames and phase, rather than single frames and magnitude only, for the frequency estimation of narrow-band spectral peaks.  Phase can either be stretched or compressed across frames during FFT resynthesis, or the frequency of higher resolution oscillators can be adjusted to match the estimated frequency during synthesis.

With increases in DSP performance, phase vocoders are now used for music as well as speech, but the "vo" prefix remains.


A few things:

Not all instruments have partials neatly defined by the $h_n = nf$. Kettle drums would be an obvious example, but bells are another. In this case, you will be way off looking in the "usual places" +/- 5%.
Many instruments have only even or odd harmonics. Others are designed to omit some harmonics. The piano, for example, is designed to have almost no 7th harmonic. This may not matter to you, but it's worth pointing out in case it does.
It is common for higher harmonics to be "out of tune", that is, not in the "usual places". An engineer who once worked as a piano tuner once admonished me for using harmonics to tune my upright bass (a technique that works well for lighter-stringed instruments), because, he said, the partials are more out of tune with heavier stringed instruments. He said it's especially true told of higher harmonics, explaining that the math that gives you the simple formula for partials assumes that strings have no mass. 5% strikes me as a little much (almost a semi-tone!), but I don't think it's out of the question.
Many instruments have vibrato. This throws all your harmonics out of wack:


Kiss fft does not work after giving it more than 32 samples

It is possible that the program expects a power of 2 samples. Does it work with 64 samples? The FFT by definition usually refers to power of two series.


Avoid complicated calculations, let $X$ and $Y$ be i.i.d. standard normal random variables, your random variable $Z$ has the same distribution of $V$
$$V \triangleq \left(\frac{X}{\sqrt{X^2+Y^2}},\frac{Y}{\sqrt{X^2+Y^2}} \right)$$
(easy to see $\lVert V\rVert=1$ and the angle of $V$ is equivalent to the angle of a circularly symmetric Normal hence uniform).

This kind of $V$ is one of the constructions of a point uniformly distributed on circle (which can be generalized to $(n-1)$-sphere, see Sphere Point Picking and for example this answer).

Thus the PDF of $Z$ is simply the reciprocal of the unit circle's circumference.
For $Z_\rho = \rho e^{j\Theta}$ with fixed $\rho$ and uniform $\Theta$,

in polar coordinates (where infinitesimal area is $r dr d\theta$),
$$f_{R,\Theta}(r,\theta)=\frac{1}{2\pi} \delta(r - \rho)$$


I've done multi-resolution spectrograms before, so they're not impossible (a unicorn).  One method is to zero-pad the shorter data window (for better time locality resolution) to the same length FFT as the longer data window (for greater frequency resolution).  Then you can simply scale for window length, and cross-fade the two sets of FFT magnitude results over some frequency range.

Another possibility (instead of zero-padding) is to use interpolation of the lower frequency FFT or lower time resolution FFT results to create more possibilities for a higher density of spectrogram plot points.

One academic paper suggested computing the local vertical and horizontal contrasts at each spectrogram point and dynamically picking the FFT set with the highest contrast (over the vertical for narrow frequency spectra or across the horizontal for short time events) to plot at each point.

If you keep both (or more!) sets of (overlapped) STFTs, your transform isn't lossy, but redundant, and the resulting image can use that redundancy to obscure less information, usually due to graphical blurring/plotting with typical linear interpolation and/or non-linear color mapping.

Even for reconstruction you can use both sets of data, for instance feeding the short IFFTs to a tweeter and the longer IFFTs to a sub-woofer, with a cross-over in the middle (be careful of phase).


Measuring the open loop transfer function in closed loop: what's the better approach?

You have two signals to use, the reference and the measurements. You generate the error from these two signals. There is no additional information in that signal.

Closed-loop identification is more difficult than open-loop because the feedback makes the input and the noise correlated. Consider the standard feedback system

where $\bar{C} = C H$, $G_{0} = G$, $r_{2} = 0$, $r_{1} = C r$ to make it equivalent to your block diagram. A direct identification of a model $G_{m}$ means minimizing
$$
V(\theta) = \frac{1}{2\pi} \int_{-\pi}^{\pi} \left| G_{0}(\text{e}^{j\omega}) - B(\text{e}^{j\omega}) - G_{m}(\text{e}^{j\omega},\theta) \right|^{2} \frac{\Phi_{u}(\omega)}{\left|H^{*}(\text{e}^{j\omega})\right|^{2}}
$$
where the bias term is given by
$$
B = \frac{\bar{C} S_{0} H_{0} \Delta H {\sigma_{\bar{e}}}^{2}}{\Phi_{u}}
$$
and $\Delta H = H_{0} - H^{*}$, where $H^{*}$ is an assumed, fixed noise model, and $S_{0}$ is the input sensitivity function of the diagram. The bias in the model estimate $G_{m}$ will be small if

The noise model $H^{*}$ is accurate.
The signal-to-noise ratio at the input $u$ is large.

The estimation can be improved by using more elaborate identification methods. A great survey is found in Chapter 9 here.

The choice of excitation signal in order to minimize your estimation variance is discussed in Chapter 7 in the above reference, and another great paper can be found here.


How to determine windows size of moving average filter

Well, you'd look at the spectrum of your noise, and deduct an appropriate filter frequency response from that. When doing so you'll realize that the moving average is not a good filter (having sinc-shaped frequency response), and that you should be using a better filter design, aligned to what your signal actually needs.
So, you'd do three things:

you come up with a physical model of how fast the acceleration changes
you verify that with a look at the spectrum of your recorded signal, and look for where the noise in your recording is
you design a filter that has a passband where your expected signal is, and especially strong stopband suppression where only noise happens (hint: many accelerometers are swinging micromachined things, and there might be significant noise from that at subharmonics of the swing frequency)

Also, this is the classical example for which Kalman filters (or nonlinear variants of it, depending on how you plan to model your state) are used, instead of plain double integration.


How to get some law parameters from its FFT figure?

If you feed an FFT a finite set of samples from a function of theoretically infinite range, you have implicitly windowed your function with a rectangular window before the FFT.  Thus, the FFT result will be the convolution of the spectrum of your function and the transform of a rectangular window (a periodic Sinc).  This convolution will be an identity operation only if the samples happen to span exactly an integer multiple of periods of a purely periodic function.  So how much your FFT result will change or differ will depend on the offset and length of your FFT aperture/window with respect to some reference point on and the scale of your function.

In the case of a single sinusoid in zero noise, I reference 2 papers by Clay Turner on solving for the 3 unknowns using only 3 or 4 non-aliased samples, somewhere on my web page: http://www.nicholson.com/rhn/dsp.html   No FFT needed.


Beam width of linear array using MVDR beamformer

I think you're mixing time domain and frequency domain, your signal $x$ is in time and you're multiplying it by a complex exponential, generating a complex signal that is not the Fourier transform of $x$, if you want to delay $x$ in order to simulate the angle of arrival of the signal you have to do $x_{m}(n)=x(n-\tau_{m})$ (in the time domain) or $X_{m}(\omega)=e^{-j \tau_{m}}X(\omega) $ being "m" the microphone number.

Also you have to include the wave number in the steering vector $k=\dfrac{2\pi}{\lambda}$.


How the Fourier transform of a cosine signal is existed?

The Fourier Transform is defined only for functions that are absolutely integrable or in other words $f\in\mathcal{L}_1$, where $\mathcal{L}_1$ is the set of all absolutely integrable functions.

If you want to be mathematically rigorous, then you should assume that the Fourier transform of the cosine does not exist, as it is not integrable. Apart from it, the Dirac's delta does not fulfill the axioms of a function in modern mathematics, so you cannot say that a Fourier Transform of a function is something that is not a function. Instead of that, the Dirac's delta is defined as a measure that can be used in the Lebesgue integral, that satisfies certain axioms(https://en.wikipedia.org/wiki/Dirac_delta_function). However, the use of the cosine is valid for the common modulation operation, and is mathematically consistent (meaning that its Fourier transform exist):

$\int_{-\infty}^{\infty} |x(t)cos(2\pi t)e^{-i2\pi ft}|dt\leq \int_{-\infty}^{\infty} |x(t)|dt < \infty $

However, if you are just interested in the engineering or physics point of view, you can use the Dirac delta "function" as a trick that works as the measure in order to get the operations you are interested to do. For the sinc function and square integrable functions, a Fourier Transform over the $\mathcal{L}_2$ space, which is the space of square integrable functions, is defined so that you can get the Fourier Transforms for those functions(http://math.mit.edu/~jerison/103/handouts/fourierint1.13.pdf), but if you are just interested in the engineering point of view, it is not important to be so rigorous mathematically.


Can I study continuous time Fourier Transform and treat the rest as special cases

Yes you can, somehow. But that won't be the easiest road.

When teaching basic DSP, I often see students afraid to "learn all the formulae" and properties pertaining to the Fourier transform, Fourier series, discrete-time Fourier transform, $z$-transform and the Discrete Fourier transform with the fast
 Fourier transformation or FFT, up to time-frequency transformations.

They share a lot of common traits: energy conservation, sorts of orthogonality, scalar product formulae, convolution turning into multiplication, Parseval-Plancherel formulae, etc. but also important (and often subtle) differences.
.
One tool to unite such concepts is sometimes called the Pontryagin duality. For a basic overview, the wiki page is ok:

  Pontryagin duality places in a unified context a number of
  observations about functions on the real line or on finite abelian
  groups

For instance, you can view the four primal and dual domains for Fourier "transformations" as in the table:

\begin{array}{lcc}
\textrm{Transform} & \textrm{Original domain}& \textrm{Transform domain}\\
\textrm{Fourier transform} & \mathbb {R}  &\mathbb {R} \\
\textrm{Fourier series}& \mathbb {T} & \mathbb {Z} \\
\textrm{Discrete-time Fourier transform (DTFT)} & \mathbb {Z}  & \mathbb {T} \\
\textrm{Discrete Fourier transform (DFT)}& {\displaystyle \mathbb {Z} /(n)} & {\displaystyle \mathbb {Z} /(n)}
\end{array}

Its provides several ways to derive properties from one domain to another.
But this tool dives deep into complicated concepts such as category theory, see for instance Notes 2: The Fourier transform (Terry Tao) or Category theory applied to Pontryagin duality.

Knowing properties about one transformation, you can derive properties of others, but for that you will have to master some advanced mathematics. There is no free lunch, unfortunately.

You would better learn how to remember a few formulae for a handful of transformations


The sensible way here is to realize that each complex sample has a magnitude $|z|=\left\lvert \Re z + j \Im z\right\rvert = \sqrt{\left(\Re z\right)^2+\left(\Im z\right)^2}$.

In your case, the maximum non-clipping amplitude is $2^{15}$ (the maximum-power samples are on a circle around the origin).

Therefore, you'd calculate the magnitude of a sample, and compare that to $2^{15}$ to get your relation to full scale.


Neural networks in system identification - What type of activation functions?

Looks like you've done a lot of work on your projects. As @MarcusMüller said, by far the majority of people start with ReLU and go from there. It doesn't have the "vanishing gradient" problem that tanh has for example. All your questions are open ended but common for designing neural networks. There are so many "nobs to turn" to try and make your network be the best it can be.

You could adjust some parameters, train, test, and compare. But then days or weeks have gone by and you'll still be tuning your parameters. One popular method to choosing "good" parameters is called hyper parameter optimization. This blog post does a good job explaining it: https://blog.floydhub.com/guide-to-hyperparameters-search-for-deep-learning-models/. The idea is that you randomly sample your parameters space (for example, number of neurons, a few different activation functions, different learning rates, number of epochs, etc.) and in search for the a good set of parameters. I hope this helps!


You usually use a DC notch filter. This is a recursive filter that should kill a very narrow band of frequencies around DC. One possible implementation is a first order DC notch with transfer function

$$\frac{1 - z^{-1}}{1 - \lambda z^{-1}},$$

where $\lambda$ is some number very close to 1, say 0.99. A to demonstrate this in MATLAB,

b = [1 -1];
a = [1 -0.99];
freqz(b,a);

Notice that the phase response around DC gets pretty non-linear. Such is the nature of these recursive filters.

The closer you make $\lambda$ to 1, the narrower the notch gets and the more distorted your phase response becomes around DC.

For ECG signals, this will usually do. Experiment with $\lambda$ values that will suit your needs. Sometimes, if you make the notch too narrow, you won't get all off the drift, so 0.99 is usually a good starting point.


Suggested Preprocessing methods for OCR on Circular Images

Here is what I did for a client (What you are asking is the same).

Assuming that you have access to certain type of a pattern on the image (or the center of the hole), you could always detect the template to obtain the location of a possible unwarp:

Note that in the transformed image, two region of interests are defined and the region within which we would like to read the characters is the torus lying between two circles. Using this information, we could warp this image using an inverse polar transform to obtain:

I also provide the MATLAB code for this transformation at this file-exchange post.

In the next step, thanks to the template matching, the characters to be read always lie in the local region defined in the image above. Now, if you have appropriate illumination the character segmentation is nothing but a trivial set of image processing operations:

The connected components obtained eventually are fed into a neural network or an SVM, first to train and then to classify. This part is composed only of a basic machine learning task, which I won't go into the details, here.

Don't forget to take good care of the illumination, as it is almost solving the problem half way.


Where does $\frac{N}{2}$ came from in approximating an N-point DFT?

The common formulation of the forward DFT preserves energy (Parseval's theorem).  This means that a longer constant magnitude sinewave input to a DFT, which has proportionally more energy, must be represented by a proportionally larger value in the DFT result, preserving that greater energy.  Thus a factor of N, which scales with the length of the input.

For strictly real input, this energy is divided evenly between the positive frequency bin (k) and the negative frequency bin (N-k) of a DFT result (for a DFT of length N). These two complex result bins are always complex conjugates of each other (thus equal in magnitude) for strictly real input to a DFT.  Thus another factor of 1/2 is included to do this even energy split.

N and 1/2 results in a factor of N/2 which has to appear somewhere in the result of an energy preserving formulation of a forward DFT.

Note that the 1/2 does not appear for bin 0 or for bin N/2 of even length DFTs, as these have no complex conjugate image (or rather they include both of their complex conjugates summed together into a single result bin).

Also note that this Sinc function result is an approximation.  The actual result is closer to the sum of two Dirichlet functions (also called periodic Sinc functions).  See here for more details.


Kalman Filter - Optimal Way to Handle "Derived" Measurements?

In ideal world you'd have the correct model and use it.
In your case, the model isn't perfect.
Yet the steps you're suggesting are based on a knowledge you have about the process - which you should incorporate into your process equation using your dynamic model matrix:

The classic and correct way given F matrix is built correctly according to your knowledge.

"Extra" prediction step will yield nothing, since $ {F}_{ik} = {F}_{ij} {F}_{jk} $ and if you reduce the time frame you should alter $ Q $ and $ R $ accordingly which should get you at the end of the chain of small steps the same $P_{k \mid k - 1}$.

If you don't measure V you'll must "Estimate" it somehow. Yet by definition, if your case falls under Kalman's assumptions using Kalman's filter would yield best results.

All in all, stick with the "Classic".


"Is there a way to measure frequency (detect pitch) better than FFT, that is, with better resolution in less acquisition time?"

yes there is.  or are.  there are multiple better ways to do musical pitch detection in real time that are far, far better than running an FFT.

consider :

Average Magnitude Difference Function (AMDF)

$$ Q_x[k] = \sum_n |x[n] - x[n-k]| $$

Average Squared Difference Function (ASDF)

$$ Q_x[k] = \sum_n (x[n] - x[n-k])^2 $$

Autocorrelation Function (AF)

$$ R_x[k] = \sum_n x[n] x[n-k] $$

note that i am playing fast-and-loose with the limits to the summation.

note also there there are no assumptions made about waveform shape or zero-crossings or other threshold crossings.  the only assumption is that when the lag $k$ is approximately a period (or two periods or some other integer multiple of the period length), $x[n]$ looks a lot like $x[n-k]$.  so the only assumption is that pitch is related to fundamental frequency of a periodic or nearly periodic (what i like to call "quasi-periodic") function.

my favorite is ASDF (and that is a thinly-veiled trade secret i just announced to everyone, but folks on comp.dsp knew that already).  these are all time domain, AMDF and ASDF look very similar and ASDF looks like an upside-down version of AF.  you are looking for nulls in AMDF or ASDF or peaks in AF which would correspond to potential period lengths of the quasi-periodic input.

here are a couple other tricks:

you can always correlate the most current $N$ samples against some $N$ samples delayed by $k$.  that way you are dealing with the most current data possible in the real-time application.
you need not calculate the correlation for every integer lag $k$.  in fact, since you like log-frequency, the spacing for larger $k$ might be bigger than the spacing for smaller $k$.
when a potential null (AM_F) or peak (AF) is found, you can compute the correlation for adjacent integer values of $k$.
between adjacent integer values of $k$, you can do interpolation to determine the peak location to a fractional-sample precision.  i won't tell you how.  use your imagination.
the whole trick (and this is the secret sauce where trade secrets and IVL patents apply) is to choose the correct peak or null when there are multiple candidates.  choosing the incorrect peak or null will result in an "octave error".  i'm not gonna tell you how to do that.  use your imagination.

send me an email and we can discuss contracting terms if you want me to design you a kick-ass pitch detector.  much better than YIN, which, in my opinion, works like shit.


Problem with simulation of the desired shape of the signal

For a square wave as used in your example phase don't do anything.
Since the fourrier series is just based on a sine wave without phase.
The fourrier series are:
$$
S_n (x) = \frac{a_0}{2}+ \sum_{n=1}^{N} A_n \cdot \sin (\frac{2 \pi \cdot n \cdot x}{P}) \qquad \text{for integer } N \geq 1
$$

Since the square wave have only the odd frequencies. you can choose for the sinewaves: 1, 3, 5, 7, 9, 11, 13 ,...

For other functions it is different. In that case you have both sine and cosine signals. Combination of both give as result the phase.
if you have an series of 7 or 8 frequencies it is allready close to a square wave.


relationship of frequency with range in Underwater acoustic channel

What determines the effective range of frequencies for a given frequency and distance in underwater acoustic channels are the frequency and distance dependent path loss, and the power spectral density (PSD) of the noise.

The path loss is given by

$$A(f,d)=d^{\kappa}\left[\alpha(f)\right]^d$$

where $d$ is the distance between the transmitter and receiver in Km, $f$ is the frequency in KHz, $\kappa$ is the spreading factor which is usually assumed 1.5 for practical consideration, and $\alpha(f)$ is the absorption coefficient which usually calculated using the Thorp's formula as (measured in dB/Km)

$$10\log_{10}\alpha(f) = 0.11\frac{f^2}{1+f^2}+44\frac{f^2}{4400+f^2}+2.75\times 10^{-4}f^2+0.003$$.

The PSD of the noise is the sum of the PSD of the noise sources in UWA channels, which are: turbulence, thermal, shipping, and wave noise. You can find the equation of each noise in the literature.

Having the path loss and the PSD of the total noise, you can find the 3-dB bandwidth as following: the narrow-band signal-to-noise ratio (SNR) is given by

$$\gamma(f,d)=\frac{P/A(f,d)}{N(f)\Delta f}$$

where $N(f)$ is the PSD f the total noise, $P$ is the transmit power, and $\Delta f$ is the frequency over which the channel is frequency-flat. For a given distance $d$, there is an optimum $f$ that maximizes the SNR (which basically minimizes the term $\frac{1}{A(f,d)N(f)}$). Let that frequency denoted by $f_d$. Then you find the range of frequencies for a given distance such that

$$\gamma(f,d) \geq \frac{1}{2}\gamma(f_d,d)$$.

This way you can find the range of frequencies that can can be used for UWA communication for a given distance.

See this reference for more information.


The "spikes" (they are quite small in magnitude, but large compared to the nearby filter taps) at the beginning and end are part of how the equiripple property is achieved.  In his book "Multirate Signal Processing for Communication Systems" (a book that I recommend quite highly), fredric harris indicates that it is sometimes advantageous to eliminate the spikes because it causes the frequency response to continue to decrease after the cut-off frequency instead of staying constant.  The price you pay is that the filter is no longer equiripple which appears to cause (see example below and comments by JasonR) the transition band to be not as sharp.

EDIT:
I created a low-pass filter in Octave using remez.  I intentionally made the transition-band very small to get it to create a decent spike.  The filter taps are plotted below.

Here is the frequency response of the filter.

As you can see, the frequency response is pretty poor due to the short transition-band, but the point is that it is an equiripple filter.  Below I plot the exact same filter, only I have zeroed out the two spike taps.

Below is the frequency response of the "no spikes" version of the filter.

The frequency response shows the continued drop beyond the cutoff frequency.  It also clearly shows the non-equiripple response in the pass-band.  This is entirely expected.  The one thing that is unexpected is that the "no spike" filter appears to be better than the "spike" filter at every point.  I assume that this is because of some flaw in Octave's remez algorithm.


Can CV or DS be seriously considered DSP roles?

in my opinion

your question seems based on a combination of thoughtful reflection and possibly, pre Alzheimer’s cognitive biases. ;)

There was an article in the WSJ not too long ago , I’m outside the paywall so I can’t provide a link, that looked at IBM’s efforts using Watson at major Cancer research centers.

The conclusion was that Watson wasn’t offering much in terms of originality. It wasn’t revealing those deep unifying subtle links in the masses of data it was being feed. It tended to find what was already known.

hard problems are still hard particularly if the way the data is represented ( labeled) isn’t relevant to what might be hidden. The threshold of what might be considered for brute force search might be advancing but NP hard hasn’t gone away.

I don’t believe I answered your question but I don’t think your question has a clear answer.

When I was an undergrad, the EE department had just the year before my junior semester had dropped land surveying as a core requirement. skill mixes change.  not many years ago i was calculating viewsheds to estimate RF coverage. Land Surveying would have probably had some use.

Are patterns truth?


If Digital Signal has frequency components then why we don't use it

Analog and digital
First off, you need to understand what an "analog" signal is, and what a "digital" signal is, how they are different, and how they are similar.
The term "analog" comes from the old distinction between "analog" and "digital" computers.  A "digital computer", even a very primitive one of decades ago, has always been more or less what we know now -- a box that manipulates "bits" (or, in very old systems, decimal digits), where the information is represented by numbers encoded in the bits.
An "analog computer", on the other hand, used voltage values in the computer to represent numeric quantities.  Ie, where the digital computer might represent the value "6" as the binary bits 110, the analog computer might represent the value "6" as a voltage of 6 volts.
(I could go on here for several paragraphs explaining more about how an analog computer works, but it suffices to say that, except in very limited applications, analog computers are no longer used -- where needed their function can be simulated with digital computers more cheaply and accurately.)
But the point is that an "analog signal" is called that because the voltage levels in the signal are analogies ("analogs") for some "real-life" phenomenon.  Eg, in an analog audio signal the voltage level is an analog for the sound pressure that would create the sound -- a higher voltage level means more sound pressure, etc.  Further (and it will be important later), the "X-axis" -- the time dimension -- is also "analog".  That is, in most cases (except where, eg, a recording has been sped up or slowed down) one second of time in the analog signal corresponds to one second in "real life".
In the image below, peaks in the signal represent higher air pressure (and valleys represent lower air pressure).

On the other hand, the ups and downs of a digital signal have no direct relationship to "real life".  Rather, they represent binary digits (hence "digital") in a numeric code used to represent whatever "real life" phenomenon the signal corresponds to.  In the diagram below, the bits (each line represents a separate bit stream) might represent characters of ASCII text, or they might represent readings from a heart rate monitor.  Looking at the data you can't even tell if it might be a pleasant sound (like the guitar), or a discordant sound (like the cymbal).
Also note that the X dimension is not an analog for anything -- an entire symphony might be transmitted as a digital signal in ten seconds, or a high-resolution photograph that was snapped in an instant may be transmitted, bit-by-bit, over a period of seconds or minutes.  (The old Mariner space probe would take hours to transmit a single high-resolution image, because the bit rate was so slow at that distance.)

But note that both types of signals can be represented by a fluctuating voltage on a wire, and, in fact, you can hook a digital signal (if the data rate is slow enough) directly up to an ("analog") audio amplifier and listen to it.  (It will come through as a noisy squeal).  So at this fundamental level -- a voltage that fluctuates over time -- analog and digital signals are basically the same.
(It also needs to be noted that, from the standpoint of "information theory", the maximum amount of information that can be encoded into a signal occupying a given bandwidth is the essentially the same, regardless of whether it's digital or analog, or how it may be multiplexed.  However, while one could spend days discussing information theory and entropy I won't go into it further, because it's for the theoreticians, and because it makes my head hurt to think about it.)
Multiplexing
Next we need to understand "multiplexing".  In the US, at least, a "theatre multiplex" is a group of motion picture theatres in a single building.  Similarly, when you "multiplex" signals you combine several signals into one all-encompassing signal.  But just as you wouldn't want to go to a motion picture theatre and see 6 movies all projected onto a single screen simultaneously, when we multiplex signals we do so in a way that the receiving station can separate the combined signals back into the individual ones.
The earliest form of signal multiplexing(*) was common radio.  Several radio stations in the same area can all broadcast their signals into the same "ether", but each uses a different assigned frequency, so that, by tuning your radio receiver to a specific frequency (actually a frequency band) you can receive the "hard rock" station you want to hear, and not have to simultaneously hear the "easy listening" station your parents listen to.
This, in essence, is Frequency Division Multiplexing (FDM) -- the available spectrum is divided into frequency bands and each signal is combined with a "carrier" frequency to occupy the appropriate frequency band.  (The original signal, before multiplexing, is referred to as a "baseband" signal, because it has not been combined with a carrier.  This sort of contrasts with "broadband", where multiple signals have been multiplexed, but unfortunately the term "broadband" has multiple confusing definitions.)
FDM used to be used extensively in the phone system, for long-distance lines.  A single pair of wires (and associated amplifiers), which might have a frequency response from 20Hz to 100KHz, carried a dozen or so audio signals, each at an assigned frequency just like broadcast radio or TV.  Later coaxial cables were used, increasing the available bandwidth to several MHz and allowing several hundred audio signals per cable.
(A significant drawback of this scheme was that many amplifiers were required -- one every ten miles at least, and one every mile for some schemes.)
Now, it should be clear that, given that there's no fundamental difference (electrically) between an analog signal and a digital signal, digital signals can indeed be broadcast using FDM.  In fact, if you have a cable modem, the cable, with over a hundred analog channels, is frequency multiplexed (just as over-the-air TV is).  Your cable modem signal occupies (sharing with your neighbors) one of those "analog" channels.  So your cable modem signal is being transmitted via FDM (though with several added twists).
Now for Time Division Multiplexing (TDM).
Time division multiplexing involves what it name implies -- the time axis of the combined signal is divided into discrete segments, and different individual signals are transmitted in the different segments.
TDM actually goes back to before regular digital data.  The Bell Telephone company in the US experimented in the 50s with chopping up analog signals into lots of thin slices and interleaving the slices from several different analog signals together, allowing the several signals to be sent over a single pair of wires without having to resort to FDM.  So in this case analog signals were transmitted via TDM.
But doing this was quite complicated.  The analog signal needed to be "sampled" rapidly enough that the signal didn't change much between samples.  For an audio signal with a maximum frequency of about 3.5KHz (the range for common telephone lines) sampling must be done at at least twice that frequency (so about 7 thousand "slices" per second) and then those slices must be interleaved with others (maybe five or ten).  (This was known as "pulse amplitude modulation".)  At the other end the slices need to be separated, and the missing space between slices filled in by "smoothing" the signal.  And all this needs to be done without seriously changing the amplitude (sound volume) of individual slices.
Of course, the above scheme was very difficult to get right and didn't really provide much advantage relative to FDM over the same wires, so when true digital capabilities came around they dropped it like a hot rock.
Why digital is different
So far I've been showing how there's really no difference between analog and digital.  Now I'll reverse course -- digital is really different in some very important ways.
First off, digital signals don't degrade like analog signals.  Old analog telephone networks introduced a lot of distortion into the sound as it was transmitted cross-country and reamplified many times (as often as once every mile).  A digital signal does develop distortion in transmission, but since it's digital it can be "regenerated" at the end of each hop.  This alone made TDM much more practical for digital, since there was no longer a need to maintain accurate signal amplitude, and the cables and amplifiers become simpler and less expensive (and farther apart).
But that's not all.
Consider, for instance, if you use Skype to talk to someone on the Internet.  Your voice is first converted to an analog signal by the microphone in your computer.  But it's then immediately digitized, and the digitized bits are "buffered" -- stored up awaiting transmission.  Algorithms in Skype will examine the digitized sound and discard periods of silence (since there's no need to transmit silence), and they'll otherwise process the sound to minimize the number of bits needed to represent it.
Eventually (after a seeming eternity of computer instructions, but in an instant from your ear's point of view) a group of sound bits will be formed into a block, have the appropriate internet headers appended to the front, and be queued for transmission over your Internet connection.  But, since you happen to be simultaneously using your web browser, the sound data may have to wait a few more milliseconds for your web request to be transmitted.
On the receiving end, Skype accepts the packet of transmitted data from the receiving computer's Internet connection and sticks it back in a queue, taking care to insert it in the right order in the queue based on a "timestamp" in the data packet.  (This extra care is necessary because data packets may be delivered out of order in some cases.)  Then another part of the Skype program removes packets from the front of the queue and "schedules" them to play back the recorded sound at the appropriate time based on the timestamps.  So the received data is not played immediately, but instead waits (perhaps a few more milliseconds) until "the time is right" to be played, so that the time dimension of the original data (and hence the fiction of continuous voiced transmission) is maintained.
So what happened here?  The voice data from Skype was time-division multiplexed with your web browser requests, and, once out on the Internet trunk lines, both were time-division multiplexed with thousands of other messages.  But, critically, the multiplexing was not based on a fixed schedule (eg, dividing the time dimension into 10 millisecond slices and "dedicating" the channel to transmitting each signal for 10 milliseconds), but instead an "adaptive" algorithm was used, where each signal only occupied the channel for as long as it needed.  Even though this technique requires that each data packet contain an address to identify it's destination (whereas with a fixed multiplexing scheme the destination would be defined by the slice location), it's still far, far more efficient than a scheme that relies on fixed slices, either in time (primitive TDM) or frequency (FDM).
Summary
Now note the differences between TDM (as used in networking systems) and FDM.  FDM has the advantage of being relatively simple (witness standard AM radio).  It can carry analog or digital signals equally well, but, because it does either equally well, it offers no great advantages to digital.
TDM is only marginally practical for (undigitized) analog signals -- the complexity and potential loss of signal quality are major negatives.  But digital, it seems, was "made for" TDM -- the data is already "sliced", and accurately maintaining amplitude or time reference is not required.  And since "dead air" is represented by a complete absence of data packets, not by simple silence on a phone line, "adaptive" multiplexing, where only "useful" data is transmitted, and not on a rigid schedule, is much more practical.  For these reasons much more digital data can be transmitted in a given bandwidth using (adaptive) TDM vs FDM.
So it's not that digital CAN'T be transmitted using FDM and analog CAN'T be transmitted using TDM, but rather there are practical reasons why doing so would be more expensive, less efficient, and less reliable.
(I need to add that even given all the above FDM is used for digital data in some circumstances.  For instance, the link to/from communications satellites is very "broadband", and no single data stream could practically use the entire bandwidth.  So the overall channel width is sliced into frequency bands, with each occupied by a separate data stream.  In this case, since the digital data stream is already time-division multiplexed and there is very little "dead air", the advantage of further TDM would be slight.)

(*) There were some crude TDM schemes used for telegraph signals as far back as the 1850s, but that's mostly lost in the mists of history.


Bandpass butterworth filter in python is not working


What are the effects on filtering upsampled data without interpolation?

How "bad" is this?

It's not too good.  What you are doing is equivalent to upsampling by a factor of five (i.e. inserting four zeros in between every sample), and then filtering it with a filter that has the impulse response [1, 1, 1, 1, 1].  The picture below shows the frequency response of this filter.

The signal aliases introduced by the upsampling are centered at 0.4 and 0.8, so the good news is that the filter has a null at the center of the aliases.  The bad news is that unless the signal has a really narrow bandwidth, the signal aliases will be much wider than the nulls, and a lot of alias energy will distort your signal.

  Is my only option to design an interpolation filter during upsampling, and if so, can anyone lead me to good references on how to design such a filter?

If you want to understand how to design a low-pass filter then I would suggest learning about windowed sinc functions.  If you just want to get a filter that works then I would recommend installing Octave and using the fir1 function to design the filter.


Is moving median always less sensitive to outliers? Sometimes. It will work if you have a very short spike (preferrably shorter than the median/average sample size). However, if you have a large spike, then taking the median won't help eliminate the spikes.

Long spike

I have illustrated this using a dataset (with sample size = 5, taking into account 2 left entries and 2 right entries):

Orig. Median Avg.
Data  Filter Filter
-----------------
1
2
3      3    3
4      4    4
5      5    11.6
6      6    19
40     40   26.2
40     40   33.2
40     40   33.4
40     40   27
7      9    20.8
8      9    14.8
9      9    9
10     10   10
11
12

Notice how the the median of the all the 40s is 40. For example, take the 1st 40. The left values are 5,6 and the right values are 40,40, so we get a sorted dataset of 5,6,40,40,40 (the bolded 40 becomes our median filter result).

Short spike

You also wanted an example for the median filter to work. So, we will have a short spike. Try this:

Orig. Median Average
Data  Filter Filter
-----------------
1
2
3      3    3
4      4    10.8
5      5    11.6
40     6    12.4
6      7    13.2
7      8    14
8      8    8
9
10

You can see that the median filter did filter out the single large spike, while the spike skewed the average filter's results for entries 4,5,40,and,6.


Calculating cutoff frequency for Butterworth filter

Everythng is alright. You made a mistake in the final calculation.

$$\Omega_c^{12} = \frac{{(0.2\pi)}^{12}}{0.25893} \implies \Omega_c = \sqrt[12]{\frac{{(0.2\pi)}^{12}}{0.25893}} =\frac{0.2\pi}{\sqrt[12]{0.25893}} = 0.7032$$

which is the desired result.


Measuring microphone Signal to Noise ratio (SNR) with A-weighting

IEC 801-21-20 Defines sound pressure as the RMS of of the instantaneous pressure, and as far as I know all meters use RMS.  The problem you have in your implementation is that you are normalizing the peak value of the active region but measuring the RMS of the inactive part. As you described, for sine waves the RMS is 3dB less than the peak.  You could take the Difference between the RMS in the active to the inactive regions.  Or you could assume that the active RMS is -3dB and use that instead.  Otherwise you should be good to go.


Confusion in the derivation of convolution formula

Your confusion stems from abusing the time invariance property.

Let's restate it for a simple 1D system:

$$ \text{if} ~~ T\{ x(t) \} = y(t) \implies T\{x(t-d) \} = y(t-d) $$ is correct but the following

$$ T\{ x(-t) \} = y(-t) ~~~\text{is wrong in general} $$

Now in your derivation you replace the input $\delta(t-\tau)$ with time reversed input $\delta(\tau-t)$ and claim the same on the output as $h(t-\tau) = h(\tau-t)$ which is not true as shown above.

The solution comes simply by recognizing the fact that $\delta(t-\tau) = \delta(-(t-\tau)) = \delta(\tau-t)$ and then proceed only with the $\delta(t-\tau)$ case as Marcus have stated.


How to verify if my cross-correlation algorithm works?

There are two issues / concerns with this approach in that you may not be getting what you want. The primary one is the equivalent of post detection averaging by using the post-processed results from a spectrum analyzer, and the second is that the result of your complex conjugate multiplied FFT is actually in the time domain since you started in the frequency domain (so the FFT processing is actually converting to the the time domain where you perform the multiplication.-- I suspect but am not certain that your goal is to observe the frequency domain spectrum after cross-correlation of your signals?).

Most significantly, spectrum analyzers "detect" the signal in that the result is dB magnitude, and regardless of the dB conversion which could easily be undone, all phase information is lost such that any processing gain through coherent averaging can no longer be done. This is equivalent to the difference in adjusting the resolution bandwidth (RBW) versus the video bandwidth (VBW) of the spectrum analyzer . Since you are already using a spectrum analyzer, playing with those knobs will actually provide a great demonstration of pre-detection versus post-detection averaging, which will give you much further insight into what you are trying to accomplish. For that reason, I will first go into that in more detail of that, starting with a simple functionally equivalent block diagram of a spectrum analyzer:

In this simplified but functionally equivalent view of a spectrum analyzer an input signal (in this diagram a 70 MHz tone) is frequency translated by the voltage controlled oscillator (VCO) and mixer to be swept through the Bandpass Filter with resolution bandwidth RBW. It may be helpful to envision the sweep rate as being extremely slow in comparison to the inverse of the filter bandwidth, such that the operation for explanation here is similar to stepping the VCO through each frequency one at a time and statically processing each result (pixel on the display) before stepping to the next frequency. This will avoid us getting into constraints on the sweep rate which exist. Thus at any given moment, the power detector is converting the entire power that is within the bandwidth RBW to a power level. This power level passes through a low pass filter with bandwidth VBW, in converted to dB and then used to control the vertical position of the pixel on the screen representing the power level of our signal under test. Similarly the ramp rate controls the horizontal position representing the frequency at that moment in time.

By adjusting the bandwidth RBW we reduce the total power that would be presented to the power detector. If the signal is centered on our tone (in this case 70 MHz), the power of the tone would dominate and we would see no change in the power detector output regardless of the bandwidth of RBW, as long as the bandpass filter remained centered on our signal under test, and assuming the signal under test was a pure tone. In this case the apparent bandwidth of that signal, as displayed in the diagram above is actually the bandwidth of RBW shown as the single tone is swept through our filter during the horizontal trace. This would be very apparent by adjusting RBW and observing the result. More notably you would see the noise floor everywhere else go down at $10Log(N)$ where N is the decrease in RBW. (If you half RBW, the noise floor goes down by 3 dB). I believe this is a direct demonstration of the noise floor reduction you seek to that the extent correlation has mathematical similarity to averaging, as in this case you would observe that the signal level would remain unchanged. In this case, this is the result of "Pre-detection" averaging which is exactly what a bandpass filter is doing (with the averaging operation mathematically translated to the center frequency of the filter). In the plot shown in the figure above, the noise floor is approximately -90 dBm. Both the noise floor and the width of the tone suggest that the RBW is about 1 MHz, so I will assume that was the case. If we adjusted RBW to 500 KHz, we would see the noise floor drop to -93 dBm while the peak of the tone would remain unchanged.

Similarly if you adjust the bandwidth of the lowpass filter (VBW), you are simply smoothing (averaging) the noise power already measured but not decreasing it! In the plot above we would see the "noise on the noise" get smoother, but it will remain at the -90 dBm level. Thus we are simply averaging the result of the noise without reducing it.

To accomplish what you want, I would consider cross-correlating the signals directly in the time domain. The complex conjugate multiplication of the FFT of these time domain signals would then represent the spectrum of the cross correlation, if that is what was ultimately desired.


Fitting an IIR filter to a complex transfer function

As I mentioned in a comment, I think that the equation error method is a very good starting point for designing an IIR filter in the frequency domain with prescribed magnitude and phase responses. For an explanation of this method have a look at this answer and the references there.

Two problems associated with the equation error method are

stability of the designed filter is not guaranteed
the frequency domain error is weighted with the squared magnitude of the denominator polynomial of the designed filter, which means that the approximation is better in certain frequency regions than in others, but this weighting is a result of the design process and is usually undesirable.

The first problem can be tackled by adapting the total delay of the desired frequency response such that it matches with what a causal and stable filter with the chosen orders of numerator and denominator polynomials can achieve. So the desired complex frequency response $D(\omega)$ is changed to

$$\tilde{D}(\omega)=D(\omega)e^{-j\omega\tau}\tag{1}$$

with some delay parameter $\tau$ (which can also be negative). Optimizing the value of $\tau$ will help designing a stable filter with the smallest possible error. The only difference with the original desired response is an added (or subtracted) delay (i.e., a linear phase component).

One option to address the second problem is to solve a sequence of equation error problems with an adapted weight function, which is the basic idea of the Steiglitz-McBride method.

I've implemented the equation error method for the frequency domain design of IIR filters, and I've also added the possibility to iteratively approximate the optimal least squares solution, as described above. Note that convergence cannot be guaranteed, but in practice convergence is usually achieved after a few iterations. Also note that the resulting filter may be unstable, but this problem can be solved by adjusting the total delay of the desired response, as explained above. The Matlab/Octave code can be found at GitHub.

Here is a simple filter design example. The desired response is a low-pass filter with a linear phase response in the pass band. The delay is desired to be an integer number of samples, but the exact value is to be optimized to get the best approximation by a stable filter. The numerator and denominator orders are chosen to be $M=N=10$. Experiments showed that the best value for the desired pass band group delay is $28$ samples:

w = pi*[linspace(0,.2,20),linspace(.25,1,75)];   % frequency grid
D = [exp(-j*w(1:20)*28),zeros(1,75)];           % desired response
W = [ones(1,20),10*ones(1,75)];                  % error weighting
[a,b] = eqnerror(10,10,w,D,W,6);                 % 6 iterations of equation error method

The figure below shows the design result. Note that the pass band group delay approximates the desired $28$ samples. The designed filter is stable and has a maximum pole radius of $0.98$.


Conceptualising the continuous time unit impulse function as derivative of unit step

I didn't look at how the function $\delta_\Delta(t)$ is defined in that lecture, but the convention is that it's a rectangular pulse of area being $1$. The definition is :
$$\delta_\Delta(t) = \begin{cases}{ \frac{1}{\Delta} ~~~, ~~~ 0 < t < \Delta \\ ~0~ ~~~,~~~\text{otherwise} } \end{cases} $$
As you can see the area of that pulse, $A(\Delta) = \frac{1}{\Delta} \times \Delta = 1$ for any finite value of $\Delta$. According to calculus, that area is still $1$ even when the parameter $\Delta$ takes on inadmissible values of $0$ or $\infty$ through a limiting process. In particular the continuous Dirac delta function is defined as $\delta(t) = \lim_{\Delta \to 0} \delta_{\Delta}(t)$ ; when the limit of $\Delta$ goes to zero.
The rectangular pulse of unit area, $\delta_{\Delta}(t)$,  is only one of many other possible ways of defining $\delta(t)$ through a limit of functions process.


Zero Phase Filter: Determining Initial Conditions for Forward Backward Filtering

For anyone who is interested, i coincidentally found a paper describing the method implemented in matlab's filtfilt.m. A link to the paper is attached. At least to my understanding matlab's filtfilt.m doesn't implement the Gustafson algorithm.

Sadovsky, P.; Bartusek, K: Optimisation of the Transient Response of a Digital Filter, Radioengineering Vol. 9, No. 2, 2000