How do I interpret the DCT step in the MFCC extraction process?

You can think of the DCT as a compression step. Typically with MFCCs, you will take the DCT and then keep only the first few coefficients. This is basically the same reason that the DCT is used in JPEG compression. DCTs are chosen because their boundary conditions work better on these types of signals.

Let's contrast the DCT with the Fourier transform. The Fourier transform is made up of sinusoids that have an integer number of cycles. This means, all of the Fourier basis functions start and end at the same value -- they do not do a good job of representing signals that start and end at different values. Remember that the Fourier transform assumes a periodic extension: If you imagine your signal on a sheet of paper, the Fourier transform wants to roll that sheet into a cylinder so that the left and right sides meet.

Think of a spectrum that is shaped roughly like a line with negative slope (which is pretty typical). The Fourier transform will have to use a lot of different coefficients to fit this shape. On the other hand, the DCT has cosines with half-integer numbers of cycles. There is, for example, a DCT basis function that looks vaguely like that line with negative slope. It does not assume a period extension (instead, an even extension), so it will do a better job of fitting that shape.

So, let's put this together. Once you've computed the Mel-frequency spectrum, you have a representation of the spectrum that is sensitive in a way similar to how human hearing works. Some aspects of this shape are more relevant than others. Usually, the larger more overarching spectral shape is more important than the noisy fine details in the spectrum. You can imagine drawing a smooth line to follow the spectral shape, and that the smooth line you draw might tell you just about as much about the signal.

When you take the DCT and discard the higher coefficients, you are taking this spectral shape, and only keeping the parts that are more important for representing this smooth shape. If you used the Fourier transform, it wouldn't do such a good job of keeping the important information in the low coefficients.

If you think about feeding the MFCCs as features to a machine learning algorithm, these lower-order coefficients will make good features, since they represent some simple aspects of the spectral shape, while the higher-order coefficients that you discard are more noise-like and are not important to train on. Additionally, training on the Mel spectrum magnitudes themselves would probably not be as good because the particular amplitude at different frequencies are less important than the general shape of the spectrum.


Cross-correlation of time compressed 1D signal

With known time scaling, you can resample to a common rate and then cross-correlate.  It is natural to think you might need to go to the higher of the two rates, but that does not increase information beyond the lower rate, it only interpolates the information. This is analogous to zero-padding an FFT to interpolate the resolution of a peak.

To determine unknown time scaling, cepstral processing may provide some insights.  I'm afraid it is not my area, but it sounds related to what you are describing.

Essentially, you've got two dimensions of uncertainty (start time, rate) across which you'd like to search to find the best correlation.

You can reduce it to one dimension by considering the magnitude of the two signals' spectral data. The lower rate data should then look like a stretched version of the higher rate data.


2D Convolution in MATLAB Causes Artifacts (Boundary Issues)

If I understood right this is a duplicate question but I cannot find the link let me re-type the answer. Given that you want to implement a 2D discrete convolution of two images of sizes $N_1 \times N_2$ and $M_1 \times M_2$ by the method of DFT multiplication, then in order to avoid time (space) domain aliasing due to circular convolution implied by the inverse DFT stage, the sizes of the forward DFTs should be properly selected (in other words original time/space domain signals should be padded with zeros to increase thier sizes to overcome the aliasing).

Particularly the 2D forward and inverse DFT sizes should be selected as:
$$ L_1 \geq N_1 + M_1 -1$$ and $$ L_2 \geq N_2 + M_2 -1$$ in order to avoid circular artifacts and get the exact convolution.

Select the central portion of the resulting convolution to get the final image the same size of the original image.


When modelling noise,how do they decide whether additive or multiplicative?

It will depend on what you are modelling in your system. Consider an active radar system. There are both additive noise components and multiplicative noise components. What distinguishes them in this case:

If I increase the transmit pulse power by a factor of $K$ and the SNR improves by a factor of $K$ then the noise is additive
If I increase the transmit pulse power by a factor of $K$ and the SNR does not change then the noise is multiplicative

To give some examples to clarify:

Thermal noise at the input of the ADC is generally additive - noise level does not generally change with power of the received signal, so SNR will change with transmitted/received power.
Clipping in the ADC (distortion) is non-linear - it is more multiplicative than additive. Increasing the receive power won't affect the SNR (assuming clipping still occurs).
Phase noise in the ADC (timing errors) is multiplicative
Speckle noise - Presence of a lot of random scatters in radar resolution cell is usually modeled as multiplicative - increasing the transmit power does not reduce the speckle effect.

Another common example is in active Sonar. In near range there is a lot of reverberation (volume reverberation), but as the transmitted pulse travels away from the sonar the reverberation drops and the system performance becomes limited by the ambient noise. In these two cases the sonar is operating in reverberation limited mode and noise limited modes respectively.

Reverberation - Multiplicative - Increasing transmit power offers no improvement in SNR.
Noise Limited - Additive noise - Increasing the transmit power improves the SNR.

This assumes a very simplistic model of underwater sound propagation. You will have other wave propagation effects mostly due to the sound velocity profile (Reverberation at surface/bottom, sound channels, convergence zones, etc).

So within a single system you can have both. In what mode of operation will determine which type of noise dominates. It should also be noted that the noise may not be entirely multiplicative or entirely additive. You can have noise which contains elements of both. In this situation you will get some SNR improvement but it may not be linearly related to the increase in transmit power.


Complex sampling does not "break" Nyquist.  IQ quadrature sampling produces twice as many bits per second of information (at the same sample rate for real or complex samples), and the 90 degree phase offset between the I and Q channel in those bits provides extra information about the spectrum.
One typical example to demonstrate aliasing is that the samples of DC and the samples of a sinusoid of a frequency at the sample rate look the same.  A single channel of samples at 40 MHz would alias 0 Hz and a 40 MHz sinusoid together.  Only if you strictly bandlimit the input spectrum to below Fs/2 in bandwidth would you prevent this aliasing.  If you feed a strictly real data vector to an DFT, the result would be conjugate symmetric, half the FFT result would be redundant, and the upper half can be called the negative frequency image.
But with quadrature sampling (using a 40 MHz sampling clock plus a 90 degree phased shifted copy), a DC signal and a 40 MHz sinusoid would be clearly differentiated.  I and Q would be the same for DC, but different for 40 MHz, e.g. if I is the same, the Q would be different, or vice versa.  Therefore the spectrum between half the IQ sample rate and the sample rate will not be aliased with spectrum below, as is the case with single channel sampling.  If you feed a DFT using the IQ data as a complex input vector, the lower and upper half of the result can be non-symmetric, thus represent independant information, thus twice the bandwidth (between 0 and Fs), compared to a strictly real DFT of data at the same sample rate.


Just looking at the Google results I find:

FGComp2013 which states:

  This challenge will target fine-grained classification, i.e. classification among categories which are are both visually and semantically very similar. This is a very difficult regime which is even challenging for humans without careful training, and is critical for establishing a more detailed understanding of the visual world.

The Stanford Dogs Dataset has its primary reference which states:

  First, being a fine-grained categorization problem,
  there is little inter-class variation. For example the basset
  hound and bloodhound share very similar facial characteristics
  but differ significantly in their color, while the
  Japanese spaniel and papillion share very similar color but
  greatly differ in their facial characteristics.

So the upshot is that "fine-grained" means that the classes that are trying to be distinguished between are visually very similar.

In terms of a single image, it means that it is hard to see the differences across the entire image.  In that case, your top image is fine-grained. The dog and bird images are not.  Though if you zoomed into the hair of the dog or the feathers of the bird, I expect they would be.


I don't think it's possible to deduce the sampling rate from this information. The sampling rate could be anything, since FFT is sampling rate agnostic. Sampling rate is pretty much independent of all other parameters which scale directly with it.

Concerning the bandwidth, it's not clear either. You're specifying pilot distance, but we don't know what the units are and we never see how those variables get used. Is there more to the code you're studying?

Looked at your full code, and the sampling rate is 8 kHz, since that's what frequency is used in the wavwrite command. Will have to get back to bandwidth.


Reconstructing audio from a given spectrogram and phase data

A spectrogram is a graphical time-scroll representation of the frequency components in a signal. Sometimes it's called a "waterfall", because that's what it looks like as it is computed and plotted in real time. The core math behind this is the FFT.
A spectrogram generally shows the strength (or magnitude) of a range of frequency components in a signal (e.g. DC to 20kHz). The FFT has a complex output, one block of values for each block of input. The magnitude of these complex numbers is likely (somehow) stored in the "spectrogram" file, and the phase is stored in the "phases" file.
The trick to reconstructing the original signal from a spectrogram is to inverse the spectrum, and the math behind that is the IFFT.
You will need the magnitude and the phase data from the files to recreate the complex spectrum and then apply it to the IFFT to get the time domain signal back.
You will have to find out

what math was applied to calculate the spectrogram: FFT block size, any log-scale (dB?), any windowing, any sub-sampling etc... and how the data is formatted and organized in the file.
how the phase is calculated and stored (block size, degrees/radians, number format etc...)

As a commenter pointed out, you will have to find out all these details from whoever created the file. But beware that even then it might not be possible to recreate the original signal.
In that case, if you can specify how to generate the spectrogram with different parameters (no windowing, no sub-sampling etc...), and then (re)generate these files, you may stand a chance.
Good luck. Since this is a bit like forensics, you will need it.


How to extract pixel-specific time-stream from mp4 video

In video, pixels are not usually encoded independently; you will need to decode full frames. If you are going to run several experiments on a single video, I would suggest that you decode the video once frame-by-frame and write the border pixels into a file in your own raw format for easy access. This should give a manageable file size perhaps of a few gigabytes for a movie-length video.


How to check the FFT results of a sine wave?

You need to plot the magnitude of the output of the FFT.  I'm not familiar with your programming language, but in Python you would use something like plot(abs(fft(a))).  For a silent input, the output should be all zeros.  For a sine wave input, you should see two spikes:

For a real signal, the spikes will be symmetrical from left to right.  If you're doing a real FFT, though (which is more computationally efficient) you'll only get the left half of the plot as your output, since it ignores the redundant mirror image.

If the frequency is higher, the spikes will be closer to the center.  If the frequency is perfectly in sync with the chunk size, the spike will only be one point wide and everything else will be exactly 0.  Otherwise it will have a tapering "skirt" like above.


FFT of a $N$-length real sequence via FFT of a N/2-length complex sequence

The following Matlab Excerpt corrects your mistake. You had a problem with the DFT length N of the composed sequence... Note that I used a modulus operator for simplicity, you should replace it with indexing for efficiency.

%% Perform Butterfly

%k = 0:(N/2 -1);
%W = exp(-2*pi*1i*k/N);
%X_even =     0.5*(X(k+1) + conj(X(N/2 - k )));
%X_odd  = -1i*0.5*(X(k+1) - conj(X(N/2 - k )));
%G = X_even + W .* X_odd;

k = 0:N-1;
W = exp(-2*pi*1i*k/N);
X_even =     0.5*(X(mod(k,N/2)+1) + conj(X( mod(- k,N/2)+1 )));
X_odd  = -1i*0.5*(X(mod(k,N/2)+1) - conj(X( mod(- k,N/2)+1 )));
G = X_even + W .* X_odd;

fN = sample_rate*(0:N-1)/N;   % new frequency vector data length N...

Mathematically speaking, the line where you made the implementation error is the Even-Odd combination stage:

$$ X[k] = X_e^{N/2}[k] + e^{-j\frac{2\pi}{N}k } X_o^{N/2}[k] ~~~,~~~\text{ for } k=0,1,...,N-1 $$

Where $X_e[k]$ and $X_o[k]$ are $N/2$ point DFTs of the even and odd indexed samples of the original length $N$ real sequence $x[n]$. And $X[k]$ is the N-point DFT of $x[n]$. Note that frequency index $k$ ranges through $0,1,..,N-1$. But those $X_e[k]$ and $X_o[k]$ are half length (but periodic!) sequences. Hence for $k> N/2-1$ they repeat. Remeber that DFT seqeuences are inherently periodic: $X_e[k+N/2] = X_e[k]$. In order to use the base length of those half-length periodic sequences, you can utilize the modulus operator:

$$ X[k] = X_e^{N/2}[(k)_{N/2}] + e^{-j\frac{2\pi}{N}k } X_o^{N/2}[(k)_{N/2}] ~~~,~~~\text{ for } k=0,1,...,N-1 $$

where $(k)_{N/2} = \mod(k , N/2)$ is the modulus of $k$ wrt $N/2$.
Argument of every DFT sequence should circularly wrap around...


Signals and systems book with linear algebra approach

I'm not particularly current, my junior level signals and systems text was

  Lathi, Bhagwal. Signal, Systems, and Controls. Intext, 1973

which had a some linear algebra, primarily because it had a control engineering bent.  The last chapter on mixed continuous-discrete systems was a lot to digest. I have come to appreciate it more over the years.

I used an older version of Chen in graduate school and it is very much from a linear algebra perspective, but is not aimed at the same audience Lathi is.

  Chen, Chi-Tsong. Linear system theory and design. Oxford University
  Press, Inc., 1998.

Kailath is a bit harder than Lathi,

  Kailath, Thomas. Linear systems. Vol. 156. Englewood Cliffs, NJ:
  Prentice-Hall, 1980.

and Luenberberger is very readable.

  Luenberger, David G. Introduction to dynamic systems: theory, models,
  and applications. Vol. 1. New York: Wiley, 1979.

In signal processing,

  Oppenheim, Alan V., and Ronald W. Schafer. Discrete-time signal processing: Pearson new International Edition. Pearson Higher Ed,
  2013.

Is sort of the "standard" text book on signal processing but light on linear algebra.

Strang is good but assumes you have some linear systems background.

  Strang, Gilbert, and Truong Nguyen. Wavelets and filter banks. SIAM,
  1996.

and his book on Linear Algebra is a very good introduction.

Its hard to recommend a book without knowing your back ground, but If I were to ask the same question you did,  I would go with Chen.


Wikipedia would be a better place to get a full answer, but in short: Discrete Fourier Transform is defined over a finite set of $N$ samples $x_n$ by the formula
$$
X_k=\sum\limits_{n=0}^{N-1} x_n e^\frac{-i2\pi kn}{N}
$$
The resulting DFT coefficients contains information about the frequency content of the original signal.


Best way to resample an OFDM signal without interpolating + filtering?

I can't resample this signal using an interpolation plus filtering approach, as it should be computed in real time (I would need to interpolate the signals up to 480 MHz...).

Fallacy! Using a rational resampler doesn't require going to the common multiple of in- and output rate. The magic of polyphase filtering allows to operate at a lower rate.

In fact, I'm just trying out the rational resampler GNU Radio comes with:

and on my slightly dusty AMD Phenom PC, it achieves 14 to 20 MS/s, depending on the fractional BW I specify (0.1 to 0.48).

However, I'm on your side with the FFT approach – only that you already have an OFDM system and hence, your OFDM modulator can do the job for you.

Now, 5 MS/s would be 333.333… carriers of spacing 15 kHz. So, if we can't live with 0.1% rate mismatch at the transmitter, there's no trivial way to hit exactly 5 MHz at a carrier spacing of 15 kHz.

The solution is easy:

Generate a 15 MHz wide OFDM signal (with still only 200 carriers occupied) by using a 1000-FFT. (no, I hear that regularly, it's not necessary that FFTs have power-of-two lengths. Only small prime factors.) That'll be roughly an increase in computational power needed for the FFT in the OFDM transmitter, but it should be in an order below 10, and in my experience, the FFT isn't even the worst part of a software-based OFDM transmitter (channel coding is).

I happen to have made an approximate graph of how much energy is in sidelobes of sincs over full subcarrier distances for another question:

If you can live with all the energy between the 56. and the ca 700. place in that graphic aliasing over into your used carriers, you don't need to filter at all when downsampling by a factor of 3.

You probably would still want to do some filtering, but this graph means that you can use a lot less great filters. Basically, anything with roughly a 1 MHz transition width and maybe 23 to 36 dB out-of-band attenuation (go for minimum OOB energy, ie. a Dolph-Chebychew window-based design, if you want to really optimize). That's not going to be a long filter, and hence not put much load on your system. In my test, 11 taps might suffice, and 31 looks pretty good.

Even 31 taps can be handled pretty well with FFT filtering at 5 MHz output rate (again, polyphase maths allows us to run the decimation filter at the target rate rather than the source rate).


Auto-correlation of the sum of two generic signals

This is correct. You can use $R_{xy}[l]=R_{yx}^{*} [l]$ to simplify a little further to

$$R_{ss}[l] = R_{xx}[l]+R_{yy}[l]+2 \cdot  \Re\{ R_{xy}[l] \} $$


Is information lost due to Hanning window perturbations?

The STFT itself is invertible, which means by means of the STFT no information about the original signal is lost.

For common implementations of the discrete STFT, we do not calculate the FFT of each possible shift of the window (i.e. every sample), but instead we have a bigger step size (of a few samples in time domain). This sampled STFT (because we evaluate the STFT(t,f) only at t=k*T, where T is the spacing between the windows) is actually a discrete Gabor Transform (of a 1D signal, most references online refer to Gabor Transform of images). In superficial terms, the Gabor transform is invertible as long as the sampling of the STFT is not too coarse. A necessary condition for invertablity (in continuous time) is that $TF\leq 1$, where $T$ is the sampling distance in time and $F$ is the sampling distance in frequency (since in continuous-time, we have infinite frequency resolution).  The invertability becomes the easier, the smaller $TF$, i.e. the finer the sampling of the STFT is done.

So, to answer your question: Despite sounding intuitive that some portion of the signal is lost, if the window is positioned at "unlucky" positions, the STFT does not lose information about the original signal, as long as it is sampled fine enough. So, the signal itself does not determine, if information is lost, but the sampling of the STFT does.


How to remove heart beats interference from pectoralis major electromyogram?

One possible solution is to use a notch filter to remove frequencies commonly associated with heartbeats somewhere between 40 to 100 bpm or .66-1.66Hz

Or if you dont need any low frequencies just use a "highpass" filter where the cutoff is around 2Hz

These solutions assume that anything near the frequency of a heartbeat is unimportant

edit

using the sample data you provided I applied a simple HPF with a cutoff at 2Hz (technically it should be a little past 2, but you will see even as it is it reduced the low frequencies by several orders of magnitude). This may have removed a significant amount of the signal power but it didn't perfectly replicate the cardiac noise. I suppose you have to think the electrical pulse associated with contracting the separate chambers of the heart are much more complex than the resultant heartbeat. In any even here is the data and code, just so you can see the process

x = importdata('pectoralis Major sample.txt');
samples_per_second = 2000;
Ts = 1 / samples_per_second;

%filter design
cuttoff_hz = 2;
cutoff_freq = 2*pi*cuttoff_hz;
tau = 1 / (2*pi*cutoff_freq);
a = Ts / tau;

length_x = length(x);
nfft = 2^nextpow2(length_x);
X=fft(x,nfft)/length_x;

%single pole hpf
x_hpf = filter([1-a a-1],[1 a-1], x);
X_HPF = fft(x_hpf,nfft)/length_x;

figure(1);
subplot(2,2,1);plot(2*abs(X(1:nfft/2+1)));title('original fft')
subplot(2,2,2);plot(2*abs(X(1:nfft/2+1)));title('original fft')
xlim([0 5])
subplot(2,2,3);plot(2*abs(X_HPF(1:nfft/2+1)));title('filtered fft 0-5Hz')
subplot(2,2,4);plot(2*abs(X_HPF(1:nfft/2+1)));title('filtered fft 0-5Hz')
xlim([0 5])

figure(2);
subplot(2,1,1);plot(x);title('original samples')
subplot(2,1,2);plot(x_hpf);title('filtered samples')


Understanding Matlab's Automatic Gain Control function

From the diagram in the Algorithms section of the documentation you can see how the different quantities are computed:

Note that $z$ in the diagram is an estimate of the output power.$^1$ The error signal $e$ is computed by comparing the reference value $A$ to $\ln(z)$. So if you choose $$A=\ln(P)\tag{1}$$ then the average output power will be adjusted to the specified value $P$.

Concerning the maximum gain, note that the input is multiplied by $\tilde{g}=e^g$. If $G$ is the maximum gain in dB, you have the following relationship:

$$G=20\log_{10}(\tilde{g})\tag{2}$$

i.e.,

$$\tilde{g}=10^{G/20}\tag{3}$$

and, consequently,

$$g=\ln(\tilde{g})=G\cdot \frac{\ln(10)}{20}\approx G\cdot  0.11513\tag{4}$$

1. There is an error in the diagram: the output of the detector must be multiplied by the square of the gain (because it's an estimate of the power). This is correctly represented by the equations below the diagram in the documentation.


The key concept that you are missing is that you are not just minimising the difference between input and output signals. The error is often calculated from a 2nd input. Just look at the Wikipedia example related to the ECG.

The filter coefficients in this example are recalculated to change the notch frequency of a notch filter according to the frequency extracted from the mains signal. One could use a static notch filter, but you would have to reject a wider range of frequencies to accommodate the variability in the mains frequency. The adaptive filter follows the mains frequency and so the stop band can be much more narrow, thus retaining more of the useful ECG information.

EDIT:

I have looked at this again and I think I understand your question a little better. The LMS algorithm needs an error term in order to update the filter coefficients. In the ECG example that I paraphrase above, I give the error term as a second input from a mains voltage. Now I'm guessing that you are thinking, "Why not just subtract the noise from the signal-plus-noise to leave the signal?" This would work fine in a simple linear system. Even worse, most examples given online tell you (correctly but confusingly) that the error term is calculated from the difference between the desired signal and the output of the adaptive filter. This leaves any reasonable person thinking "If you already have the desired signal, why bother doing any of this!?". This can leave the reader lacking motivation to read and comprehend the mathematical descriptions of adaptive filters. However, the key is in section 18.4 of Digital Signal Processing Handbook, Ed. Vijay K. Madisetti  and  Douglas B. William.

where:

x=input signal,
y=output from filter,
W=the filter coefficients,
d=desired output,
e=error

  In practice, the quantity of interest is not always d. Our desire may
  be to represent in y a certain component of d that is contained in x,
  or it may be to isolate a component of d within the error e that is
  not contained in x. Alternatively, we may be solely interested in the
  values of the parameters in W and have no concern about x, y, or d
  themselves. Practical examples of each of these scenarios are provided
  later in this chapter.

  There are situations in which d is not available at all times. In such
  situations, adaptation typically occurs only when d is available. When
  d is unavailable, we typically use our most-recent parameter estimates
  to compute y in an attempt to estimate the desired response signal d.

  There are real-world situations in which d is never available. In
  such cases, one can use additional information about the
  characteristics of a “hypothetical” d, such as its predicted
  statistical behavior or amplitude characteristics, to form suitable
  estimates of d from the signals available to the adaptive ﬁlter. Such
  methods are collectively called blind adaptation algorithms. The fact
  that such schemes even work is a tribute both to the ingenuity of the
  developers of the algorithms and to the technological maturity of the
  adaptive ﬁltering ﬁeld

I will keep building on this answer when I get the time, in an attempt to improve the ECG example.

I found this set of lecture notes to be particularly good also: Advanced Signal Processing Adaptive Estimation and Adaptive Filters - Danilo Mandic


Generally, one wouldn't gain much from such parameter tuning or workarounds.
The accuracy heavily depends on the image features that you are using, together with KNN algorithm. If you are using an approximate ANN, try using FLANN. In FLANN, increasing the number of trees and the iterations contribute to the accuracy.

To improve the matching accuracy, enhancements over descriptors are also possible. Zisserman proposes RootSift, which is very easy to implement and enhances SIFT with an intuitive modification:

http://www.robots.ox.ac.uk/~vgg/publications/2012/Arandjelovic12/arandjelovic12.pdf

There are also other good details in this paper.

Other than that, if you are already using BoW, you end up with sparse features. For a more meaningful representation, try VLAD (http://www.vlfeat.org/api/vlad.html). Yet, a quantization through a vocabulary tree as in Nister and Stewinius (http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.61.9520) would help.

Imposing spatial neighborhoods isn't straightforward, but can be done. However, I recommend you to first try the aforementioned methods with post processing (geometric verification). A nice one is found here: http://web.stanford.edu/~sstsai/selected_papers/2010_GeometricReranking.pdf

Cheers,


Background Removal and Object Tracking - Track Fish Jumping Out of the Water

I`d recommend you to go through any of this courses (they include lectures, references and task for labs)

https://cw.felk.cvut.cz/wiki/courses/ae4m33mpv/start - course from Czech Technical University
http://cs.brown.edu/courses/cs143/ from Brown University


How to choose a microphone and make sure not to overcharge (clip) it?

Choose a commercial Sound Pressure Level meter to roughly measure the noise level close to the machine
Look at the data sheet of the MEMS and see if the specification matches the expected level with enough head room (at least 10 dB). IF not, chose a MEMS that meets this spec
Maker sure the MEMS is mechanically de-coupled from the machine, i.e. picks up only sound and NOT vibration.


Rough estimate of how close a human speaker is to a microphone?

Here is one possible suggestion: use a second order microphone (e.g. cardioid microphone). The human voice is a point source and has strong reactive velocity near-field component which velocity mics are sensitive too. In plain English this means you get substantially more relative bass when the speaker is close by. You could leverage to compare the energy in the 500-1kHz Band to the energy in the 100-200Hz band. With some reasonable calibration you ought be to able to tell whether the speaker is close or far.


FM: Do rational ratios of carrier to modulator signals really produce harmonic spectra?

That paper (1973) is about more natural, simple and rich audio synthesis using FM (which later became Yamaha's very populer FM synthesis audio chips).

In that respect, the carrier frequency is taken very low compared to the modulating tone frequency, unlike in communication applications. Hence the the side lobes cross the zero frequency boundary and are reflected back into the positive band using simple trigonometry.

Another alternative to see this is to use positive and negative frequency version (exponential) of Fourier series analysis, to see that negative frequency side lobes penetrate into the positive (vice versa) band when the carrier frequency is insufficiently low.

Note that the resulting Bessel spectra of the single tone modulated FM signal predicts that the side lobes about the carrier will be placed at the frequencies of:

$$f_n = f_c \pm n f_m $$
where $f_c$ is th carrier frqeuency, and $f_m$ is the single tone modulator. As you can see, if $f_c$ is sufficiently low, then for some $n$, the component at $f_n = f_c - n f_m$ will be less than $0$ Hertz and will be reflected into the positive frequency $f_r = n f_m - f_c$.

So the FM side lobes that reflect from negative frequencies and the already existing ones do mix up to create the final audio spectrum which can be harmonic or inharmonic depending on the carrier frequency and modulating tone frequency ratios.

This was a genuine way of producing harmonically rich audio synthesis in legacy analog synthesizers.


White noise and its phase after narrowband bandpass filtering

...And now for a differing opinion....

The OP's representation of bandpass white noise as
$$n(t) = n_I \cos(2\pi f_ct) - n_Q \sin(2\pi f_ct)\tag{1}$$ is inadequate; because each sample path of this noise process is a pure sinusoid of fixed frequency $f_c$ Hz which is not noise-like at all. Why so? Well, a sample path is what one gets when all the random variables are assumed to have taken on appropriate values, and so with $n_I = \frac{\sqrt{3}}{2}$ and $n_Q = \frac 12$ (say), we get the sample path
$$\frac{\sqrt{3}}{2}\cos(2\pi f_ct) - \frac 12 \sin(2\pi f_ct) = \cos\left(2\pi f_ct+\frac{\pi}{6}\right)$$ which is a pure sinusoid as claimed.  A better representation of $n(t)$ is
$$n(t) = n_I(t) \cos(2\pi f_ct) - n_Q(t) \sin(2\pi f_ct)\tag{2}$$
where $\{n_I(t)\}$ and $\{n_Q(t)\}$ are low-pass uncorrelated white noise processes with identical power spectral densities that have constant value for $|f| \leq \frac{B_T}{2}$. Note that a typical sample path from $(2)$ is "sort of" sinusoidal-looking at frequency approximately $f_c$ Hz,  but the instantaneous amplitude and instantaneous phase (as well as instantaneous frequency (derivative of instantaneous phase) are "slowly varying" functions of time. Here, "slowly varying" means in comparison to the approximate time period $f_c^{-1}$ of each cycle of the sinusoid.

All this is fine and dandy but has nothing to do with the question being considered: the distribution of the instantaneous phase $\operatorname{atan2}(n_Q(t),n_I(t))$. The WSS theory tells us nothing about the distributions of the random variables $n_Q(t)$ and $n_I(t)$, except, of course, when $n(t)$ is a Gaussian process in which case $\{n_I(t)\}$ and $\{n_Q(t)\}$ are independent Gaussian processes and we know for sure that $\operatorname{atan2}(n_Q(t),n_I(t))$ is indeed uniformly distributed on $[0,2\pi)$. But, I contend that since we are discussing white noise processes, any distribution other than $\mathcal U(0,2\pi)$ for the phase would mean that the noise exhibits a preference for some phases and a dislike for some other phases, and this idea sticks in my craw as being totally unrepresentative of white noise. Indeed, even the answer by @StanleyPawlukiewicz that $\arctan(n_Q(t)/n_I(t))$ is uniformly distributed on $[-\pi/2,\pi/2]$ is disturbing because it suggests that it might be the case that $n_I(t)$ is always positive while $n_Q(t)$ has a more symmetric distribution. We really need $\operatorname{atan2}(n_Q(t),n_I(t))$ for defining the phase; $\arctan$ is inadequate for the job.

So, there you have it, folks, The phase is uniformly distributed on $[0,2\pi)$ because it offends my sensibilities to have any other distribution for the phase, and that's all there is to it; ymmv.


Adding Salt n Pepper noise into an Image manually

It seems that the final image is in the variable "b". Observe that the max (salt) and min (pepper) values are respectively 1 and 0. This indicates that your original image needs to be an intensity image with graylevels normalized to [0,1]. This may be the reason why you don't see a sensible result when you display "b" as an image.

In case your image is grayscale {0,...,255} i.e. 8-bit unsigned integer format (unit8), you can change your code to assign 0 and 255 instead of 0 and 1.

I suggest you read the documentation on various image storage types here: http://www.mathworks.com/help/toolbox/images/f14-13543.html

It can get pretty confusing at times, mainly because there are so many names (intensity image, grayscale image, RGB, HSV, binary image, indexed image, etc.)


In its bare classical form the RLS algorithm, recursively (for every new iteration), solves the classical problem of least squares; by computing the optimal FIR transversal filter coefficients $w[n]$ which minimizes the sum of error squares between the produced filter output $y[n]$ and a desired response $d[n]$ where the error sum is based on the entire data set beginning from the first sample $y[0]$ to the last sample (the current sample) $y[n]$. Therefore as time progresses, $n \to \infty$, the range of data set also goes to infinity, which means that the memory of the algorithm is infinite; i.e., the current value of the output $y[n]$ and the FIR coefficients $w[n]$ depend (indirectly via recursion) on the entire sample set.

A consequence of the infinite memory is that the estmation becomes more and more exact and refined given that the inputs were WSS random process. On the other hand a clear drawback of this infinite memory dependency is the fact that if the inputs are not WSS then the filter response to reflect the new conditions is slow (the entire past which is irrelevant to the current computation still effects the output due to infinite memory) and tracking is lost.

In order to improve this last situation, and therefore to make the algorithm track non stationary inputs, its memory should be shortened; i.e., only a finite amount of last local samples should effect the current computations. This is achieved by a forgetting factor $\lambda$ embedded into the recursion equations. The range of $\lambda$ is between $0$ and $1$ where $\lambda = 1$ converts the algorithm to classical form with infinite memory with no tracking and $\lambda \to 0$ has less and less memory and more responsive with better tracking...


none of them are showing WHAT TO DO AFTER I PERFORM DTW!.

There are two quotes, at least in the page that you have linked that hint at what you do next, these are:

"We assume that you are familiar with the algorithm and focus on the application. Further information about the algorithm can be found in the literature..."

The linked reference is not a bad start for a book on audio signal processing. The second quote is towards the end:

"...you could apply time scale modification algorithms, e.g. speed up the slower signal to the tempo of the faster one."

And this provides a way towards the answer to your secondary question too.

The "cost matrix" that Dynamic Time Warping returns is discrete with units of samples (e.g. match sample-to-sample), or time.

With reference to this diagram from the Librosa documentation:

You can see that this "picture" has two axes: The $x$ and $y$ axis. On $x$ axis, are the time instances that samples occur on one waveform. On the $y$ axis, are the time instances that the same samples occur on the second waveform.

If the two waveforms were identical, their samples would be played at the same instances and the cost matrix would describe a straight line which means that sample 0 on $x$ is found on sample 0 on $y$, sample 1 on $x$ is found on 1 in $y$ and so on until sample $N$.

But, most of the times you apply DTW, you are interested in aligning the two waveforms.

In that case, the "cost matrix" can be interpreted as a speed-up or slow-down factor, so that the samples of waveform $x$ catch-up with those of waveform $y$.

  ...why interpolation is done at all in audio?

When you think of the cost matrix as modifying the speed by which samples are reproduced, you start thinking about how to speed up or slow down a waveform.

Especially in audio, there are many different methods by which you can achieve this. But the easiest way to think about it is as if changing the speed by which your "sample pointer" rolls through the waveform.

The most trivial example is this: At the waveform's $Fs$ (Sampling frequency), the sample pointer rolls through the waveform at a rate of $Fs$ samples per second. To make the recording sound twice as fast, you can replay it at $2 \cdot Fs$. But, what happens when, to make the two waveforms catch up with each other, you have to run the "record" at $1.056 \cdot Fs$?

When the waveform is played back at normal speed, the sample sequence is $0, 1, 2, 3, 4, 5, 6, \ldots$. When the waveform is played back at twice its normal speed, the sample sequence is $0, 2, 4, 6, 8, 10, \ldots$.

But for non-integer factors, you would have to "land" in between two samples. For example $0, 1.056, 2.112, 3.168, \ldots$.

This is the reason for using interpolation. For example, to work out the sample value at time $1.056$ when you know the sample values at $1$ and $2$, you could use plain simple linear interpolation. This assumes, that when the Analog-to-Digital Converted (ADC) is not sampling the signal, this signal evolves smoothly (as a line) between its values. For high enough $Fs$, this assumption might be valid.

So, suppose that your signal is in $x[n]$ and you want to work out the signal value somehwere (let's call this $k$) in between $n .. n+1$. So, $n < k < n+1$. Then, $v = x[n] + (k-n) \cdot (x[n+1] - x[n])$.

This method is implemented in numpy 's interp() function or even MATLAB's interp1().

Because of this need of having to read in between samples, interpolation is used in such operations as reverberation, chorus, pitch shifting and others. This linear interpolation method will work in general for high enough $Fs$ but as the speed-up/slow-down factor becomes bigger and potentially, applied over a longer block of audio, the artifacts become more audible. This is what more complex methods of applying pitch shifting (such as the phase vocoder), are trying to address.

Hope this helps.


ARMA & MA methods: how do you know the error terms?

In the context of discrete-time statistical signal processing, an ARMA-(p,q) random process is defined as (assuming zero initial conditions)

$$ \sum_{k=0}^p a_k ~ x[n-k]  = \sum_{k=0}^q b_k ~ v[n-k] $$

or equivalently

$$ x[n] = - \sum_{i=1}^p a_i ~ x[n-i]  + \sum_{i=0}^q b_i ~ v[n-i] $$

where $v[n]$ is a white-noise (WSS) random process with variance $\sigma_v^2$ and $x[n]$ represents the resulting ARMA process.

In this context your coefficients are related as $$\phi_i = - a_i $$ and $$\theta_i = b_i$$ The process reggresses over its past values $x[n-i]$ ($X_{t-i}$ acc) and is also a moving average of the input noise $v[n-i]$ ($\epsilon_{t-i}$ acc).


You should take the FFT on Is+j*Qs. If you have say 1000 IQ samples, you can form 1000 I+jQ complex values and take 1000 point FFT to get another set of 1000 complex numbers. If you plot magnitude of these output values, it will show you the frequency content of those 1000 IQ samples.

Since you took 1000 point FFT, each of those complex output values correspond to frequencies at multiples of $\frac{F_s}{1000}\text{Hz}$ where $F_s$ is the Sampling Frequency. For example, the first value will be corresponding to frequency at $0\text{Hz}$, next value is the frequency content at $\frac{F_s}{1000}\text{Hz}$, 3rd value corresponds to frequency $\frac{2F_s}{1000}\text{Hz}$ and so on..

If you are using numpy fft, the documentation should be a good enough guide to you.


What does Fast and Slow Time Weighting signify for real-time audio signal processing?

By "time weighting" in acoustics we are applying an RC (resistor–capacitor) circuit to a time signal.

See: RC circuit on Wikipedia

The time constant for FAST is 0,125 s and 1 s for SLOW.

You can also apply this digitally on your time signal by using my code below - it is basically a exponential moving average filter. I think you understand the steps by just looking at the comments in the code below.

  public class FilterTimeWeighting {

private int nBands = 29; //number of bands
private float initValue;

FilterTimeWeighting() {
    initValue=0;
}

/**
 * Apply the LAmax RC-filter on the signal. Valid for fsamp=44100 Hz.
 * Returns the signal in Volts after applying filter.
 *
 * //y(n)=a*y(n-1)+(1-a)*x(n)
 *
 *
 *
 * @param signalIn
 * @return
 */
public float[] applyLmaxFilter(float[] signal) {

    int n_samples = signal.length;

    float[] x_signal = new float[n_samples];
    float[] y_signal = new float[n_samples];
    float[] ret_signal = new float[n_samples];

    //Square the signal
    for (int i_sample = 0; i_sample < n_samples; i_sample++) {
        x_signal[i_sample] = signal[i_sample] * signal[i_sample];
    }

    //Exponential moving average parameters
    float tau=0.125f;
    float fs=44100f;
    float alpha=tau/(1f/fs+tau);

    //for the first sample
    y_signal[0] = initValue*alpha+(1-alpha)*x_signal[0];

    //for the rest of the samples
    for (int i_sample = 1; i_sample < n_samples; i_sample++) {
        y_signal[i_sample] = y_signal[i_sample-1]*alpha+(1-alpha)*x_signal[i_sample];
    }

    //remember last sample for next time step;
    initValue=y_signal[n_samples-1];

    //Square root
    for (int i_sample = 0; i_sample < n_samples; i_sample++) {
        ret_signal[i_sample] = (float) Math.sqrt(y_signal[i_sample]);
    }

    return ret_signal;
}

  }


I highly recommend taking a closer look on great book:

  Zolzer U. - DAFX Digital Audio Effects

As well as his publication on this specific topic:

  Disch S., Zolzer U. - Modulation and Delay Line Based Digital Audio
  Effects

Leslie (Rotary) speaker can be simulated using tremolo and Doppler effect for both change in signal amplitude and pitch.

Basic signal flow diagram that can be found in the article:


Frequency Domain Filtering with big kernel size

You dont have sigma specified for the gaussian filter. I typically use something like

fspecial('gaussian',6*sigma,sigma);

as the gaussian kernel is roughly 6 x sigma wide. (see graph here but note they plot against variance https://en.wikipedia.org/wiki/Gaussian_function)

Because of that, your padding should be 3 x sigma. But padding with zeros creates an edge in the image that might not have been there. You may get better results by replicating the outer pixels instead of padding with zeros:

padarray(I,[3*sigma,3*sigma],'symmetric');

Also you might want to do an ifftshift before doing the fft2 of the filter (note the i)

You could up this to 8 x sigma (and 4 x sigma) if you really want most of the kernel. The gaussian kernel has infinite support so there will always be some circular convolution effect - you can't pad to infinity - but with 3 or 4 sigma padding it will be small.

btw MATLAB recommends using imgaussfilt http://au.mathworks.com/help/images/ref/imgaussfilt.html


You should position the board at roughly the same distance from the camera as the distance between the camera and your objects of interest. The board should be big enough to cover a good portion of the field of view at that distance.

Another point: ideally, your board should have an even number of squares along one side, and an odd number of squares along the other.  Then the orientation of the board can be detected unambiguously.

Finally, since you have a "matlab" tag on your question, please check out the Camera Calibrator app, which is a part of the latest release of the Computer Vision System Toolbox for MATLAB.


a neural network approach for FIR filter

Yes. The FIR filter model you're used to is a series of Neurons with weighted inputs, and a linear activation function.
In other words, a standard FIR filter is a neural network.
I mean, it's called "CNN". The C is exactly the operation a filter does.


How to determine precision needed for sin approximation used for sound synthesis?

This is very good question.
Fortunately for you, I have a very good answer.
Assuming standard CD quality, your sound level has just over 4 significant digits accuracy, thus this is the level you need to be truly indistiguishable.
Is this level necessary for auditory purposes?  Let's assume so.
Your fastest solution, by far, is going to be linear interpolation on a segmented domain.  You need two tables.  The first hold the sine value of the center of each interval.  The second holds the first derivative, aka cosine.
To calculate your sine value, multiply your angle by the spacing you are using (so a unit of one corresponds a segment).  The integer portion will give you the index into your table and use the fractional portion (-0.5 to 0.5) for the interpolation.
Now the question becomes "How many entries do I need in my table?"
I'll let you take a crack at that.

I have a better answer (I think, you need to test it).
$$ \sin( x + d ) = \sin( x ) \cos( d ) + \cos( x ) \sin( d ) $$
Back to two tables.  One for the broad range ($x$) , and one for the fine range ($d$).
If the range of $d$ is small enough, then you can use (Taylor or find the Remez instead):
$$ \cos(d) \approx 1 - x^2/2 + x^4/24 $$
$$ \sin(d) \approx x - x^3/6 + x^5/120 $$
But that will require more computation.
For uber accuracy, which you don't need, you could do the interpolation thing on the fine table value.

I was curious, so here you go:

import numpy as np

#==========================================================
def main():

        N_coarse = 128
        N_fine   = 128

#---- The coarse table

        sc = np.zeros( N_coarse )   # Sine Coarse
        cc = np.zeros( N_coarse )   # Cosine Coarse

        theSlice = np.pi * 0.5 / N_coarse
        theAngle = 0.0

        for n in range( N_coarse ):
          sc[n] = np.sin( theAngle )
          cc[n] = np.cos( theAngle )
          theAngle += theSlice

#---- The fine table

        sf = np.zeros( N_fine )   # Sine Fine
        cf = np.zeros( N_fine )   # Cosine Fine

        theSlice /= N_fine
        theAngle = 0.0

        for n in range( N_fine ):
          sf[n] = np.sin( theAngle )
          cf[n] = np.cos( theAngle )
          theAngle += theSlice

#---- The test

        theFactor = N_coarse * 2.0 / np.pi

        for a in range( 157 ):
          theAngle = a * 0.01
          n = theAngle * theFactor

          nc = np.floor( n )
          nf = np.floor( ( n - nc ) * N_fine )

          sine = sc[nc] * cf[nf] + cc[nc] * sf[nf]

          print sine, np.sin( theAngle ), sine - np.sin( theAngle )

#==========================================================
main()

This old fart got it right on the first try!
Here are the first few rows of the output, the rest is comparable:

0.0 0.0 0.0
0.00997070990742 0.00999983333417 -2.91234267487e-05
0.0199404285515 0.0199986666933 -5.82381418187e-05
0.0299081647675 0.0299955002025 -8.73354349791e-05
0.0399687249608 0.0399893341866 -2.06092258625e-05
0.0499294807897 0.0499791692707 -4.96884810173e-05
0.059885272753 0.0599640064794 -7.87337263996e-05
0.0699307504776 0.0699428473375 -1.20968599355e-05
0.0798735744039 0.0799146939692 -4.11195652507e-05
0.089808457497 0.089878549198 -7.00917010057e-05
0.0998298073783 0.0998334166468 -3.60926852212e-06
0.109745746461 0.109778300837 -3.25543758203e-05
0.119650774894 0.119712207289 -6.14323952304e-05
0.129543907942 0.12963414262 -9.02346778173e-05


Be careful, a median filter cannot be expressed as a convolution, and thus is not considered a kernel in this respect. This is because the median filter is based on order statistics of an image patch, and the resulting pixel at the output of a median filter is not a linear combination of other pixels within a patch.

Otherwise, you are right, kernels are generally thought of as just other 'small images', that one would convolve with the image to be analyzed. Those kernels can be chosen for a variety of reasons. Gaussian kernels are used for smoothing, laplacians are used for edge detections, etc.

There are as many kernels as you can imagine. However usually we start with a particular application in mind, "I need to blur the image", or "I need to detect edges in the x-direction", and then find kernels that magnify those features we are after.


How do I perform a time domain phase shift in the frequency domain?

Let’s be clear on what we will refer to as time delay and phase shift. Due to the common association of individual frequencies as sinusoids many confuse delay and phase shift as being equivalent. However a phase shift in time is simply multiplying a time sample by $e^{j\phi}$ while a time delay is displacing the sample in time such as done with $x(t-\tau)$.

The same thing occurs in frequency, where a phase shift of a particular carrier is done by multiplying that carrier by $e^{j\phi}$, while a frequency displacement is done by moving the carrier to a new frequency such as $X(\omega-\omega_c)$.

As you showed, a delay in time is a linear phase versus frequency in the frequency domain.

Similarly due to the bidirectional properties of the Fourier Transform, a displacement in the frequency domain would be a linear phase versus time in the time domain, specifically $e^{j\omega t}$. Multiplying a time domain signal by $e^{j\omega t}$ with $\omega t$ as a linear phase ramp, results in a frequency translation, while multiplying a time domain signal by $e^{j\phi}$ with $\phi$ as a constant, results in a static phase shift (or if it makes it clearer: a phase displacement, or a phase rotation). This phase displacement is not directly a time displacement.

Further, any multiplication in the frequency domain would be equivalent to a convolution in the time domain (and vice versa). If you want to apply a phase change in the time domain; if the intention is equivalent to a multiplication of the time domain samples- this would need to be done by a convolution effectively in the frequency domain.

To provide further details we need to understand more precisely what is meant or intended by a “phase shift in the time domain”, specifically if this means all time samples are rotated by the same phase (which is done by convolving the frequency domain function with a complex impulse at f=0 with magnitude = 1 and associated phase such as $\delta(\omega)e^{j\phi}$ or if the time samples have a linearly increasing or decreasing phase versus time (which is done by convolving the frequency domain function with an impulse at a frequency offset $\delta(\omega-\omega_c$).


Frequency response given poles and zeros

The MathWorks Documentation of the function zp2tf says that

  The zeros and poles must be real or come in complex conjugate pairs.

which is not the case for your example. With the given values of the vectors z and p you can do the following:

num = poly(z);
den = poly(p);

which gives

1.00000 + 0.00000i  -0.37545 - 0.77876i   0.09714 - 0.02978i

for num, and

1.00000 + 0.00000i   0.83665 + 0.07610i   1.00092 + 0.10802i

for den. Note that these are different from the coefficients of your 'expected transfer function'. However, these coefficients correspond to the given pole and zero locations.

Also note that one of your poles lies outside the unit circle, so your filter is unstable.


What does the stability margin have to do with the quality of a filter?

Even though you might have designed your filter to be stable, Stability Margin is a concept which is used to estimate how far we are away from becoming Marginally stable/unstable both in terms of Gain and Phase.

This would mean then at Unity Gain(0 dB) any phase rotation of -180 degrees or ODD multiples of it is bound to make your system unstable and thus comes Phase Margin, the amount of rotation on the Pole Zero plot which would bring the Root Locus plot to intersect the point -1(or rather how far you're away from the 180 degree line on a Bode Plot when your Gain Magnitude is Unity).
Likewise we can state the same in terms of Gain. The amount of Gain that'd be necessary for the same Nyquist plot to expand and hit the point -1 (or how far you're away from the Unity magnitude gain when the phase response is -180 degrees).

Practical consequences of system implementation always require you to have some Margin because you're never guaranteed that the Gain and the Phase plots are exactly replicated when the system is actually working. So these parameters will tell you if you're safe or not. Theoretically, even being away from 1 degree away from -180 Degree line at Unity gain is stable; but practically there are things that can happen which will push it away easily.

You can learn a lot about these at

:http://en.wikipedia.org/wiki/Phase_margin
http://www.facstaff.bucknell.edu/mastascu/econtrolhtml/Freq/Nyquist4.html (a good simulation is shown here).

Stability Margin doesn't have any impact on the Stability of the Filter, it'll just indicate if your filter is Stable or not.

The Phase plot is actually shown in Degrees. Your filter looks like it is Marginally stable because at a few multiples of -180 degrees, the Gain is 0dB(oscillatory), and that too within the frequencies you've designed the filter to operate, it looks not promising at all. The Black Lines on the gain plot should appear at every multiples of -180 degrees. Also coming to the phase plot, there are no Black lines in this case because at least one Unity Gain crossover point has a Phase rotation that is an Integral odd multiple of -180 Degrees.


Is linear convolution same thing as aperiodic convolution?

What they call noncircular or aperiodic convolution is more commonly called linear convolution. Linear convolution describes the input-output relation of linear time-invariant (LTI) systems. Linear convolution can be computed via circular convolution if we use sufficient zero-padding. Circular convolution is important because it can be computed using fast algorithms (FFT).


There's a lot of different domains of knowledge coming together here, so I'll split my answer into multiple sections, each answering an implicit question that you raise in your explicit question. Hope that helps!

Can your RTL-Dongle actually receive at 1.72 GHz?

So, first the bitter pill: There's a lot of sellers out there that offer RTL dongles and claim that they've got an R820T tuner, but actually ship an Fitipower tuner that covers far less spectrum. I find that the most likely explanation here, especially considering the product description and the "tech specs" table disagree. Said table only claims:

  Receiving Frequency:48.25-863.25 MHz Tuning Range

So, possibly you've got an RTL-Dongle with a fitipower FC0012, not a R820T (which went out of production a while ago). RTL-SDR comes with the rtl_test utility, but possibly your software can show you the tuner installed in that device, too.

However, considering your display shows the right frequencies, let's assume it actually tuned to ~1.720 GHz when you instructed the device to tune there.

Are you using an antenna that's suitable for these frequencies?

Just mentioning that, because I've supported far too many people that told me after a couple of hours that they were trying to receive FM broadcast radio at 100 MHz using a Wifi antenna, and vice versa.

Are the powers in these different applications related?

Then, first thing to notice: The powers displayed in your SDR software are not related to the power your App shows for reasons:

The power shown is the downlink power. You're observing uplink. Uplink should probably have significantly more power than downlink at your handset.
I don't know where the app gets that data. But, I'll go out on a limb and say it's not really an RX power estimate, at all. LTE specifies a whole set of quality metrics, all of which are useful to assess the error probability of a receiver. None of these is a pure power measurement (because that doesn't matter to a device, actually. What matters are SNR measures that incorporate noise, but also things like adjacent channel interference, clock errors, …).
An RTL dongle doesn't even remotely give calibrated values. What it sends to the computer are simply 8 bit numbers, relative to the ADC's full scale, whatever that means at any given gain, frequency, rate, filter setting, temperature, component tolerances, moon phase, amount of black cats in neighborhood… . I don't know why your program shows you "dBm", but unless you have run some calibration on your dongle for exactly your gain, frequency, and sampling rate settings, these numbers mean nothing. All you can say is that 3dB more means twice as much power.

Can I analyze LTE signals using an RTL Dongle, at all?

LTE uses a total bandwidth much larger than 3 MHz, so, if at all, you're looking at a small excerpt of your uplink; possibly even one where your handset doesn't even get an allocation of resources.
LTE is very flexible about assigning resources to different user equipment in the same cell, and very quick – unless you have a very aggressive means of holding spectral peaks, it might be pretty impossible to spot your own contribution in a cell.

Especially the downlink uses OFDM, and usually also power control, so if your handset is the main downlink user in a cell, all that you should be seeing if you had an SDR that could actually see the whole LTE downlink channel, would be what would look like a "flat-top" section of spectrum that is basically elevated white noise, possibly exhibiting non-flatness due to frequency-selectivity of the channel.

If you look at a section of OFDM, all you'll see is a noise floor, unless you manage, by luck or synchronization, use an FFT that exactly "hits" (integer factors) of the subcarrier spacings (in which case you'd be lucky to see staggered sinc spectra of the OFDM subcarriers). You'd need a very long FFT and very good SNR (which you don't seem to have) to see the sinc carriers unsynchronized.

Are you actually observing your phone's LTE signals in your pictures?

Doesn't look like it.

Your 800 MHz spectrum shows something 200 kHz wide with a rectangular shape; that might be a bunch of LTE subcarriers (those have 15 kHz spacing), but it's definitely not coming out of the noise floor far enough to be your phone's resource blocks (also, questionable whether you'd get a resource block allocation that's static for that long).

Considering the periodic spurs and the intereferer at ~1719.35 MHz that looks a lot like a mixed-up rectangular wave, I'll go with no, you're likely not looking at LTE here; there should be far more dynamics in the waterfall.

Is your ultimate goal achievable?

  My goal here is to be able to detect my phone's transmitted signal and then, to make sure I can distinguish it before I invest in a more expensive SDR to be able to monitor higher bands.

Detecting activity of any one phone in vicinity

As said, your TX signal should be much stronger close to your handset than it's RX (which really is just the simple physics of spatial power density going down with the square of distance). Especially when your phone displays "few bars" (i.e. it has low receive SNR), it will have to use a lot of power when transmitting (which is why I'm pretty sure you're not observing your phone in either screenshot).

A simple power detector with a high-enough threshold should do an approximately good enough job. You'd still not be sure whether it's your phone, or the phone in the next room. LTE relies on power control, i.e. your phone will only use as much power as necessary to achieve the rate it needs (that's not something that your phone is allowed to choose – it's necessary to keep the dynamic range of the receiver on the base station side manageable and to keep interference at bay).

Detecting activity of a specific phone

Now, this get's a lot more complicated, as you'll have to analyze the LTE signal. Not quite sure from the top of my head, but I'd go with "no, from the uplink signal alone, you can't say which phone is currently transmitting".

LTE is a cellular infrastructure network, not an ad-hoc burst network; hence, stations are assigned times when they should be transmitting by the base station using the downlink. Inherently, that means that the station doesn't need to indicate who it is when using the uplink; the base station knows because it assigned the uplink resources, who's using the uplink. Thus, you'd have to observe the downlink more than the uplink to know when you're phone's talking. Luckily, LTE is generally pretty well-encrypted, so you shouldn't be trivially able to extract very specific information from the downlink about your phone, but you might at least be getting unencrypted local IDs of the different stations in the cell from the downlink, if your SDR receiver software is able to demodulate the signal.


Key detection Dance Music inaccuracies - matlab

A 12-dimensional chromagram may not estimate pitch very well, given synthesized semi-tonal sounds that include lots of energy at multiple odd harmonics, as those harmonics will often end up in the wrong pitch class slots.  If there's lots of tonal bass, low-pass filtering just above the primary range of the bass instrument/synth's fundamental might help getting rid of some of the harmonics and overtones.

If the drums are not tuned to the key (probably not), then they will interfere with any key centered pitch class information.  Perhaps time-domain filtering out of all percussive events, and only analyzing the tones between the drum beat attacks might help in this case.  Even a simple low-envelope amplitude gate function might be suitable (assuming the sound isn't so compressed that there is no lower amplitude in any extent portion of the envelope!).

And, as per the previous answer, are you sure the sounds of your "dance music" actually contains a musical key? (as opposed to producing some cognitive or auditory illusion of tonality).


Difference of Gaussian filter but using Gaussian CDFs not PDFs

You need to ask yourself why do we use the difference of Gaussians from the first place? The reason is because the difference will give us a measurement for the change in value around the point we apply it to as a function of the variances of the Gaussians. If we have a big change the difference between the Gaussians it means that we have some frequencies that were low passed by the Gaussain with the bigger variance. If the change is low it means that there i no information lost. So we can think of it as a band-pass filter.

A simplified reason why it acts like that is because of the shape of the Gaussian PDF as can be seen in the following figure.

If we would use the CDF of the Gaussains and not the PDF it will give us a positive difference on one side and a negative difference on the other side.


What is hidden in DFT below frequency resolution?

Computing your DFT with a 0.5 Hz spacing is equivalent to appending your signal with 1000 zeros and then computing a 2000 point DFT. It should be clear that this adds no real information to the signal. Yet it can sometimes be useful.
What really happens by this zero padding is that you're seeing an interpolated version of your spectrum. This means that between each two "valid" points (at 1 Hz spacing) you're seeing one more that produces an interpolate, a potentially smoother version of your spectrum. The interpolation kernel is essentially a sinc function (periodified sinc to be accurate).
Depending on what you want to do with the spectrum, such interpolated results can be visually more pleasing and sometimes more clearly show things we would otherwise have missed.
In particular, if you have one harmonic signal and its frequency is not on the DFT grid, your DFT will not show one line but many. This is also known as the picket fence effect. This picture may be confusing and in cases like this, an interpolated version would show us more clearly what is going on and thus be easier to interpret. Also, this allows us to find the actual frequency with a better precision (NB: resolution is not improved by this, but precision is).
Here are two examples. In each example I considered N = 30 for clarity. For the first example we have a frequency of 5 Hz, falling exactly on the DFT grid. The 30-point DFT (depicted in red) gives just one line). An interpolated version is added with the dashed line just for clarity. You can clearly see the sinc-like interpolation kernels and how in the on-grid case we're sampling exactly the zero crossings, which is the reason we're seeing only one line.

In the second example, the frequency was chosen to be off the grid (between 5 and 6 Hz). Now the 30-point DFT shows many lines and it may be hard to see that these stem from just one harmonic if we only have those. Adding the interpolated DFT once more shows a very similar picture to before: We can see that it looks like one harmonic and we can much better spot its maximum.


What should be the correct scaling for PSD calculation using $\tt fft$

There is only one correct way of scaling DFT when calculating PSD with RMS values. Given input signal $x$ and its DFT $X$, the exact formula is:

$$\mathrm{PSD}=\frac{2\cdot \hat{X}}{f_s\cdot S} $$

where:

$\hat{X}=|X|^2=X\cdot X^*$ - squared spectrum magnitude

$S=\sum_{i=1}^{N} w_{i}^{2}$ - scaling factor defined as sum of squared samples of window function

$f_s$ - sampling frequency

Scaling coefficient $2$ takes into account removal of energy at negative frequencies (we drop this side of PSD).

Sampling frequency $f_s$ happens to be in the equation since it is a part of the Effective Noise Bandwidth, and PSD is simply PS (Power Spectrum) divided by ENBW. So in the end,  the unit of PSD is $x^2/\mathrm{Hz}$, where $x$ is the quantity that you want to measure, i.e. Volts.

Scaling factor $S$ is compensating for any loss of energy when applying a window function. Obviously for rectangular window (or no window in fact) this doesn't change anything and one, might simply use $N$ instead (length of input sequence).

You did also mention that you were planning to use high number of DFT points ($2^{20}$) since it "makes more sense". Well, it doesn't. If your input sequence has length $N$ (assume 1024), then you can't get better frequency resolution than:

$$\Delta f = \frac{f_s}{1024}$$

Even if you increase number of DFT points to $2^{20}$ by zero-padding, those spectral samples in between will be only interpolated window spectrum - useless garbage. For more detailed explanation please see my other answer.

Another thing to remember. Zero-padding doesn't change the scaling factor! Samples are added after the input sequence had been multiplied by the window, hence the total energy will not change.

Since you wanted to match the calculation with MATLAB's periodogram function, here comes the code. It is a modified version of this one. I tried to comment the code.

% Parameters and signal
Fs = 1000;
dt = 1/Fs;
t = 0:dt:1-dt;
x = cos(2*pi*100*t)';
N = length(x);
% Add a bit of noise
x = x + 0.001*randn(N, 1);
% Use window function of your choice. Rectangular window means no window.
win = rectwin(N);
% The actual number of points used when calculating the DFT might be higher
% this can help in speeding up the computation (next power of 2), but won't
% increase the frequency resolution (only frequency spacing is affected).
% See: https://dsp.stackexchange.com/a/32161/8202
NFFT = N; % N=2^20
xw = x.*win; % Apply windowing here (periodogram does it on it's own)
xdft = fft(xw, NFFT);
psdx = xdft.*conj(xdft); % psdx = abs(xdft).^2
psdx = psdx(1:NFFT/2+1); % Take only first half of the PSD
% Calculate the scaling factor.
% Note that in case of rectangular window this is simply N.
S = sum(win.^2);
psdx = (2/(Fs*S)) * psdx;
freq = 0:Fs/NFFT:Fs/2;
% Calculate the corresponding PSD using MATLAB's in-built
[pxx, ~] = periodogram(x, rectwin(N), NFFT, Fs);

% Plot the PSD's
plot(freq, 10*log10(psdx), 'r')
grid on
hold on
plot(freq, 10*log10(pxx), 'b--')
title('PSD calculation comparison')
xlabel('Frequency [Hz]')
ylabel('Power/Frequency [dB/Hz]')
legend({'DFT', 'periodogram'})

% Plot error between PSD's
figure
diff = pxx - psdx;
plot(freq, diff)
grid on
xlabel('Frequency [Hz]')
ylabel('Difference')
title('Difference between PSDs')

Please note the perfect match between both PSD's:

And here is the difference:

Probably you can notice that there is a bigger error for 1'st and last bin (Nyquist in fact). That is because we multiplied all frequency bins by scaling factor $2$ to compensate for negative frequencies. This shouldn't be done for DC component (1'st bin) and for Nyquist bin (if it's present) - both of them are unique.

But in a real life, Nyquist bin will be removed by anti-aliasing filter and people will remove mean from the signal anyway. Also DC scaling for RMS values is pointless, especially when doing windowing.


System characterization given pole-zero mapping

I will not do the homework for you, but I'll make a few comments on your observations to clear up some of your misconceptions. Concerning the first plot (top left I suppose) you wrote:

"the system is stable since we have a region of convergence containing the unit circle"

This would be correct if the system were non-causal (i.e., with a two-sided impulse response), but since - according to the first sentence in the question - the system is causal, it can't be stable if it has poles outside the unit circle.

" ... and FIR because it is stable"

This is wrong. What you're saying is that any stable system is FIR, which is definitely not the case. Of course there are stable IIR filters. Here you'll have to do some more studying to learn what stability means and to understand that it is not necessary (only sufficient) for a stable system to be FIR.

"... and minimum phase because we have a pair of complex conjugate pairs"

That's also wrong. The fact that a system has complex conjugate pairs of poles and zeros has nothing to do with the minimum-phase property. You need to read up on minimum-phase systems. They can be inverted by a causal and stable system. What does this tell you about its zeros?
Any system with a real-valued impulse response has either real-valued or complex conjugate pairs of poles and zeros. That's because $H(z)$ is a real-valued function of a complex argument.


What is the purpose of the log when computing the MFCC?

The logarithm serves to transform a multiplication into an addition. It is part of the computation of the cepstrum. The basic idea is as follows:

Assume a source signal $x$ is convolved by some impulse response $h$. The resulting magnitude spectrum is

$$|Y(\omega)| = |X(\omega)||H(\omega)|$$

By applying the logarithm we get

$$\log |Y(\omega)| = \log |X(\omega)| + \log |H(\omega)|\tag{1}$$

If we want to equalize / undo the effect of filtering by $H(\omega)$ we can hope that this task is easier if we transform the convolution into additive noise. This is exactly what happens by taking the logarithm in (1).


2-d circularly symmetric low-pass filter

Excerpted from Jae S.Lim 2D signal and image processing ch.1, as an example of $2$-D circularly symmetric lowpass filter with a cutoff frequency of $\omega_c$ radians per sample, whose impulse response is given by:
$$h[n_1,n_2] = \frac{\omega_c}{2\pi \sqrt{n_1^2 + n_2^2} } J_1 \big( \omega_c \sqrt{n_1^2 + n_2^2} \big) $$

where $J_1$ is the Bessel function of the first kind and the first order...

Interested readers may consult the book for a derivation which is not untitive but nevertheless tractable; familiarity with Bessel functions is required but also provided as is.[Derivation is added.]

[Olli below]

At $n_1 = n_2 = 0$ the limiting value must be used:

$$h[0, 0] = \frac{\omega_c^2}{4\pi}$$

A slice through the middle of $h[n_1,n_2]$ with $\omega_c = \pi$:

Python source for 2-d filter kernel (you may want to apply a 2-d window function):

from scipy import special
import numpy as np

def circularLowpassKernel(omega_c, N):  # omega = cutoff frequency in radians (pi is max), N = horizontal size of the kernel, also its vertical size, must be odd.
  kernel = np.fromfunction(lambda x, y: omega_c*special.j1(omega_c*np.sqrt((x - (N - 1)/2)**2 + (y - (N - 1)/2)**2))/(2*np.pi*np.sqrt((x - (N - 1)/2)**2 + (y - (N - 1)/2)**2)), [N, N])
  kernel[(N - 1)//2, (N - 1)//2] = omega_c**2/(4*np.pi)
  return kernel

Example with $\omega_c = \pi$:

import matplotlib.pyplot as plt

kernelN = 11  # Horizontal size of the kernel, also its vertical size. Must be odd.
omega_c = np.pi  # Cutoff frequency in radians <= pi
kernel = circularLowpassKernel(omega_c, kernelN)
plt.imshow(kernel, vmin=-1, vmax=1, cmap='bwr')
plt.colorbar()
plt.show()

Example with $\omega_c = \pi/4$:

kernelN = 41  # Horizontal size of the kernel, also its vertical size. Must be odd.
omega_c = np.pi/4  # Cutoff frequency in radians <= pi
kernel = circularLowpassKernel(omega_c, kernelN)
plt.imshow(kernel, vmin=-np.max(kernel), vmax=np.max(kernel), cmap='bwr')
plt.colorbar()
plt.show()

[fat32 below]

Ok, for those who would benefit from a (not so intuitive) derivation, here I make a (n almost) verbatim copy of it from the same book.

First, lets simply write the 2D discrete-time inverse Fourier transform to define the impulse response as:

$$ h[n_1, n_2] = \frac{1}{(2\pi)^2} {\int \int}_{\omega_1^2+\omega_2^2< w_c^2} 1 \cdot e^{j(\omega_1 n_1 + \omega_2 n_2)} d\omega_1 d\omega_2 \tag{1} $$

Let's make a change of variables $\omega_1 = r \cos(\theta)$ and $\omega_2 = r \sin(\theta)$ (effectively write the integral in (1) in the cyclindrical (or polar) coordinates) :

$$ h[n_1, n_2] = \frac{1}{(2\pi)^2} \int_{r=0}^{\omega_c} \int_{\theta=a}^{a+2\pi} e^{j r (\cos(\theta) n_1 + \sin(\theta) n_2)} r ~ dr d\theta \tag{2} $$

Now, make a further change of variables $n_1 = n \cos(\phi)$ and $n_2 = n \sin(\phi)$, with $n = \sqrt{ n_1^2 + n_2^2 }$ and obtain

$$ h[n_1, n_2] = \frac{1}{(2\pi)^2} \int_{r=0}^{\omega_c} r ~dr  \int_{\theta=a}^{a+2\pi} e^{j r n \big( \cos(\theta) \cos(\phi) + \sin(\theta) \sin(\phi) \big) } d\theta \tag{3} \\$$

which is now:
$$ h[n_1, n_2] = \frac{1}{(2\pi)^2} \int_{r=0}^{\omega_c} r ~dr  \int_{\theta=a}^{a+2\pi} e^{j r n \cos(\theta -\phi) } d\theta \tag{4} \\$$

defining $f(r) = \int_{\theta=a}^{a+2\pi} e^{ j r n \cos(\theta-\phi) } d\theta $ , then we get:

$$ h[n_1, n_2] = \frac{1}{(2\pi)^2} \int_{r=0}^{\omega_c} r f(r) dr  \tag{5} $$

Now, exploring $f(r)$ we can apply the Euler's identity :

$$ f(r) = \int_{\theta=a}^{a+2\pi} \cos(r.n.\cos(\theta-\phi)) d\theta + j \int_{\theta=a}^{a+2\pi} \sin(r.n.\cos(\theta-\phi) ) d\theta \tag{6} $$

And we can notice that the imaginary part is zero (can check an integral table for that) and with $a =  \phi$  then $f(r)$ becomes

$$f(r) = \int_{\theta=0}^{2\pi} \cos(r.n.\cos(\theta) ) d\theta  \tag{7}$$

Now the integral in (7) is recognized as the the Bessel function of order zero, kind one $J_0(x)$ which is given as:

$$ J_0(x) = \frac{1}{2\pi} \int_{\theta =0}^{2\pi} \cos( x \cos( \theta) ) d\theta \tag{8} $$

from (7) and (8) we see that $f(r) = 2\pi J_0(r n) $...

And the last identity is given as is:
$$ x J_1(x) |_a^b = \int_a^b x J_0(x) dx \tag{9}$$

Now putting the impulse response into the form

$$h[n_1,n_2] = \frac{1}{(2\pi)^2} \int_{r=0}^{\omega_c} r 2\pi J_0( r n) dr \tag{10}$$

applying (9) onto (10) with the recognition that $x = r n$ and $dr = dx/n$ and $n= \sqrt{n_1^2+n_2^2}$ yields the result:

$$\boxed{ h[n_1,n_2] = \frac{\omega_c}{2\pi \sqrt{n_1^2+n_2^2}} J_1( \omega_c \sqrt{n_1^2+n_2^2}) }  \tag{11}$$


so MJ, my approach to using the Gaussian window for the STFT is first to use the "ordinary frequency" definition of the continuous-time Fourier Transform:

$$ \mathscr{F}\Big\{x(t)\Big\} \triangleq X(f) = \int\limits_{-\infty}^{\infty}x(t)\, e^{-j 2 \pi f t} \,\mathrm{d}t$$

$$ \mathscr{F}^{-1}\Big\{X(f)\Big\} \triangleq x(t) = \int\limits_{-\infty}^{\infty}X(f)\, e^{j 2 \pi f t} \,\mathrm{d}f$$

and start with this really cool isomorph of the Fourier Transform:

$$ \mathscr{F}\Big\{e^{-\pi t^2}\Big\} = e^{-\pi f^2} $$

and use the well known time-scaling and translation or frequency-scaling and translation theorems of the Fourier Transform to get you

$$ \mathscr{F} \Big\{ e^{a t^2 + b t + c} \Big\} = e^{A f^2 + B f + C} $$

where the constants $A$, $B$, and $C$ can be explicitly mapped from $a$, $b$, and $c$. It appears to me that the mapping is:

$$\begin{align}
A &= \frac{\pi^2}{a} \\
\\
B &= j \frac{\pi b}{a} \\
\\
C &= c - \frac{b^2}{4a} - \tfrac{1}{2}\log\left(-\frac{a}{\pi}\right) \\
\end{align}$$

and the inverse mapping (which should be self-similar) is:

$$\begin{align}
a &= \frac{\pi^2}{A} \\
\\
b &= -j \frac{\pi B}{A} \\
\\
c &= C - \frac{B^2}{4A} - \tfrac{1}{2}\log\left(-\frac{A}{\pi}\right) \\
\end{align}$$

Looks like $\Re\{a\}<0$ and $\Re\{A\}<0$ for the integrals to converge and for the $\log(\cdot)$ to be real and finite in the mapping.


Audio signal comparison for automatic singing evaluation

Please update your question to mention which kind of comparison you want to perform. You seem to want your comparison to be robust to amplitude changes and time delays, but what about:

Audio coding. Should a signal and a MP3 version of itself be considered as similar or different?
EQ and filtering. Should a signal and the same signal with the high frequencies shifted by 10° be considered similar or different?
Dynamics change. Does applying a time-varying gain factor to one signal make it different from the original?
Time and pitch shifting. Should a signal and a slightly accelerated version of itself be considered as similar or different?
Noise and recording conditions. Should two versions of the same signal recorded on different equipment, microphone type, etc. be considered as similar or different?
Timbre. Should two versions of the same melody played by two different instruments be considered as similar or different?

Depending on your application (assessing degradations through a signal processing chain, identifying copyrighted work in an audio stream, shazam-like song identification, plagiarism detection on music...), you might want your comparison to be robust to some of these transformations.

Your approach consisting in normalizing the signals is not very robust. If you have two identical signals $x$ and $y(t)$, and if you change $x(0) = 2 \max_t y(t)$, the normalized version of $x$ is twice smaller than the normalized version $y$ and the comparison will fail.

Using cross-correlation to look for a time delay is a sound approach.


Simulating musical instruments using their impulse response

there are a lotta people doing research regarding music synthesis and DSP.  yes you can sample notes and different key velocities and play them back with different MIDI velocity parameters.  if you wanted to somehow interpolate between a note played at $mf$ to $f$ and on to $ff$ and $fff$urther, you would need to find a way to phase align the waveforms.  IMO the only decent analytical way to do that in the time domain is wavetable synthesis where the wavetables can be lined up (or it's frequency domain, it's called additive synthesis) - additive synthesis is more general, wavetable requires the overtones of the note to be mostly harmonic.  doing this stuff right is not particularly easy.

on the other hand, since you brought up "impulse response", there are folks who do physical modeling of musical instruments, where they model the instrument as a form of resonant LTI filter and they separate the impulse response of the LTI filter from the driving function, which is not an impulse but something more.  then, for different key velocity hits, they have different driving functions but the same LTI impulse response.  sometimes it's called commutated synthesis or https://ccrma.stanford.edu/~jos/pdf/svd95.pdf.


What papers or resources explain the conditions needed for filter stability in real-time tuning of parameters?

Instead of trying to fit everything in a comment, I will answer the question here.
Assessing stability varies to which method you use to describe a filter (DE, Laplace domain, State - Space, FIR, etc). Stability itself is a very elaborate and interesting topic in controller and filter design, but I must keep it short enough. In essence, if you apply an impulse response to your filter, the energy of the output must either decrease or stay constant to make the filter (marginally) stable. But this does not tell you how this is linked to the parameters. Depending on your filters description, the following things should be checked for stability:

Laplace domain: roots of the denominator $\leq0$
Discrete Z-domain: absolute value of the roots of the denominator $\leq1$
Continuous State Space: Eigenvalues of A $\leq0$
Discrete State Space: absolute value of the Eigenvalues of A $\leq1$

So, suppose you have a simple LPF in the Laplace domain: $$F(s) = \frac{\omega_0}{s+\omega_0}$$
Then it can be seen that it is stable for any value of $\omega_0\geq0$
Next, lets take a standard FIR system in Z-Domain:
$$H(z) = \frac{b_0 + b_1z^{-1} + ... + b_Nz^{-N}}{1} = \frac{b_0z^{N} + b_1z^{N-1} + ... + b_N}{z^{N}}$$
Since the parameters do not affect the denominator, a FIR filter is always stable!

Is stable to parameter tunings the same as "being continuos function"

If your filter is linear, any continuous input will yield a continuous output, regardless of the parameters used. Any discontinuous input might yield an discontinuous output, but it remains stable. If your filter is non-linear, the methods to assess stability are much harder (Lyapunov equations is one I know), especially to link them directly to parameter tuning (you could trigger a bifurcation point where the filter is stable for some values, but unstable for others).
I hope this clarifies a bit and I am not currently explaining the obvious.


How do I keep a value between known bounds in an open system?

In your plot of nb vs. image it appears that there is some high-frequency "noise" (no surprise there), and that the windowed mean is surprisingly stable.  For the purposes of this answer I am going to assume that the high-frequency noise is fairly representative of the general case, but that there will also be some slow-to-medium trends as well.  It seems like there would have to be if the images are changing.  The frame rate should generally be faster than the changes, though, so the trends should be relatively slow.

Normally I would suggest a low-frequency filter to get rid of the high-frequency "noise", but in this case I don't think that would be wise.  For any given image the next image is very likely to be highly correlated, so I think it makes sense to adjust somewhat if an image has a large change in nb.  I suggest two different solutions.

1) Use a simple moving average of the past n nb values.  Performance is obviously important in this application, and a moving average is computationally cheap.  A moving average gives some stability, which is important given the high-frequency noise, but also is responsive to changes in the data.  You can trade responsiveness for stability by decreasing or increasing n.

Once you have the nb moving average you use that to produce the "error" which is your target nb minus the moving average nb.  Thus, if the average nb is too low, the error is positive, if the average nb is too high, the error is negative.  You use the error to adjust the threshold.  Typically the new threshold would be

$threshold = threshold - \alpha * error$

$\alpha$ needs to be determined either by analysis of your problem or by empirically testing.  You substract the error result from the threshold because a positive error means you need nb to go up, which means that you need the threshold to go down.

2) Use a PID controller.  The PID controller would produce a threshold value for every image, and use the target nb vs. actual nb error (similar to the moving average error except it's an instantaneous error, not an average) to update its threshold for the next image.

The "integral" portion of the controller will give you stability, while the proportional/derivative part will give you responsiveness.  You would need to make sure that the loop itself is stable, and test the constants on different videos, but this would be a very powerful way of adjusting the threshold.  It also offers the possibility of making your app smarter by changing the loop constants dynamically, if needed.

Solution 1 is easier to implement and is a little cheaper computationally, solution 2 is more powerful but also more complicated.  I would see if the moving average is good enough.  If it is, great, if not, implement the PID controller.

EDIT: Tried to make the two solutions a little clearer.


How will the signal $\sum_{n=-\infty}^\infty \Delta (t-n)u_{-1}(t-n)$ look like?

Initially, don't worry about the sum. Just try to figure out how the signal $\Delta(t)u_{-1}(t)$ looks like. As soon as you know that, you just need to shift it by integer values of $t$, and add them all up. If the shifted versions do not overlap, you simply get copies of the basic signal shifted to integer values of $t$. From the solution in the solution manual (as given in a comment), the basic signal is the right half of a triangle, and the period equals $1$. The symbol $u_{-1}(t)$ refers to the unit step function, which is more frequently referred to as $u(t)$ (see also this question and its answer(s)).


difference equation block diagram delay or advance confusion

It depends from which side of the delay block you look at. Below is an example one sample delay from left to right in the arrow direction :

$$ x[n] \longrightarrow \boxed{ z^{-1} }  \longrightarrow y[n] $$

Then either $y[n] = x[n-1]$ or $x[n] = y[n+1]$ are both correct and equivalent expressions.

The first one is written when you trace the delay block in the arrow direction. And you get a delay effect from input x to output y. This is the normal operation direction.

The latter expression is obtained when you trace the delay block in the opposite direction of the arrow. Which gives you an advance effect from signal y to signal x.


You are talking about symmetry in two contexts :-

Gaussian
Window

I will talk about these two contexts one by one

Gaussian

A symmetric 2-D gaussian (I assume 1-D as in 1-D the question of symmetry does not arise) means that it does not have any directional preference (dictated by equal variances $\sigma_x$ and $\sigma_y$) . Now the job of Gabor Transform is to give us a joint time-frequency distribution (i.e telling us about frequency and phase content of a local image patch)
Now consider the following two cases :-

Symmetric Guassian ($\sigma_x$ = $\sigma_y$)
In this case the gaussian does not have any directional preference.

In the above image, suppose you are centered on the central pixel and consider a patch around it which is indicated in the above image by a black rectangle in the center. Following points are important for this image :

Doing a Gabor Transform over that patch would yield a time-frequency representation, which involves all the pixels in the patch accordingly weighted by the symmetric gaussian.
All the directions at a certain distance from the origin get the same weigthage.

However, if a gaussian oriented at an angle of $135^{\circ}$, then :-

the white principle diagonal pixels would get almost all the weightage while the other pixels will get little weightage.
The time-frequency representation will be mainly a time-frequency representation around the principle diagonal.
If the values of $\sigma_x$ and $\sigma_y$ are sufficiently small, the gabor transform will be very very close to the time-frequency representation for the straight line which marks the principle diagonal.

Thus, for the Gaussian function,

Symmetric gaussian is non-directional.
Non-Symmetric gaussian is directionally oriented.
Non-Symmetric gaussian will be useful in cases where the time-frequency analysis of a directionally oriented feature in a patch has to be carried out.

Window

The Window size basically represents the time-resolution (or space-resolution in the case of 2-D signals like images). A narrow window means that the focus on some local patch of the image is high and hence a more faithful time resolution of the Time-frequency analysis will be possible.

However, this is limited by the Heisenberg's Uncertainty principle. Although, Gabor Transform attains the theoretical lower bound for the uncertainty principle, it stays within the bounds of the principle. Thus, changing the window size is akin to changing your time resolution.

A Non-symmetrical window (like $3\times7$) is generally not chosen because generally we do not have a reason or a preference for changing our space-resolution in one direction (x or y) relative to other. In case there is a concern for the same, it can be done. However, I have never till now seen such a situation.


Matlab: Radix 2 FFT with 8 Channels and High Number of Inputs

Your problem is probably with twiddle factor adjustment. Below is a working code for 8 channel implementation of FFT (where the last 3 stages are explicitly computed, but not in a butterfly structure). Its derived based on the recursive structure of FFT. Where the splitting algorithm of N-point DFT into 2 N/2 point DFTs is applied for three consequtive stages.

Mathematical definition of the solution is as follows: Assume that $N$-point DFT $X[k]$ is split into two $N/2$ point DFTS according to the even and odd partitioning as usual:
$$ X[k] = X_e[k] + W_N^k X_o[k] ~~~,~~~\text{ for } k=0,1,...,N-1 $$

where the twiddle factor is $W_N = e^{-j\frac{2\pi}{N}}$ and note that the partial DFTs $X_e[k]$ and $X_o[k]$ are of length $N/2$and periodic. Then we can further decompose those $X_e[k]$ and $X_o[k]$ recursively to reach

$$\begin{align}
X_e[k] &= X_{ee}[k] + W_{N/2}^k X_{eo}[k] ~~~,~~~\text{ for } k=0,1,...,N/2-1 \\
X_o[k] &= X_{oe}[k] + W_{N/2}^k X_{oo}[k] ~~~,~~~\text{ for } k=0,1,...,N/2-1 \\
\end{align} $$

And further split those $N/4$ point DFTs into $N/8$ points DFTS which will be the channel signal DFTS.

$$\begin{align}
X_{ee}[k] &= X_{eee}[k] + W_{N/4}^k X_{eeo}[k] ~~~,~~~\text{ for } k=0,1,...,N/4-1 \\
X_{eo}[k] &= X_{eoe}[k] + W_{N/4}^k X_{eoo}[k] ~~~,~~~\text{ for } k=0,1,...,N/4-1 \\
X_{oe}[k] &= X_{oee}[k] + W_{N/4}^k X_{oeo}[k] ~~~,~~~\text{ for } k=0,1,...,N/4-1 \\
X_{oo}[k] &= X_{ooe}[k] + W_{N/4}^k X_{ooo}[k] ~~~,~~~\text{ for } k=0,1,...,N/4-1 \\
\end{align} $$

Below is an implementation of the above splitting algorithm. Note that your definition of the signals X0_0,X_0_1...,X0_7 is actually correct according to this splitting mode. But the butterfly computations seems not right.

clc;clear all; clear all

% S0 - Generate N point input data sequence x[n] of a sine wave.
% --------------------------------------------------------------
N           = 128;          % Length of signal
M           = 8;            % Number of channels
x           = sin(2*pi*0.123*[0:N-1]);

% S1 - obtain the even-odd partition signals :
% --------------------------------------------
xe = x(1:2:N);              % EVEN part of N-point x[n]
xo = x(2:2:N);              % ODD  part of N-poibt x[n]

xee  = xe(1:2:N/2);         % EVEN part of N/2-point xe[n]
xeo  = xe(2:2:N/2);         % ODD  part of N/2-point xe[n]
xoe  = xo(1:2:N/2);         % EVEN part of N/2-point xo[n]
xoo  = xo(2:2:N/2);         % ODD  part of N/2-point xo[n]

xeee = xee(1:2:N/4);         % EVEN part of N/4-point xee[n]
xeeo = xee(2:2:N/4);         % ODD  part of N/4-point xee[n]
xeoe = xeo(1:2:N/4);         % EVEN part of N/4-point xeo[n]
xeoo = xeo(2:2:N/4);         % ODD  part of N/4-point xeo[n]
xoee = xoe(1:2:N/4);         % EVEN part of N/4-point xoe[n]
xoeo = xoe(2:2:N/4);         % ODD  part of N/4-point xoe[n]
xooe = xoo(1:2:N/4);         % EVEN part of N/4-point xoo[n]
xooo = xoo(2:2:N/4);         % ODD  part of N/4-point xoo[n]

%S1 - Compute 8 of those N/8 point FFTs for the eee-ooo signals :
% ---------------------------------------------------------------
Xeee = fft(xeee,N/8);
Xeeo = fft(xeeo,N/8);
Xeoe = fft(xeoe,N/8);
Xeoo = fft(xeoo,N/8);
Xoee = fft(xoee,N/8);
Xoeo = fft(xoeo,N/8);
Xooe = fft(xooe,N/8);
Xooo = fft(xooo,N/8);

% S1 - obtain N/4 point DFT's from 8 , N/8 point FFTs of eee-ooo signals :
% ------------------------------------------------------------------------
W4 = exp(-1j*2*pi*(4/N));      % twiddle factor for the 3rd stage (Wn/4)
Xee = [Xeee, Xeee]  + (W4.^[0:N/4-1]).*[Xeeo, Xeeo];
Xeo = [Xeoe, Xeoe]  + (W4.^[0:N/4-1]).*[Xeoo, Xeoo];
Xoe = [Xoee, Xoee]  + (W4.^[0:N/4-1]).*[Xoeo, Xoeo];
Xoo = [Xooe, Xooe]  + (W4.^[0:N/4-1]).*[Xooo, Xooo];

% S2 - obtain N/2 point DFT's from N/4 point FFTs of ee-oo signals :
% ------------------------------------------------------------------
W2 = exp(-1j*2*pi*(2/N));      % twiddle factor for the 2rd stage (Wn/2)
Xe = [Xee, Xee] + (W2.^[0:N/2-1]).*[Xeo, Xeo];
Xo = [Xoe, Xoe] + (W2.^[0:N/2-1]).*[Xoo, Xoo];

% S3 - obtain N point final DFT from two N/2 point DFTs of e-o signals :
% ----------------------------------------------------------------------
Wn = exp(-1j*2*pi/N);
X = [Xe, Xe] + (Wn.^[0:N-1]).*[Xo, Xo];

% SX - Display results
% --------------------
figure,stem(abs(X))
figure,stem(abs(fft(x,N)),'g');


What is the proper way to implement a real time spectrogram?

Spectrograms are normally computed using the Short-Time Fourier Transform (STFT), and can be done in real-time. It consists of splitting the signal into overlapped frames, using a window function (typically scaled cosine like Hann) and then computing the FFT for the frame.

The overlap and windowing is key to avoid 'spectral leakage', and sounds like what is missing in your case.


Solved it!
function [t, pulse] = PAM (array)
  syms t
  u(t) = heaviside(t);
  pulse = array(1)*u(t);

  for i = [2:length(array)]
    if array(i) > array(i - 1)
      pulse += abs(array(i-1)-array(i)).*u(t - (i-1));
    else
      pulse -= (array(i-1)-array(i)).*u(t-(i-1));
    endif
  endfor
endfunction


STFT on time varying signal with good time and frequency resolution

I think that you are splitting your signal $x[n]$ into $N$ STFT intervals which you are overlapping by 50%, then on each interval finding the Power Spectral Density (PSD) with Welch's Method.

If this is correct I think this method is valid,

but I think you might be confused by the 'overlapping' because there are two at play here. When you use Welch's method you are splitting the STFT interval into smaller intervals, from which the periodograms is calculated, and then these are averaged to give you your PSD for that STFT interval.

So there are some variables which can be changed here

The length of your STFT intervals, $N$
The overlap of your STFT intervals
The overlap of your welch's periodograms

The smaller the length of the STFT the less smooth your 'time resolution' will be because you are averaging over a smaller length of time which will allow you to pick out sudden changes in frequencies, but less data is obviously going to have some kind of negative effect (maybe lower frequency resolution, but im not sure)

I think how much you overlap the STFT is quite related to the 'time resolution' here because you can think of it like, for the first interval you've windowed you are finding the PSD over $[x[0] - x[N]]$ but then you could quite easily take the next interval to be $[x[1] - x[N+1]]$ these two clearly overlap depending on $N$, but its better to just think of you moving your STFT window along $x[n]$ and calculating the PSD at each one, so in this sense your 'time resolution' can just be the same as your sampling rate, but its smoothed by the overlapping.

The overlapping of the periodograms smooths the frequency resolution but also reduces noise

I'm not sure what Jains method is but if you wanted to look at higher frequency resolution methods I would look at instantaneous frequencies using the Hilbert-Huang Transform, see this paper for some comparisons


What is Super Resolution in the Context of Windowing?

The term super-resolution is used very loosely nowadays. But I think the following problem was the original idea, or at least, the one that made the term famous.
Suppose you have a scene and several observations of that scene, i.e. several frames of a video with slight camera motion between them. Super-resolution algorithms combine those several observation into one highly resolved image. This is possible because each observation might contain unique information not present in the others.
If you consider a certain frame as the key frame, this process can be interpreted as a MAP estimation since a priori (or extra) information is present in others frames.
Additionally, we can use natural image statistics as extra information to regularize the solution. Total variation regularization is a broadly used technique.
I suggest the paper Sung Cheol Park - Super Resolution Image Reconstruction: A Technical Overview.