calculate SNR of two signals with power spectral density

Since you only say your channel is nonlinear, but not exactly which math it follows, the only answer possible to give here is

1.-5.: By filling in the signals and noise in your device equation, and transforming the math until you get an expression for what you're looking for.

"nonlinear" doesn't really say anything about what's going on (aside from the transfer function not being linear; especially, questions like determination, time-variance, boundedness, memory... are left unanswered).

If you think about it, it can be shown that with the lack specification given, the statement above is the strictest statement possible.

Even incorporating the fact that you claim the power spectral density exists doesn't do us much good – it only tells us the output of the device is weak-sense stationary, and that might as well mean the thing is just generating random tones.


How can I find the transform function, amplitude and frequency response of a digital filter in python?

Assuming each element is the filter coefficients for the R, G, and B channels (I do not know if this is indeed the case here), you can get the frequency response of each channel using scipy.signal.freqz:

scipy.signal.freqz(c, 1, v)

Where c is the filter coefficients and v is the signal data for that channel.


How to account for multiple signals in the input to a Kalman filter

The kinematic filter (which is mathematically a polynomial Kalman filter) has the ability to track any signal as long as the sampling rate is high enough, because when it's so, the complex motion of the target becomes linearized and fits into the quadratic kinematic model being employed by the Kalman filter dynamic model but the disadvantage is the processing cost due to high sampling rate.

Also the polynomial (kinematic) Kalman filter will not distinguish between two sinusoids but instead track their sum ignoring the noise only, unlike a specifically designed sinusoidal filter which tracks the chosen sinusoid and ignores all the rest as disturbance. However the model needs to be exact and known in advance as its fundamental disadvantage. The figure below shows a simulation example of a sinusoidal Kalman filter with two sinusoidal signals added together plus uncorrelated zero mean white noise as the measurement signal. It can be seen that the sinusoidal Kalman filter tracks the chosen sine wave (which is used in deriving the system dynamics matrix) and ignores all the rest as disturbance including the added white noise and the other sinewave which was at an harmonic multiple of the tracked sinusoidal frequency.

The extended kalman filter is a compromise, it has some additional computational burden, specifically the per sample update of the linearized derivatives, but a compromise can be found by lowering the sampling rate compared to a kinematic filter's requirement.


A first note: Most modern text-to-speech systems, like the one from AT&T you have linked to, use concatenative speech synthesis. This technique uses a large database of recordings of one person's voice uttering a long collection of sentences - selected so that the largest number of phoneme combinations are present. Synthesizing a sentence can be done just by stringing together segments from this corpus - the challenging bit is making the stringing together seamless and expressive.

There are two big hurdles if you want to use this technique to make president Obama say embarrassing words:

You need to have access to a large collection of sentences of the target voice, preferably recorded with uniform recording conditions and good quality. AT&T has a budget to record dozens of hours of the same speaker in the same studio, but if you want to fake someone's voice from just 5 mins of recording it will be difficult.
There is a considerable amount of manual alignment and preprocessing before the raw material recorded is in the right "format" to be exploited by a concatenative speech synthesis system.

Your intuition that this is a possible solution is valid - provided you have the budget to tackle these two problems.

Fortunately, there are other techniques which can work with less supervision and less data. The field of speech synthesis interested in "faking" or "mimicking" one voice from a recording is known as voice conversion. You have a recording A1 of target speaker A saying sentence 1, and a recording B2 of source speaker B saying sentence 2, you aim at producing a recording A2 of speaker A saying sentence 2, possibly with access to a recording B1 of speaker B reproducing with his/her voice the same utterance as the target speaker.

The outline of a voice conversion system is the following:

Audio features are extracted from recording A1, and they are clustered into acoustic classes. At this stage, it is a bit like having bags will all "a" of speaker A, all "o" of speaker A, etc. Note that this is a much simpler and rough operation than true speech recognition - we are not interested in recognizing correctly formed words - and we don't even know which bag contains "o" and which bag contains "a" - we just know that we have multiple instances of the same sound in each bag.
The same process is applied on B2.
The acoustic classes from A1 and B2 are aligned. To continue with the bags analogy, this is equivalent to pairing the bags from step 1 and 2, so that all the sounds we have in this bag from speaker A should correspond to the sounds we have in that bag from speaker B. This matching is much easier to do if B1 is used at step 2.
A mapping function is estimated for each pair of bags. Since we know that this bag contains sounds from speaker A, and that bag the same sounds but said by speaker B - we can find an operation (for example a matrix multiplication on feature vectors) that make them correspond. In other words, we now know how to make speaker 2's "o" sound like speaker 1's "o".
At this stage we have all cards in hand to perform the voice conversion. From each slice of the recording of B2, we use the result of step 2. to figure out which acoustic class it corresponds to. We then use the mapping function estimated at step 4 to transform the slice.

I insist on the fact that this operates at a much lower level than performing speech recognition on B2, and then doing TTS using A1's voice as a corpus.

Various statistical techniques are used for steps 1 and 2 - GMM or VQ being the most common ones. Various alignment algorithms are used for part 2 - this is the trickiest part, and it is obviously easier to align A1 vs B1, than A1 vs B2. In the simpler case, methods like Dynamic Time Warping can be used to make the alignment. As for step 4, the most common transform are linear transforms (matrix multiplication) on feature vectors. More complex transforms make for more realistic imitations but the regression problem to find the optimal mapping is more complex to solve. Finally, as for step 5, the quality of resynthesis is limited by the features used. LPC are generally easier to deal with a simple transformation method (take signal frame -> estimate residual and LPC spectrum -> if necessary pitch-shift residual -> apply modified LPC spectrum to modified residual). Using a representation of speech that can be inverted back to the time domain, and which provide good separation between prosody and phonemes is the key here! Finally, provided you have access to aligned recordings of speaker A and B saying the same sentence, there are statistical models which simultaneously tackle steps 1, 2, 3 and 4 in one single model estimation procedure.

I might come back with a bibliography later, but a very good place to start to get a feel for the problem and the overall framework used to solve it is Stylianou, Moulines and Cappé's "A system for voice conversion based on probabilistic classification and a harmonic plus noise model".

There is to my knowledge no widely piece of software performing voice conversion - only software modifying properties of the source voice - like pitch and vocal tract length parameters (For example IRCAM TRAX transformer) - with which you have to mess in the hope of making your recording sound closer to the target voice.


Extract breaths from respiratory signal, i.e. peak detection

I have a slick trick which may do what you want.

Do an exponential smoothing in the forward direction, call it F.

Do an exponential smoothing in the backard direction, call it B.

Take the difference of these two D = (B - F) / 2.  (Rescaling by half isn't necessary)

This will do two things for you:

1) Smooth out the noise

2) Convert peaks (and troughs) to zero crossings, no matter what their height.

You can find details in my blog article "Exponential Smoothing with a Wrinkle".  There is a "ramp up" and "ramp down" zone at the beginning and the end that you should not use.  They are fairly short and their length depends on your smoothing factor.

Hope this helps,

Ced


How to derive the transmitted signal from a constellation diagram

This is a basic question about how passband pulse-amplitude-modulation (PAM) works. It has nothing to do specifically with QAM, but it applies to any type of passband PAM (such as PSK, or any other choice of a constellation).

The question is how to obtain a transmit signal from a given point in a constellation diagram. Let's say that at time $t=kT$ we want to transmit a signal corresponding to the point $a_k$ in the constellation diagram, where $a_k$ is a complex number:

$$a_k=\text{Re}(a_k)+j\text{Im}(a_k)=|a_k|e^{j\phi_k}$$

The complex number is multiplied by a (real-valued) transmit pulse $g(t)$ and by a complex carrier, and the transmitted signal is the real part of that complex-valued signal:

$$\begin{align}s_k(t)&=\text{Re}\left(a_kg(t-kT)e^{j\omega_0t}\right)\\&=\text{Re}(a_k)g(t-kT)\cos(\omega_0t)-\text{Im}(a_k)g(t-kT)\sin(\omega_0t)\\&=|a_k|g(t-kT)\cos(\omega_0t+\phi_k)\end{align}\tag{1}$$

In practice you transmit a sequence of symbols, and the corresponding transmitted signal is given by

$$\begin{align}s(t)&=\sum_ks_k(t)\\&=\cos(\omega_0t)\sum_k\text{Re}(a_k)g(t-kT)-\sin(\omega_0t)\sum_k\text{Im}(a_k)g(t-kT)\\&=\sum_k|a_k|g(t-kT)\cos(\omega_ot+\phi_k)\end{align}\tag{2}$$

As you can see from $(1)$ and $(2)$, the transmitted signal can be represented either as two amplitude-modulated orthogonal carriers (sine and cosine), or as an amplitude and phase modulated carrier. In the first case the real and imaginary parts of $a_k$ determine the amplitudes of the two carriers, and in the second case the magnitude $|a_k|$ determines the amplitude of the carrier, and $\phi_k=\arg\{a_k\}$ is the phase of the carrier.


Signal to Noise ratio (SNR) is a ratio of powers and hence it is always greater than or equal to zero, it cannot be negative.

On the other hand, very commonly SNR is expressed in decibel (dB) notation

$$ \text{SNR}_{\text{dB}} = 10 \log_{10}\left( \frac{\sigma_x^2}{\sigma_n^2 } \right) $$

where $\sigma_x^2$ and $\sigma_n^2$ are the signal and noise powers respectively. In which case, a $0$ dB SNR means that the signal power is equal to the noise power. And when the signal power is less than the noise power, then one gets a negative SNR in dB...


Balance brightness and contrast between two images

The images you are working with are probably with gamma compensation, you should remove it first by applying inverse gamma.
More accurately, each camera has its own tone curve mapping. You can build a lookup table that converts one from another.
After you removed the non-linearity source, (whether you used method 1 or method 2), as a first order approximation, you can map colors from different cameras by using a 3x3 matrix applied on RGB.

Now you can compare your images.


How can I obtain the transfer function from these two transfer functions in this case?

The numerator and denominator of the total system are respectively the convolution of the numerators and the denominators of the filters to be cascaded.
$$
b(n)=b_1(n)*b_2(n)
$$
$$
a(n)=a_1(n)*a_2(n)
$$
You can validate it by the following code in MATLAB
b = conv(b1, b2);
a = conv(a1, a2);

sys1 = dfilt.df2t(b1, a1);
sys2 = dfilt.df2t(b2, a2);
sys = cascade(sys1, sys2);
[b, a] = tf(sys); % convert system object to coefficients


The cross shape that i see in this example seems to remain more or less orthogonal. Maybe using a houghlines transformation with a piecewise flexibility might provide more true positives, but this also increases the false positives. In this case a corner detector might also be a goodidea, thus providing a set of interest points to search from and eliminate. The problem again there would be the large number of false positives. Try this morphological corner detector (Robert Laganiere (1998): "Morphological Corner Detection") - its simple and structural. The problem in such cases is classifying your background on which the cross falls to help improve with the false positives!


I don't know if I understand exactly what you mean, because your question is very unclear (you may need to edit it), but I will try to answer what you may need in your task.

For DFT you will need this equation:

$ X[k] = \sum\limits_{n=0}^{N-1} x[n] \, e^{-j 2 \pi nk/N} $

where $k = 0, 1, . . . , N − 1$

But I guess you will be asking for IDFT more, so first you need to understand how IDFT works. I will try to explain this in my picture:

So the formula will look like this:

$ x[n] = \frac{1}{N} \sum\limits_{k=0}^{N-1} X[k] \, e^{j 2 \pi nk/N} $

where $n = 0, 1, . . . , N − 1$

Note that here, x[n] will be our complex time series, and X[k] is complex frequencies.

Now, we can get your sum, thanks to a simple code:

x = [6.29 8.11 -7.46 8.26 2.64 -8.04 -4.43 0.93 -9.15 9.29];
XF = fft(x);
k =  0:9;
YF = exp(-i*3*pi*k/5).*XF;
output = [sum(YF)];
disp(output)

And, as we can see, the sum equals 9.3.


will always periodogram produce correct frequencies?

A periodogram created from a finite number of samples using an DFT is not really a "correct" frequency estimator, as much as it is a (leaky) frequency energy binning sieve, where the width of each bin is set by the DFT length and the window function applied.

The width and spacing of each bin may or may not be narrow enough for you to consider an error of up to half that size to still result in a useful frequency estimate.  Various interpolation methods that post-process DFT results may allow more precise estimates of spectral peaks than the bin spacing.


Suppose that the sampling is done at a rate of $1000$ Hz, one sample every millisecond. Suppose also that the signal being sampled is at $3200$ Hz,
the first sample is at the peak of the sinusoid. The next sample,
will be taken one millisecond later during which time the sinusoid will
have gone through $3.2$ periods, and so the next sample will have the
same value as if the sinusoid had gone through $0.2$ periods, not $3.2$
periods. The one after that will be $0.4$ periods away from the peak, and
so on.  This is exactly the same set of samples that we would gave gotten
if we had been sampling a $200$ Hz sinusoid.  In one millisecond it would have
progressed through $0.2$ of its period of $5$ milliseconds and so on.
In other word, just by looking at the samples alone we cannot tell whether
the samples came from a $3200$ HZ signal or from a $200$ Hz signal.

If the signal being sampled was at $2800$ Hz, then we would get samples
corresponding to $0$, $-0.2$ of the period, $-0.4$ of the period and
so on.  But because sinusoids look the same in either direction in time,
these samples also look like they are the result of sampling a $200$
Hz signal.  This is the reason why the formula you are given, viz.

  Aliased frequency is the absolute difference between the actual signal frequency
  and the nearest integer multiple of the sampling frequency.

works to give you the right answer.


Translation property of 2-D discrete Fourier transform

Very informally, they use the property $e^{j\pi} = -1$ (which is correct), the property $e^{ab}={(e^{a})}^b$ (valid for certain values of $b$), and derive the final expression for $b=x+y$ when $x+y$ is an integer.

The pattern $(-1)^{x+y}$ looks like a checkerboard.


Time-domain LPF not showing expected behavior

Your math is correct but your code doesn't match it:

y = (wc*Ts*(x-xprev)+(2-wc*Ts)*yprev)/(2+wc*Ts)

It should be $x+xprev$ not $x-xprev$. The minus sign makes it a highpass and hence the attenuation.


Time-variance and causality of $y(t)=\sum_{k=-\infty}^{+\infty}x(t)\delta(t-kT)$

I agree with you on the system being a causal system since at a point in time $t$, the output is only made out of the input signal at the same point in time, weighted with an impulse (or not).

You can prove its time-variance by showing that a time-shifted input signal yields another output signal than the original input signal, for example:

$x(t) = \delta(t-\frac{T}{2})$ evaluates to $y(t)=0$ while
$x(t) = \delta(t)$ does not.


How many bits of precision are often used in fixed point DSP?

Almost all such programming is done in C, which means that 16 bits (shorts) and 32 bits (ints) make sense.  I doubt it will surprise you that most chips support those bit widths, with 16 bits being a particular favorite because it is wide enough to be useful and narrow enough to be power efficient, small, and fast.

I believe that there are a few chips that operate on 12 bits, but that is much rarer than the 16 bit chips.


Invariance of curvature in Curvature Scale Space

OK, so let's look at the curvature of some standard objects.

First, a circle:
$$
x(s) = \cos(s)\ ;\ \dot{x}(s) = -\sin(s)\ ;\ \ddot{x}(s) = -\cos(s)\\
y(s) = \sin(s)\ ;\ \dot{y}(s) = \cos(s)\ ;\ \ddot{y}(s) = -\sin(s)\\
$$
so the curvature is:
$$
\kappa(s) = \frac{\sin(s)\sin(s) + \cos(s)\cos(s)}{(\sin^2(s) + \cos^2(s))^{3/2}} = 1
$$

Rotating the circle gives:
$$
x'(s) = \cos(\theta) x(s) - \sin(\theta) y(s)\\
y'(s) = \sin(\theta) x(s) + \cos(\theta) y(s)\\
$$
so that
$$
\dot{x}'(s) = -\cos(\theta)\sin(s) -\sin(\theta)\cos(s) = -\sin(\theta + s)\\ \ddot{x}'(s) = -\cos(\theta)\cos(s) +\sin(\theta)\sin(s) = -\cos(\theta + s)\\
\dot{y}'(s) = \sin(\theta)\sin(s) +\cos(\theta)\cos(s) = \cos(\theta - s)\\
\ddot{y}'(s) = \sin(\theta)\cos(s) - \cos(\theta)\sin(s) = \sin(\theta - s)\\
$$
Even thought it looks like I got a sign wrong somewhere, this will also evaluate to 1.

So, at least for a circle, the formula is rotationally invariant.

The trouble comes when you try to implement this.

function [kappa,xdot,ydot,xdotdot,ydotdot] = curvature(x,y,delta_s)
    xdot = diff([x x([1 2])])/delta_s;
    xdotdot = diff([xdot])/delta_s;
    ydot = diff([y y([1 2])])/delta_s;
    ydotdot = diff([ydot])/delta_s;

    xdot = xdot(1:length(xdot)-1);
    ydot = ydot(1:length(ydot)-1);

    kappa = (xdot.*ydotdot - xdotdot.*ydot)./sqrt(xdot.^2 + ydot.^2 + 10^9).^3;

endfunction

The only way to get the "derivatives" of your parameterized curves will be to take differences.  For the circle, this isn't too bad, but for more complex shapes, rotations might make these values a little different --- especially if the shapes are generated from two different sets of measurements.

For example the plot below shows:

The original curvature (note problem at end of range)
The curvature of a circle rotated
The curvature of a noisy circle
The curvature of a rotated noisy circle

// 25885
delta_s = .01;
s = -%pi:delta_s:(%pi-delta_s);
N = length(s);

x = cos(s);
y = sin(s);

xn = cos(s) + 0.001*rand(1,N,'norm');
yn = sin(s) + 0.001*rand(1,N,'norm');

[k,xd,xdd,yd,ydd] = curvature(x,y,delta_s);
[kn,xdn,xddn,ydn,yddn] = curvature(xn,yn,delta_s);

theta = %pi/4;
rot = [cos(theta), -sin(theta); sin(theta) cos(theta)];

all = rot*[x;y];
xdash = all(1,:);
ydash = all(2,:);

alln = rot*[xn;yn];
xdashn = alln(1,:);
ydashn = alln(2,:);

[k2,xd2,xdd2,yd2,ydd2] = curvature(xdash,ydash,delta_s);
[k2n,xd2n,xdd2n,yd2n,ydd2n] = curvature(xdashn,ydashn,delta_s);

clf
subplot(221)
plot(s,k)
title('Curvature of a noiseless circle')
subplot(222)
plot(s,k2)
title('Curvature of a ROTATED noiseless circle')
subplot(223)
plot(s,kn)
title('Curvature of a NOISY circle')
subplot(224)
plot(s,k2n)
title('Curvature of a ROTATED NOISY circle')


Why does this transfer function has a second zero


Computer Vision - Gaussian Pyramid, Scale Space and Pattern Matching

Template matching is related to finding the maximum correlation between two signals, where one of the signals (the image being explored) is usually much bigger than the reference signal (the object being searched for). Obviously, the best correlation score can only be obtained when the size of the object in the image is equal to its size in the template (intuitively, both shapes will superpose).

  However I am not able to relate in any meaningful way the number of pyramid levels to create (for model and for target)

The number of scales is related to two things:

the size of the object of interest in the explored image
the computation time reduction that you want to obtain.

To have the best localization accuracy, you need to have one template image that matches exactly the size of the object in the explored image. Thus, if you don't know a priori this size, but you have instead a range of sizes, your template pyramid should include each of these sizes.

Since the template matching part can be expensive (one correlation computation for each possible position of the template in the image), you want to reduce the number of pixels explored. Thus, you create a pyramid of the image to explore, start by looking for the object at the coarsest scale to get a rough position that gets refined at the next scale, etc. Again, since the size of the object needs to match between the template and the explored image, you need to rescale appropriately each of the images from the pyramid that you built for the template (previous paragraph).

  For example 1 level of the pyramid reduces the resolution of the image by 2 (is this the same as saying that the image has been scaled by 0.5 ) ?

It depends... Reducing the resolution by two along each dimension divides the number of pixels by 4. This is usually referred to as scaling by 0.5 but as an image resolution reduction by 4.

A common practice is to divide the number of pixels by 2, yielding a $\sqrt{2}$ factor along each dimension.

  Than what about a pyramid level of -1 (which should make the image bigger x 2.0) : with what gaussian filter should it be treated before upsampling ?

You need to use image interpolation techniques here. OpenCV uses (in optical flow methods) a gaussian upsampling by injecting $0$'s along the odd lines/columns then applying a gaussian blur.


A discrete sine wave is defined as
$$\sin (2\pi  f n)$$
Now your index is n=2*pi*(0:1/200:1). This means your sampling frequency is 200Hz.


The problem is that the following statement- "So now for the spectrum between 0 to 100 MHz the vector s⃗ l contains amplitudes for 20000 frequencies (D=20000)"- is faulty.

A signal that lasts $5\mu s$ is captured at 4 GS/s.  That is a total of 5e-6 * 4e9 = 20000 samples.  The Nyquist frequency is half the sample rate, or 2 GHz.  They take the Fourier transform of the samples ($F(w,l)$).  At this point, assuming the samples are real (as opposed to complex), there are 20000 frequency bins representing -2 GHz to +2 GHz.  The negative frequency information is redundant since it mirrors the positive frequencies.  I believe that when they said "We then remove from F(ω,l) its DC component and the redundant part of the spectrum" they were talking about getting rid of the 0 Hz and negative frequency information.

That would leave you with 9999 frequency bins representing 200 kHz to 2 GHz.  To make life easier, though, let's bring back the DC bin and say that we have 10000 frequency bins representing 0 Hz to 2 GHz.

Hopefully it is clear now why "So now for the spectrum between 0 to 100 MHz the vector s⃗ l contains amplitudes for 20000 frequencies (D=20000)" is faulty.  There are 10000 bins and they represent 0 to 2 GHz, not 0 to 100 MHz.  They now bandpass the signal down to 4 MHz bandwidth.  $(4e6 / 2e9) * 10000 = 20$ frequency bins.  I believe they get 40 dimensions from those 20 bins by counting both the amplitude and phase of the bins as features, $20 * 2 = 40$.


How is RMS and Peak-to-peak of white noise calculated

three things to think about:

"white" describes the power spectrum, but to know the peak values, you need to know the probability density function (p.d.f.) of the noise process.  often, with "white noise", we are assuming "normal" or "gaussian" p.d.f.
even for a finite power (which means a finite variance for the p.d.f.), a gaussian random variable has infinite peak value (just that it hits whatever arbitrary high value with decreasing probability).  so for gaussian, you need to define like a "1% likelihood peak value".
true white noise (which does not exist in reality) has infinite power and infinite variance.  this is because true white noise has infinite bandwidth.  so, in reality, we deal with bandlimited "white" noise, so the bandwidth might be a necessary parameter to know, before you can think about its r.m.s. (which is the square root of the variance, assuming no DC component) and peak.


How to write time waveform of transmitted sequence with certain chip rate?

While a reference to where you got this from would be helpful, I have a theory on where the term in question comes from. Note its form for increasing $n$:

$$
\exp\left({j\frac{\pi n}{2}}\right) = [1, j, -1, -j, 1, j, -1, -j, \ldots ], \ n = [0, 1, 2, 3, 4, 5, 6, 7, \ldots]
$$

This is a discrete-time complex sinusoid whose frequency is one-fourth of the sample rate. In digital communications signal processing, this is a very common intermediate frequency selection when using real-valued sampling. It is often chosen because the center of the desired passband lies in the center of the single independent sideband that you get when using real samples. Shifting the resulting signal to complex baseband is often a subsequent processing step.

In your case, the equation that you gave may express that you have a chip sequence that is modulated onto a complex-valued carrier at one-fourth of the sample rate. Since you don't indicate what the values of the signal are at time instants other than $nT_c$, it's hard to say whether this is actually the case.


Kullback-Leibler Distance of Spectral Data

The expression in your question seems to have been written for univariate Gaussians. For multivariate Gaussians,
$KL(A,B) = \frac{1}{2}\left[\log\frac{|\Sigma_B|}{|\Sigma_A|} - d + Tr(\Sigma_B^{-1}\Sigma_A) + (\mu_B - \mu_A)^T \Sigma_B^{-1}(\mu_B - \mu_A)\right]$
where $d$ is the dimensionality of the feature vectors.

(https://stats.stackexchange.com/questions/60680/kl-divergence-between-two-multivariate-gaussians)

So

\begin{eqnarray*}
KL2(A,B) &:=& KL(A,B)+KL(B,A)\\
&=& \frac{1}{2}\left[ - 2d + Tr(\Sigma_B^{-1}\Sigma_A) + Tr(\Sigma_A^{-1}\Sigma_B) + (\mu_B - \mu_A)^T (\Sigma_A^{-1}+\Sigma_B^{-1})(\mu_B - \mu_A)\right].
\end{eqnarray*}

I don't think you compute this for the spectrograms themselves but for feature vectors (some of the features may be derived from the spectrogram).


Crop all written letters from single image

Presumably this is an academic assignment of some sort and you are expected to use a specific tool, OpenCV I am guessing. You should review the documentation for that tool and try to find a feature which matches the requirements. I suspect you would need to use the Image Segmentation feature: https://docs.opencv.org/3.4.3/d3/db4/tutorial_py_watershed.html


morphological processing : Hit or Miss shape detection

They are using the same method. The second resource is just showing the windowing function. Look at slide 12, the windowing function is shown here.

Erosion implies you use a windowed kernel, and slide it around the entire image (as seen in this gif), so the first source doesn't mention this step. The 2nd source is explicitly stating you must pad the search space, and then slide the kernel around it.Hope that makes sense


Calculating power in frequency band using power spectral density

It's less about power spectrum and more about Fourier transform of real signals. Let's take an easy signal, $g(x) = \cos(10\pi x)$.

Then:
$$ F(\omega) = \frac{\delta(\omega -10) + \delta (\omega +10)}2.$$

When you ask "what is the energy in the band $ \omega = [9,11]$ you ask what is the energy of that cosine (of course this is an example).

If you look only at the positive side of the frequency you lose half the energy.


Sound synthesis effects for natural sounds

First, grab Sonic Visualizer, it is much better than Audacity for looking at sounds.

What you see here is probably the result of the sum of two simple and stationary sounds at fundamental frequencies close to each other. This causes beating of their fundamental, causing the amplitude modulation (tremolo) you observe.

Two important factors make a synthetic sound "feel" natural:

Whether it accurately reproduce the sound of the original instrument over time. Some criteria to look at: timbre (distribution of harmonics), timbre modulations over time, amplitude modulation over time (presence of tremolo, overall amplitude envelope), pitch modulations (presence of vibrato, does the instrument instantly reaches the target note or is there a short transition like on brass instruments...).
How it responds to control by the performer. You have to list which parameters, along with note data, will be available to your system (will it synthesize sounds from a bare MIDI file? From a MIDI keyboard sensitive to velocity and aftertouch? From a dedicated hardware device emulating an accordion, with accelerometers and pressure sensors?), and understand how to map these inputs to synthesis parameters. You also have to list which controls the performer will have over the system (for example, is there an option for the user to switch on/off the different registers?).

I will mention 3 synthesis approaches you can follow for this problem:

Sampling. Record, one by one, all the notes of an accordion, under different performance parameters (how hard the below is moved for example, all different combinations of registers). Play back from this exhaustive database of sounds, with loops to make them last as long as you want. Advantages: a note taken in isolation sounds indistinguishable from the real thing. Drawbacks: might sound "fake", and will not respond to performer input (not "expressive") unless you start doing some signal manipulation on the recorded sample - at which point they won't sound as natural as they were recorded. A lot of disk/memory used by the samples. Warning: it's quite an involved job to exhaustively sample an instrument and organize all the data, and there are already companies very good at doing this kind of stuff - it'll be hard to beat them.
Classic subtractive synthesis. You could start with two square or pulse waves, slightly detuned to get beating, a tad of high-pass filtering to give it a "nasal" quality, and an EQ to color this... Then stack several of those to get different registers. Check in the patch libraries of classic "old" synths (JX8p, D50, M1)... they all have decent accordion sounds which are totally sample-free. Advantages: can be made very expressive - since you have control over the whole production of the sound, you can easily map input parameters on sound parameters. You won't have to reinvent the wheel and can base yourself on existing sound synthesis code. Inconvenient: will sound "synthetic" (but "lively & synthetic" is sometimes better than "realistic & dead"!).
Physical modeling. You could get something more realistic than vanilla subtractive synthesis by using physical modeling - try cobbling together (for example in csound, Max/MSP or supercollider) a few reed models (the 3 mentioned programs have basic primitive for physical modeling - but I'm not sure they have free reeds), add EQ or a simple resonator. Ideally, you could get something both realistic sounding and controllable, though it's not the easiest path!


Since you stated "external to the loop BW" I assume you mean specifically an output signal, not the signal fed back to the phase detector but provided to an external sink. Then yes, in this case a frequency divider would likely have phase ambiguity regardless of the division number unless special provisions are made to synchronize state (which requires an external time reference). This is because the division would start on any arbitrary edge (in typical cases the rising edge, but the salient point is any edge). Without an external synchronization of state, two identical PLL's may start the divided signal on different edges resulting in complete phase ambiguity. To eliminate this, some sort of state synchronization or reset is required.

Consider the operation of the divide by 2 as shown in the figure below. This typically operates by having the output change state on every rising edge of the input. Depending on the starting state of the local circuit, either of the two lower plots in the figure could result. Unless the circuit has a specific reference to common time there would be no way for it to remain in phase with its prior phase reference (were it to stop and start again, or considering two PLL circuits etc). This is what is meant by phase ambiguity. From this it should be clear that this would apply to any number N.


Plotting a histogram of particle intensity vs pixel

The expected distribution of intensity versus x,y position is formed by two idependent and identically distributed gaussian random variables; one for each axis. The distribution shown is for one "slice" through this 2D distribution.  The top view of the distribution can be envisioned as a 3d Gaussian "hill", and the histogram shown is for one slice through this hill as viewed from the side, where the slice can be from any direction as long as it goes through the center (origin).


For the amplitude, the inequation $|a+b|≤|a|+|b|$ is satisfied.

However, since the power includes a square calculation on the amplitude, you can find all the three ($=$,$\leq$,and $\geq$) relations.

The difference is that amplitude is a vector, while power is a scalar.


Coherent detection in radar and radio signals

In the context of wireless systems, coherent detection has the definition of detection when the wireless channel is known at the receiver and thus matched filtering is possible and infact optimal at low SNR.

The signals from different antennas at the radar receiver are coherently combined using maximum ratio combining. The MRC is nothing but the matched filter.


How to find the inverse Z transform of this function in z domain?

I think you are using the partial fraction expansion in wrong way. You have to know for $z$ transform we expand it as a function of $z^{-1}$ and not $z$. for your special problem factor out  $z^2$ from denominator then divide the nominator by that.

Then you could factor your polynomial over a real field or complex field. This means after factorization over complex field you may obtain polynomials with complex coefficient but factorization over a real field you obtain polynomials with real coefficients. so for your case the denominator is irreducible over real field but we could factor it over a complex field, which lead to a pair of conjugate poles (roots). Now you could easily obtain the inverse using the inverse of a single pole (don't worry about the complex numbers, after summing response of 2 conjugate pole the complex parts will cancel out each other and only a real part remain).
Also I have to say, for your case you will reach 2 answer considering the region of convergence (RoC) in $z$ space.


Why is interpolation a time varying system

Considering a discrete time framework, the interpolation operation has two stages; first expand the signal by zero stuffing between the existing samples, and then lowpass filter the expanded signal to get the interpolated samples. The second operation is LTI, hence it's the first stage which is not LTI.

The expansion block is shown like $$x[n]\longrightarrow \boxed{ \uparrow L } \longrightarrow y[n]$$ and mathematically defined as

$$ y[n] = \begin{cases}{ x[\frac{n}{L}] = x[m] ~~~ , ~~~ \text{ for } n = m L~~ ,~~  m=0,\pm 1, \pm 2,... \\  ~~~~~ ~~~~ ~~~ 0  ~~~~~ ~~~~ ~ , ~~~ \text{ otherwise } }\end{cases} $$

Now apply the test for time-invariance  on $y[n]$ :

First shift the ipnut. Let $x_d[n] = x[n-d]$, then call $y[n]$ as $y_2[n]$ as:

$$ y_2[n] = \begin{cases}{ x_d[\frac{n}{L}] = x_d[m] ~~~ , ~~~ \text{ for } n = m L~~ ,~~  m=0,\pm 1, \pm 2,... \\  ~~~~~ ~~~~ ~~~ 0  ~~~ ~~~~~~~~~~~~ ~, ~~~ \text{ otherwise }}\end{cases} $$

which becomes:
$$ y_2[n] = \begin{cases}{ x[\frac{n}{L} - d]  = x[m-d] ~~~ , ~~~ \text{ for } n = m L~~ ,~~  m=0,\pm 1, \pm 2,... \\  ~~~~~ ~~~~ ~~~ 0   ~~~ ~~~~~~~~~~~~~~~~~ ~ , ~~~ \text{ otherwise } }\end{cases} $$

Second, similarly shift the output $y[n]$ by $d$ which is denoted as $y_d[n] = y[n-d]$:

$$ y_d[n] = y[n-d] = \begin{cases}{ x[\frac{n-d}{L}]=x[m] ~~~ , ~~~ \text{ for } n-d = m L~~ ,~~  m=0,\pm 1, \pm 2,... \\  ~~~~~ ~ ~~~~ ~~~~ 0   ~~~ ~~~~~~~~~~ ~, ~~~ \text{ otherwise }  }\end{cases} $$

which becomes:
$$ y_d[n] = \begin{cases}{ x[\frac{n-d}{L}] = x[m] ~~~ , ~~~ \text{ for } n = m L + d~~ ,~~  m=0,\pm 1, \pm 2,... \\  ~~~~~ ~~~~ ~~~ 0  ~~~ ~~~~~~~~~~ ~, ~~~ \text{ otherwise } }\end{cases} $$

As you can see, $y_2[n]$ and $y_d[n]$ are not the same signals. You can explicitly see this by evaluating them for specific values of $n$ and $d$, such as $d=1$ and $n = L$, then $y_d[L] = 0$ while $y_2[L] = x[0]$ and hence $y_d[L] \neq y_2[L]$. And therefore the conclusion that the expansion system (or as a consequence the interpolation system) is a Time-Varying and not LTI system.


If a signal is sparse in the time domain, it has infinite support in the frequency domain, thus is not band-limited in the mathematical sense. However it might be so close to band-limited that the out-of-band spectrum disappears under the quantization or other noise floor.


Analyzing a particular discrete-time LTI system for input signal $x[n]=(1/3)^n$ for *all* $n$

This problem tests the understanding of two basic properties of (discrete-time) LTI systems. First, for an input signal of the form $x[n]=z_0^n$, $-\infty<n<\infty$, the output is given by

$$y[n]=z_0^nH(z_0)\tag{1}$$

where $H(z)$ is the system's transfer function, and where it is assumed that $z_0$ is inside the region of convergence (ROC) of $H(z)$. [Note that $(1)$ is not true for $z_0$ outside the ROC.]

The second property is simply the $\mathcal{Z}$-domain input-output relationship

$$Y(z)=X(z)H(z)\tag{2}$$

Applying $(1)$ to your example you can conclude that

$$H\left(\frac{1}{3}\right)=\frac23\tag{3}$$

and from $(2)$ you can get an expression for $H(z)$ in terms of the constant $a$. This constant can be determined from Eq. $(3)$.


Remove a known wav file from recorded file

Assuming that there are no non-linear effects between the sound N source and the microphone (such as AGC), you might try to estimate the impulse response of the channel between the source for sound N and the microphone data.  The better estimate you can make of that impulse response, the more of sound N that you can remove.  It’s unlikely that the negative of the impulse response is just the data inverted. That's due to the fact that the amplifier, speaker, microphone, and room acoustics (etc.) probably have frequency dependent phase delays and/or multiple reverb/echoes.

Channel estimation in noise is a non-trivial problem.


First, it is instructive to understand the normalised frequency axis of both the fft and the freqz.

By default, in fft, the frequency axis is normalized by the sampling frequency $f_s$. On the other hand, in freqz, the frequency axis is normalised by $f_s/2$, it seems. Having said that, we note that the radian frequency of $\pi/8$ corresponds to the normalized frequency of $\frac{1}{16}=0.0625$ ($2 \pi \frac{f}{f_s}=\frac{\pi}{8}$) in fft and $\frac{1}{8}=0.125$. In  your case, the frequency response at $\pi/8$ seems to be about $-90$ dB.

fft does not provides the mirror image of the freqz. Freqz only shows the half portion of the fft spectrum.


How can I get the power of a specific frequency band after FFT?

To get the total power across bins, sum the power in each bin. Also you need to compensate for your window loss if you want an accurate result.

For a rectangular window, the power in each DFT bin is given as:

$$P = \frac{Y[k]Y[k]^*}{N^2} = \bigg(\frac{|Y[k]|}{N}\bigg)^2$$

and $$P_{dB} = 20Log_{10}\bigg(\frac{|Y[k]|}{N}\bigg)$$

The total power over a range of such bins would be:

$$P_{k} = \frac{\Sigma_k(|Y[k]|^2)}{N^2}$$

Windowing Loss

When we apply a window, the coherent gain and non-coherent gain are affected differently. If the N samples are coherent (such as a frequency completely contained in a bin), the output grows at rate N (direct summation). If the N samples are non-coherent (spread across bins), the output grows at the rate of $\sqrt{N}$ (rms summation).

If we are only concerned with the signal power in each bin, we would use the coherent gain which is given as

$$G_c = \frac{\Sigma w[n]}{N}$$

where $w[n]$ are the window coefficients.   (For example the Blackman(30) window has a coherent gain of 0.406 or -7.83 dB):

Likewise non-coherent samples grow at $\sqrt{N}$. The non-coherent gain refers to the change in DFT output level of the non-coherent samples (such as white noise) due to the window function.

The normalized non-coherent gain is the rms value of the window weights given as:

$$G_{nc} = \sqrt{\frac{\Sigma w[n]^2}{N}}$$

The SNR change or processing gain of the window would be the ratio of the two, in dB this would be:

$$PG = 20log_{10}\bigg( \frac{G_c}{G_{nc}} \bigg)$$

If you are only interested in the signal power level and not signal to noise ratio, so you would simply compensate for the power computed by using $G_c$.

It may confusing to some how the power for a single tone can be effected differently by the window than the power in a signal that is spread evenly across multiple tones as described above. An intuitive explanation of this is that windowing the time domain function results in each bin in the FFT having an equivalent noise bandwidth (ENBW) that exceeds one bin. A rectangular window in comparison has an ENBW of exactly one bin, so when we sum the power of all the bins in the frequency domain, Parseval's theorem proves correct in that the total power from that summation is equal to the total power in the time domain signal. However if we window in time with something other than a rectangular window, the effect of the wider ENBW means that each bin in frequency is now reporting the power for that bin PLUS adjacent bins. So in this case when we sum the reported power over all bins, the total summation will over-estimate the actual noise power. In a case of a single tone, or tones that are spaced wider than the ENBW, each bin will accurately report the power for each tone when we have scaled the result by the recripocal of $\sum W_n$ where $W$ is the windowing funciton. (for the rectangular window this scaling results to $1/N$).


Problem with Covariance matrix using diagonal loading involved in calculation of eigenvalues

You can write

$$ R=YY^H $$
where $Y$ is a matrix of size $N\times N_f$ and $N$ is the dimension of $y_k$. $Y$ contains all the measured $y_k$ as its columns.

Then, the rank of $R$ is upper bounded by $N_f$. In particular, if $N_f<N$, $R$ will always be a singular matrix. So, if you have too few measurements, you will likely run into the problem you described.

For your second part, how do the Eigenvalues of $R$ and $R+\sigma I$ relate:
Consider the eigenvalue decomposition of $R$ by
$$ R=PVP^{-1}$$
where $P$ contains the eigenvectors and V is the diagonal matrix with the eigenvalues (they can be zero). Then you have

$$ R+\sigma I=P V P^{-1} + \sigma I = P(V+\sigma I)P^{-1} $$

This means that the eigenvalues are all increased by $\sigma$.


You have a typo in your definition of $R_y(k)$, and an error in the time indices when developing the equation. The proper definition of autocorrelation for complex signals (in the case of wide-sense stationary processes) is
$$R_y(k)=\mathbb{E}\left[y(n)\overline{y(n+k)}\right],$$
and setting $k=0$, we obtain
$$\begin{align}
R_y(0)&=\mathbb{E}\left[y(n)\overline{y(n+0)}\right]=\mathbb{E}\left[\left|y(n)\right|^2\right]\\
&=\mathbb{E}\left[\left|\sum_{m=-\infty}^{\infty}h(m)x(n-m)\right|^2\right]\\
&=\mathbb{E}\left[\left(\sum_{m=-\infty}^{\infty}h(m)x(n-m)\right)\overline{\left(\sum_{i=-\infty}^{\infty}h(i)x(n-i)\right)}\right]\\
&=\mathbb{E}\left[\sum_{m=-\infty}^{\infty}\sum_{i=-\infty}^{\infty}h(m)\overline{h(i)}x(n-m)\overline{x(n-i)}\right]\\
&=\sum_{m=-\infty}^{\infty}\sum_{i=-\infty}^{\infty}h(m)\overline{h(i)}\,\mathbb{E}\left[x(n-m)\overline{x(n-i)}\right]\\
&=\sum_{m=-\infty}^{\infty}\sum_{i=-\infty}^{\infty}h(m)\overline{h(i)}\,R_x(m-i)\\
\end{align}$$
where I first use the fact that the product of sums is the sum of products (and I use a different summation index to better show the cross-products), and second, that the impulse response of the system is deterministic and therefore constant for the purposes of taking expectation.
Now, we can group those terms where $m=i$ and those where $m\neq i$, obtaining
$$\begin{align}
R_y(0)&=\sum_{m=-\infty}^{\infty}\sum_{i=-\infty}^{\infty}h(m)\overline{h(i)}\,R_x(m-i)\\
&=\sum_{m=-\infty}^{\infty}h(m)\overline{h(m)}\,R_x(0)+\sum_{m=-\infty}^{\infty}\sum_{i=-\infty\\ i\neq m}^{\infty}h(m)\overline{h(i)}\,R_x(m-i)\\
&=R_x(0)\sum_{m=-\infty}^{\infty}|h(m)|^2+\sum_{m=-\infty}^{\infty}\sum_{i=-\infty\\ i\neq m}^{\infty}h(m)\overline{h(i)}\,R_x(m-i)\\
\end{align}$$
Finally, we know that the input noise is white, and therefore $R_x(k)=0$ for $k\neq0$. We also know that $R_x(0)=\sigma_x^2$. Thus, the second term vanishes, and we obtain
$$\begin{align}
R_y(0)
&=R_x(0)\sum_{m=-\infty}^{\infty}|h(m)|^2=\sigma_x^2\sum_{m=-\infty}^{\infty}|h(m)|^2
\end{align}.$$


Convolution in frequency domain, what is the output bandwidth

Following a time-domain convolution analysis which relates the output signal's domain of support to the input signals' domain of supports, we could reach a similar conclusion:

Given two frequency domain signals with domain of supports (for positive frequencies) such that $X(\omega) \neq 0$ for $\omega_1 < \omega < \omega_2$ and $H(\omega) \neq 0$ for $\omega_3 < \omega < \omega_4$, then their convolution $Y(\omega) = X(\omega) \star H(\omega)$ will have the domain of support as $Y(\omega) \neq 0$ for $$ \omega_1 + \omega_3 < \omega < \omega_2 + \omega_4 $$

Therefore the bandwidth of $Y(\omega)$ is related to the bandwidths of $X(\omega)$ and $H(\omega)$ as:

$$ BW_Y =  (\omega_2 + \omega_4) - (\omega_1 + \omega_3) $$
$$ BW_Y =  (\omega_4 - \omega_3) + (\omega_2 - \omega_1) $$
$$ BW_Y =  BW_H + BW_X $$

This assumes that the domains were compact and finite. If they are either noncompact or non-finite the result should be modified accordingly.


Periodicity of Fourier series's coefficients

The reason why Fourier series coefficients of continuous functions are generally not periodic is because of the continuity of the function. The discrete Fourier series coefficients are periodic because the analyzed signals are discrete. Note the duality relationship of the Fourier transform

$$\text{periodicity}\Longleftrightarrow\text{discreteness} $$

Periodicity in one domain leads to discreteness in the other domain. E.g., a periodic function has only discrete frequency components (as shown by its Fourier series). If these discrete frequency components are periodic, then, consequently, the periodic function must be discrete, i.e., it must be non-zero only at discrete times.

In his answer, Fat32 gave an example of a periodic signal with periodic Fourier series coefficients. The most general form of such a signal is

$$x(t)=\sum_{m=-\infty}^{\infty}b_m\delta\left(t-m\frac{T}{N}\right)\tag{1}$$

where the coefficients $b_m$ are essentially the discrete Fourier series coefficients of the $N$-periodic Fourier coefficients $a_k$ of the $T$-periodic signal $x(t)$:

$$b_m=\frac{T}{N}\sum_{k=0}^{N-1}a_ke^{j2\pi km/N}\tag{2}$$
$$x(t)=\sum_{k=-\infty}^{\infty}a_ke^{j2\pi kt/T}\tag{3}$$

As mentioned above, the signal $x(t)$ is of course discrete in the sense that it is non-zero only at discrete time instances $t_m=mT/N$. Any periodic signal $x(t)$ with $N$-periodic Fourier coefficients $a_k$ must be of the form $(1)$.

A derivation of Equations $(1)$ and $(2)$ can be found in this answer.


The spectrum is that of the base pulse used in the modulation, so in this case a rectangular pulse. A single rectangular pulse in time, as given by the Fourier Transform, is a Sinc in Frequency with the first nulls spaced at $1/T$ away from the center of the Sinc, where $T$ refers to the duration of the pulse in time. (Just as that shown by the OP). If the pulse is repeated at a consistent rate, then then only spectral components would be weighted impulses spaced by the repetition rate, and the weights would be the same Sinc function shape shown. Thus a 50% duty cycle square wave would have its even harmonics in the nulls of the Sinc function shown, resulting in the well known odd-harmonic spectrum of impulses. If the data is random, such that a patter never repeats, then the spectrum will be continuous as shown here.


Determine type of filtering applied between two images

You can use Discrete Fourier transform. One of it's interesting qualities is that convolution in the time domain is translated into multiplication in the frequency domain.

Applying a filter, represented by the function $h(t)$, on a signal\image $x(t)$ is equivalent, in fact, to convolving the two in the time domain. You can learn more about convolution here. We will denote the output of the operation $y(t)$, and it will be represented by $y(t)=x(t)*h(t)$. Note that $*$ in this context is a convolution operation.

If we will apply a Discrete Fourier transform $\mathcal F$ on both sides we get:
$$\mathcal F[y(t)]=\mathcal F[x(t)*h(t)]=\mathcal F[x(t)]\cdot \mathcal F[h(t)]$$
Where the $\cdot$ denotes multiplication. This is again derived from the fact that convolution in the time domain is translated into multiplication in the frequency domain.

You can now use simple math to get:
$$\mathcal F[h(t)]=\frac{\mathcal F[y(t)]}{\mathcal F[x(t)]}$$
And hence:
$$h(t)=\mathcal F^{-1}\left[\frac{\mathcal F[y(t)]}{\mathcal F[x(t)]}\right]$$ where $\mathcal F^{-1}$ represent the inverse of DFT.

You can now use $h(t)$ for any new $x'(t)$ by applying:
$$y'(t)=x'(t)*h(t)\implies \mathcal F[y'(t)]=\mathcal F[x'(t)]\cdot \mathcal F[h(t)]$$

Sure DFT is a bit more complex on an image and you will have to learn how to apply it, but the principle is the same.


suppose if we send AM modulated message signal through space how the receiver will be knowing that we applied AM modulation?

Communications systems are always designed under the assumption that both emitter and receiver know what "language" they will be speaking to each other.

AM modulation comms are standardized in frequency, channel width for example. So each receiver is materially designed to demodulate AM signals with the appropriate hardware.

When people start trying to find which modulation was applied to a signal its generally involving electronic warfare or reverse engineering.


Why the correct probability for a square M-QAM is equal to the square of correct probability of ($\sqrt M$)-PAM?

The decision in M-QAM is actually two independent sub-decisions:

(i) what is the horizontal coordinate of the constellation point nearest to the received signal point?

and

(ii) what is the vertical coordinate of the constellation point nearest to the received signal point?

Consequently, the decision is correct if and only if both
sub-decisions are correct. But the probability of a sub-decision being correct is just $(1-P_{\sqrt M})$, and so the probability that both
are correct is $(1-P_{\sqrt M})^2$ by independence of the events.  If
the I and Q branches of the receiver had mismatched gains, the independence would still hold but the two probabilities being multiplied
would be different (corresponding to the different SNRs in the two
branches).


Comparing multiple signals for similarity

...best results come from a weighted ensemble of techniques...

Maybe they do, depending on the application. But each one of the similarities mentioned, is equivalent to the other at least when we are talking about signals originating from linear systems.

Cross correlation provides very good results especially if you are trying to figure out if a signal is noise or not and a very good application of that is in the rake receiver. And via the correlation's link to the Discrete Fourier Transform (DFT) that enables us to carry it out in the frequency domain we can start to see the equivalences between the given statements.

A more informative division would be between linear "metrics" and non-linear metrics. But, the field of "similarity metrics" is vast, so, it is best to start with what is known from the problem and then work towards a solution.

What I am going to do here is to talk a little bit about some metrics of similarity which you might find useful but then also write a little bit about other ways of doing the same thing that don't require the assessment of similarity as I think that that second part is closer to what you are trying to achieve (but not what you are asking directly).

On Similarity Metrics:

Cross-correlation will work fine if your signals have been consistently sampled to each other. That is, each device's sampling frequency remains absolutely stable but the relative points in time that the recording started can be varied. This will insert a delay between waveforms which is very easy to recover by doing multiple cross-correlations at different delays.

If you want to look at specific frequencies, you can either filter the signals prior to cross-correlation or use coherence.

Another thing you can try is Principal Component Analysis. To do this, pack all the signals that you receive from $N$ devices into a matrix and send it to PCA. PCA will try to decompose these signals into a smaller set of "principal" (uncorrelated) signals and will return a matrix with those elementary signals and their strengths. It does this using covariance, which is very closely related to correlation. The point here is that if all of the signals are the same, there are no "components" and the vector of strengths is going to be zero. But, if there is significant variation within the signals, the vector of strengths will start exhibiting some structure. You can then use the vector of strengths in a simple classifier that discriminates whether the indicated variation is negligible or worthy of an alarm.

Now, it might be that the sampling period of each GPS device is not stable (very likely). And this now means that if you were to put "one signal on top of another" (as if they were printed on transparency slides), it would be possible to align some peaks, but not others.

In that case, you could use Dynamic Time Warping (DTW). Imagine that you listen to a vinyl record with a lot of wow. You can remove that wow, by speeding up the recording when the record player slows down and vice versa. This is exactly what DTW does. You pass to it two signals and it returns a matrix that  tells you where to slow down and where to speed up. After you apply this transformation to one of the signals, the peaks would align perfectly. The total amount of deviation in this matrix is a metric of similarity. But, DTW is not exactly fast...particularly if you have to consider $\frac{N*(N-1)}{2}$ unique combinations (because all of these metrics are symmetric) between all 2 out of $N$ signals.

On the specific application:

So, the application is a GPS interference alert. There are two use-case scenarios that I can see at the moment:

GPS Jamming

Where the receiver cannot access the signals from the satelites because a local source is "flooding" its input with "noise".

GPS Spoofing

Where the receiver thinks they are listening to a satellite, but actually they are listening to a loud black box that behaves like a satelite.

Since it is impossible to raise an interference alarm based solely on information available to a device, the idea is to incorporate information from other data sources (nearby devices / data) in order to confirm the "story" that the GPS receiver on this device is listening to.

How can we get a gauge on the quality of the fix on a given device?

There are two obvious choices:

The Dilution of Precision
The GPS Quality Indicator Field.

Between these two, the Dilution of Precision is a "softer" indicator (wider range of values) of the quality of the fix.

So, assuming that there is an infrastructure for collecting GPS data from a set of devices (like Waze or Glypmse do for example), then the following alarms can be raised:

To counteract GPS Jamming:

Given the timestamp and the GPS Quality Indicator for a set of devices

Select a point in time
Select the subset of devices which have provided "bad" GPS quality data at or around a reasonable time interval around the selected point in time.
Cluster the provided GPS locations in space

This can be done ridiculously fast
A Threshold is required to define which two clusters are considered apart enough. This can be discovered automatically or be set to a reasonable value. Bunches of receivers approximately 1km apart should be considered separate, for example.

Assess the "size" of the "patch" of bad reception.

An obvious choice here would be to take the centroid of the cluster and find the furthest data point (so, one radius), or find the mean radius of the points or take their convex hull or use another more accurate algorithm that describes the points.

Decide if the patch is big enough to warrant raising the alarm.

Obviously, lots of receivers in a short tunnel or underpass or close to tall buildings are not supposed to raise alarms.

To counteract GPS Spoofing:

This is much more challenging because the problem there is to identify an unjustifiable deviation. This assumes that you have some way of knowing what is "correct" and that is the key difficulty. The easiest thing to tackle here is the GPS being spoofed over a wide area (more than one device).
In this case, you want to identify a cluster of devices that start behaving consistently all together. As you can understand, if 10 sailing autopilots were listening to the same spoofed signal and we start rotating the frame of reference then all of them will be taking consistent decisions at the same time. The above algorithm more or less remains the same but this time you want to cluster devices whose bearing (derived bearing) deviated by a consistent value at the same time. That is, they all appear to be adjusting at the same time and by the same amount. (As you can see here, no reason to assess the similarity of specific signals, we work with location directly).
If you are trying to understand if a single device's GPS is being spoofed, then that is even more challenging and there you would need to have access to other data as well (again to try to confirm the story the GPS is telling you).

Hope this helps


Standard deviation of spectral flatness -- so what am I measuring, conceptually?

Since you are interested in the "flatness" of your spectrum, in fact, you are interested in how close your signal is to a white noise (that by definition has a flat spectrum + random phases). If you step back, one measure would be the "distance" of your observation to the white noise reference.

The obvious measure in terms of information theory is Kullback-Leibler divergence. You do not need to understand every part of it, but it measures in bits (if you use log base 2) the distance between both distributions.

The good thing in your case is that your reference is flat, so that what remains is the entropy of your spectrum. There are many existing implementation (for instance in scipy).

Note that you are still on the safe side: if your distribution is approximately gaussian, both measures (entropy and std) will be proportional. The entropy is however more general and more principled. As an extension, you will be able to generalize to other types of noises (1/f for instance).


Frequency domain samples to time domain without knowing the sampling rate?

Generally the relationship between the sampling frequency $f_s$ and the frequency spacing of each bin $\Delta f$ is given by:

$$\Delta f= f_s/N$$

For example if you have 1000 bins and the sampling rate is $f_s = 1$ KHz, then each bin is spaced by 1 Hz given by $f_s/N$. So if the frequencies $f_1, f_2, ..., f_N$ (using the OP's indexing) were associated with bins 0 to $N-1$, each spaced by 1 Hz. In this case the frequencies would be $0, 1, 2... 999$ Hz and the sampling rate is $999$ Hz + 1 bin or $999$ + $1$ Hz = 1 KHz.

Update due to additions to the question:
In the specific case that the OP is using, the center unused bin would most likely represent the carrier frequency which would map to DC (and it is unused due to the possibility of interference due to carrier feed through). So here the center is likely to be bin 0, the frequencies to the right of center are likely to be bins 1 to 26, and the frequencies to the left of center are likely to be bins 27 to 52 (bin 27 is to the far left, representing the negative frequencies as would be mapped with fftshift in MATLAB/Octave). Regardless due to the cyclical nature of the FFT we can still use the guidance I provided to determine the sampling rate such that it would be N=53 total bins (including the nulled bin) and with a frequency spacing of 312.5 KHz the sampling rate would be 312.5 KHz x 53= 16.5625 MHz.


Determining which Filter from a Z-Plane Plots?

just to add a little to Matt's answer.  if you know where your poles, $p_n$, and zeros, $q_m$, are your transfer function is

$$ H(z) = H_0 \frac{(z-q_1)(z-q_2)(z-q_3)...(z-q_M)}{(z-p_1)(z-p_2)(z-p_3)...(z-p_N)} \quad \quad M \le N $$

the complex frequency response (which contains both magnitude and phase information) is

$$ H(e^{j\omega}) = H_0 \frac{(e^{j\omega}-q_1)(e^{j\omega}-q_2)(e^{j\omega}-q_3)...(e^{j\omega}-q_M)}{(e^{j\omega}-p_1)(e^{j\omega}-p_2)(e^{j\omega}-p_3)...(e^{j\omega}-p_N)} $$

and the magnitude response (which is where you determine if it's a LPF or BPF or HPF or something else like that) is:

$$ |H(e^{j\omega})| = |H_0| \frac{|e^{j\omega}-q_1| \ |e^{j\omega}-q_2| \ |e^{j\omega}-q_3|...|e^{j\omega}-q_M|}{|e^{j\omega}-p_1| \ |e^{j\omega}-p_2| \ |e^{j\omega}-p_3|...|e^{j\omega}-p_N|} $$

and $e^{j\omega}$ is a point on the unit circle.  what this means is, as your frequency increases from DC (or $\omega=0$), that corresponds to moving along the unit circle from $z=e^{j0}=1$ to other points on the unit circle.  then as your point on the unit circle gets close to a zero $q_m$, a factor of $|e^{j\omega}-q_m|$ gets small.  $|e^{j\omega}-q_m|$ is literally the distance your point on the unit circle is from that particular zero.

likewise for the poles, but your poles are in the denominator, so as your point on the unit circle gets close to a pole $p_n$, a factor of $|e^{j\omega}-p_n|$ gets small, but since you're dividing by that small number, the effect on the gain magnitude is to increase it.

so a LPF will have poles closer to the DC point of $z=e^{j0}=1$ and, perhaps, zeros closer to Nyquist $z=e^{j\pi}=-1$ so when you're close to DC you're multiplying by big numbers and when you're close to Nyquist, you're multiplying by small numbers.

likewise for a HPF, you'll have poles closer to Nyquist and zeros closer to DC.

for a BPF, you might have zeros close to both DC and Nyquist and poles close to the unit circle at other places where they will boost the gain for those bandpass frequencies.

that's how you look at the pole-zero plots on the $z$-plane to judge if it's a LPF or whatever.


Why does fft require a power of 2 sample numbers?

It doesn't.  An FFT is just a method or implementation of computing a DFT that runs faster if the length (N) can be factored.  It runs relatively faster (compared to a DFT of the same length) given a larger number of small prime factors of N, and is simpler for smaller numbers of unique prime factors of N (such as: only one).

Picking just one prime factor, such as the smallest, 2, makes textbook derivations and examples of FFT computations use less ink in the textbook and/or less chalk on the chalkboard.

But many modern FFT libraries are capable of computing an FFT for an N that is a composite product of several other (fairly small) prime factors, not just 2.


so I am wondering where can I go to form an understanding of DSP from the ground up.

There is a great list of resources in this answer What Resources Are Recommended for an Introduction to Signal Processing (DSP)?  You should spend half an hour with each one and then decide which fits your own learning style the best.

My skills in calculus are also quite bad.

That's a problem. DSP is by nature fairly math heavy and if you lacking some fundamental skills, you will we miserable the entire time. It's like trying to build a house but all you have is a hacksaw and a bent hammer. Get the right tools first. At the very least make sure you are familiar with complex number and the associated math. There are plenty of resources to do that, including Khan Academy.


The structure tensor consists of first derivatives of the image. If the first derivative is high in one particular direction (one large eigenvalue of the structure tensor), then you have an edge. It the first derivative is high in two directions (two large eigenvalues), then you have a corner.

The Hessian matrix consists of second derivatives. Think of the Laplacian. It is the trace of the Hessian, and it is a center-surround filter, which gives you a high absolute value for a bright blob surrounded by dark background or vice versa. That's why it is blob detector.


Write $h(z) = -1 + 3z -2z^2 + z^3$ and compute $h(z) \bmod (z^2 - 1)$, that is,
divide $h(z)$ by $z^2 - 1$ and take just the remainder.  While this looks
very complicated, if you think about it a bit you will see that all you
are doing is dividing $[-1\quad3\ -2\quad 1]$ into $[-1\quad 3]$ and
$[-2\quad 1]$ and
adding the shorter vectors to get $[-3\quad 4]$.  Repeat for
$x = [1 \ -1\ -2\quad 1 \quad 3 \quad 2\quad 1 \quad 2]$ to add
four vectors of length $2$ to get
$[3\quad 4]$.  This is presumably not too hard to do in MATLAB though I am
not familiar enough with the syntax to suggest
specific commands.  Next, compute the cyclic convolution
of $[-3\quad 4]$ and $[3\quad 4]$ preferably without invoking MATLAB functions.
The result is
$$[\{(-3)\times 3 + 4\times 4\}\quad \{(-3)\times 4 + 4 \times 3\}]
= [7\quad 0]$$

Mathematically, what you are doing is computing $h(z)x(z) \bmod (z^2-1)$
which can be done fancily by first finding $h(z)x(z)$ using FFTs and what
have you followed by the $\bmod (z^2-1)$ computation (this effectively
chops up the long vector into short pieces and adds them), or more simply
by first computing $\hat{h}(z) = h(z) \bmod (z^2-1)$ and
$\hat{x}(z) = x(z) \bmod (z^2-1)$ (the chopping up into shorter vectors
and adding them) and then computing the cyclic convolution
$\hat{z}(z)\hat{x}(z) \bmod (z^2 - 1)$ which is easy to do.

  Chop-add-convolve is easier than convolve-chop-add


I think that the question here is "how is the dominant direction of the local slope actually estimated?" with the outlook of using it in a shock filter.

A Shock Filter is applied iteratively and each time it propagates grayscale values along the direction of an edge but not across the edge, therefore, progressively, the edge is preserved but the area towards the broad direction of the structure is "filled". This is similar to applying dilation or erosion but taking into account the local contrast gradient. The edge can be defined at various scale levels, which results to the average direction of an edge defined over a bigger area surrounding a pixel of interest $(x,y)$.

The code seems alright, in that it sets up the normalised Gaussian g and then derives its first order derivative in gd. It uses that to obtain the gradients Ix,Iy but sx, sxy should be pointwise, rather than matrix, operations. Otherwise, Ix, Iy, "filter" each other. By controlling for the variance parameter of g (and by extension, gd too), we control how big the area around the pixel of interest $(x,y)$ is (or, in terms of dilation and erosion, how big the disk -or structuring element- is).

At this point, Ix, Iy could be used to estimate the local contrast gradient but if we come across areas with roughly opposite slopes, then averaging their components would result to 0 and the (false) perception that the area is flat. Instead, the Structure Tensor is used and more specifically the eigenvector that corresponds to the largest eigenvalue (and there is one of those for each $(x,y)$ pair). For more information please see this link and specifically sections 2, 3 and 4.

So, with that vector indicating the direction of the gradient of an edge, we find the perpendicular vector (with that sign flipping) and propagate values in that direction.

Hope this helps.


Implementing Gaussian random variable by using a uniform random variable

Your initial algorithm creates a random variable that's uniformly distributed between 0 and 1. The variance of that is 1/12. If you sum NUM_GAUSSIAN_SUMS instances of that the variance will be NUM_GAUSSIAN_SUMS/12. In order to get to a target variance, V, you need to multiply the summed random variable with sqrt(V*12/NUM_GAUSSIAN_SUMS).

As a side note, a template will work reasonable well for floats and doubles but there will be significant numerical problems with any fixed point type.


Can the FFT reduce its computation load for specific sparse data patterns?

Sort of but only of it's very sparse. You can always NOT use an FFT but implement the inverse FFT directly and omit the 0 terms from the sum.
We can guestimate the break even point. An FFT takes about 10 operations per butterfly and for a length 256 FFT you have 8 stages with 128 butterflies each, so you end up with about 10000 operations per FFT
For a direct inverse DFT you'll have for each real output sample about 3 operations per non-zero frequency coefficient. In this case the break even point comes out to be 13 or thereabouts.
Your case of 16 is probably close to a wash, so I wouldn't bother.


Savitzky–Golay filter vs. IIR or FIR linear filter

Since the discussion in the existing answers and comments has mainly focused on what Savitzky-Golay filters actually are (which was very useful), I will try to add to the existing answers by providing some information on how to actually choose a smoothing filter, which is, to my understanding, what the question is actually about.

First of all, I'd like to repeat what has become clear in the discussion spawning from the other answers: the categorization of smoothing filters in the question into linear and time-invariant (LTI) FIR/IIR filters on the one hand, and Savitzky-Golay filters on the other hand is misleading. A Savitkzy-Golay filter is just a standard FIR filter designed according to a specific criterion (local polynomial approximation). So all filters mentioned in the question are LTI filters.

The remaining question is how to choose a smoothing filter. If computational complexity and/or memory are an issue, IIR filters may be preferable over FIR filters, because they typically will achieve comparable noise suppression (i.e., stopband attenuation) with a much lower filter order than FIR filters. But note that if real-time processing is necessary, one possible disadvantage of IIR filters is that they cannot have an exactly linear phase response. So the desired signal will suffer some phase distortions. For offline-processing, the phase distortions can be avoided, even with IIR filters, by applying zero-phase filtering.

Apart from the considerations discussed in the previous paragraph, it is mainly the design criterion that matters, not so much if the filter is FIR or IIR, because any (stable) IIR filter can be approximated with arbitrary accuracy by an FIR filter, and any FIR filter can be approximated by an IIR filter, even though the latter can be much more difficult. The appropriate design criterion obviously depends on the properties of the data and the noise. When it comes to smoothing we usually assume sufficiently oversampled (i.e., smooth) data. If the noise has mainly high frequency components, i.e., if there is little spectral overlap between the data and the noise, we want to maximize stop band attenuation, or minimize stop band energy, while preserving the desired signal as well as possible. In this case we could choose a linear-phase FIR filter designed according to a minimax criterion using the Parks-McClellan algorithm. We could also minimize the stop band energy (i.e., minimize noise power in the stop band) by choosing a least squares method. A mix between the two criteria (minimax and least squares) is possible by choosing a constrained least squares design, which minimizes stop band energy while constraining the maximum approximation error in the pass band.

If the noise spectrum significantly overlaps with the signal spectrum, a more careful approach is required, and brute-force attenuation will not work well because either you leave too much noise (by choosing the cut-off frequency too high) or you distort the desired signal too much. In this case Savitzky-Golay (S-G) filters may be a good choice. The price to be paid is mediocre stopband attenuation, but one advantage is that some signal properties are preserved very well. This has to do with the fact that S-G filters have a flat pass band response, i.e.,

$$\frac{d^kH(e^{j\omega})}{d\omega^k}{\huge |}_{\omega=0}=0\quad k=1,2,\ldots,r\tag{1}$$

where $r$ is the order of the approximating polynomial and $H(e^{j\omega})$ is the filter's frequency response. Property $(1)$ guarantees that the first $r$ moments of the input signal are preserved in the output, which means that the width and height of peaks in the desired signal are well preserved.

Of course there's also a compromise between the two types of smoothing filters discussed above (high stopband attenuation and S-G). We could design an FIR filter with a certain degree of flatness at $\omega=0$ and use the remaining degrees of freedom to maximize stop band attenuation, or minimize stopband energy. In the case of FIR filters, the resulting design problem is sufficiently simple (and convex), and general optimization routines available in several software packages can be used to obtain the optimal filter for the given application.

For the ones interested in the theory of S-G filters, the most relevant references I can recommend are the following:

the original paper by Savitzky and Golay
a tutorial paper by Ronald W. Schafer
section 8.3.5 of Orfanidis' DSP book


Remove unwanted object of image in MATLAB?

What you are asking is probably hopeless.

There are two quite difficult aspects:

detecting the trees that you want to remove. Given the very poor contrast, this will be hard to achieve by automated means. You need to implement some image segmentation method, then select the regions corresponding to the trees.
reconstructing the background image under the occluded parts. This is called image inpainting and requires sophisticated approaches.

All of this is highly technical and the standard tools in Matlab will probably not suffice.

If you just have one image to process, photoshopping will be much more efficient.


Envelope of time domain signal // Hilbert transform

Hilbert envelope, also called Energy-Time Curve (ETC), only works well for narrow-band fluctuations. Producing an analytic signal, of which you later take the absolute value, is a linear operation, so it treats all frequencies of your signal equally. If you give it a pure sine wave, it will indeed return to you a straight line. When you give it white noise however, you will likely get noise back.

Example 1:

$$
\lvert\cos(\omega t)\rvert^2 + \lvert i\sin(\omega t)\rvert^2  = \cos^2(\omega t) + \sin^2(\omega t)= 1
$$

Here, you have $\cos(\omega t)$ as you input signal, and its analysis counterpart is $\cos(\omega t) + i\sin(\omega t)$. In it easy to see from expression above that its envelope is always a constant. Let's look at another example.

Example 2:

Let $x(t) = \cos(\omega_1 t) + \cos(\omega_2 t)$. Then its analytic counterpart is
$$
a(t) = \cos(\omega_1 t) + \cos(\omega_2 t) + i\sin(\omega_1 t) + i\sin(\omega_2 t).
$$
Let's take the absolute value of this signal.

$$\begin{align}
\lvert a(t)\rvert^2 &= \lvert\cos(\omega_1 t) + \cos(\omega_2 t)\rvert^2 + \lvert i\sin(\omega_1 t) + i\sin(\omega_2 t)\rvert^2\\&=\left(\cos(\omega_1 t) + \cos(\omega_2 t)\right)^2 + \left( \sin(\omega_1 t) + \sin(\omega_2 t) \right)^2 \\ &= \cos^2(\omega_1t) + \cos(\omega_1t)\cos(\omega_2t) + \cos^2(\omega_2t)\\& \quad + \sin^2(\omega_1t) + \sin(\omega_1t)\sin(\omega_2t) + \sin^2(\omega_2t) \\&= \cos^2(\omega_1t) + \frac{1}{2}\left[ \cos(\omega_1t-\omega_2t)+\cos(\omega_1t+\omega_2t)\right] + \cos^2(\omega_2t)\\&\quad + \sin^2(\omega_1t) + \frac{1}{2}\left[ \sin(\omega_1t-\omega_2t)-\sin(\omega_1t+\omega_2t)\right] + \sin^2(\omega_2t) \\ &= \left[ \cos^2(\omega_1t) + \sin^2(\omega_1t)\right] + \left[ \cos^2(\omega_2t) + \sin^2(\omega_2t)\right] \\ &\quad + \frac{1}{2}\left[\cos(\omega_1t-\omega_2t) + \cos(\omega_1t-\omega_2t) \right] + \frac{1}{2}\left[\sin(\omega_1t+\omega_2t) - \sin(\omega_1t+\omega_2t) \right] \\ &= 1 + 1 + \cos(\omega_1t - \omega_2t) + 0 \\ &= 2+\cos\left( \left[\omega_1-\omega_2\right]t\right)\end{align}$$

As you can see, this is not a constant as one would expect. Note that if $\omega_1$ and $\omega_2$ are close to each other, only low-frequency oscillations are produced as an artifact. hope this helps.


How to generate a sinusoid with STFT at a frequency which is between bins

A sinusoid with a frequency that is between bins in the FFT frequency domain is circularly discontinuous in the time domain.  So you can't use the same IFFT results back-to-back without the noise from this discontinuity between each IFFT window, as the end of one window will have a value too far from the beginning of the next window (unless your frequency is bin centered, e.g. exactly integer periodic in the FFT's length).

If you calculate and vary the starting phase for each successive (complex input) IFFT, then you might be able to reduce or get rid of noise from the window discontinuity.

Note that the spectrum of a "between bin" frequency sinusoid is not just a cross-fade between the two nearest bins, but a Sinc shaped spectrum, spanning the full width of the FFT.  (Or a windowed Sinc function for a decent approximation.)


I have:

The function $h_1(t)=\cos(2\pi\frac14 t)$ sampled by $|||(t)$ has fourier transform $H_1(f)=4||(4f)*|||(f)=2|||(2(f-1/4))$.
The function $h_2(t)=\cos(2\pi\frac12 t)$ sampled by $|||(t)$ has fourier transform $H_2(f)=2||(2f)*|||(f)=|||(f-1/2)$,

When I bandlimit $h_1$ with $\Pi(2f)$ I obtain $4||(4f)$.
If I try to bandlimit $h_2$ by the same window, exploiting the Heaviside half maximum convention $\Pi(\pm1/2)=1/2$, I obtain $2||(2f)$.

Both sinuoids are perfectly recovered, though the Heaviside limit in the second case is not applicable in practice for dozens of reasons.

In here, I am using the pulse, sampling and sinusoidal symbols. $H$ is the Heaviside unit step function:
$$
\Pi(f)=H(f+1/2)-H(f-1/2),\\
|||(t)=\sum_{t_n \in \mathbb{Z}}\delta(t-t_n)\\
||(t)=\frac12\delta(t+1)+\frac12\delta(t-1)
$$


relation between luminance and eigen values of an image

Focusing too much on the eigenvalues neglects the eigenvectors which are more important here.

In a way, it is natural to associate eigenvalues with brightness but what is actually achieved with such techniques is contrast enhancement. This is not a mere scaling of the magnitude of the pixels.

For example, given an image $A$, its pixel brightness can be changed by $2 \times A$. $A$ will now appear brighter but the relationships between the values of the pixels will have been preserved. Therefore, any objects hidden in shaddows (small pixel differences relative to their background) will remain hidden in the shadows.

To see those objects, you would have to change the contrast, that is, the relationship BETWEEN pixel values. The most simple way to achieve this is with an image window (nothing to do with spectral windows). For more information about this please see this and this link. There are other ways to modify contrast, often equivalent to a complex window, such as equalisation, normalisation and others but these are not the only ways to modify contrast.

We usually tend to think of some matrix $A$ as a whole (e.g. an image). However, for the purposes of eigendecomposition, it is more useful to look at some $A$ with dimensions $M \times N$ as a way of grouping together a set of vectors from a multidimensional space.

So, in our mind, we now abandon the view of $A$ as a spatial index of intensity (which is what an image is) and we start thinking about a huge "box" into which we see lots and lots of little vectors that are pointing towards all sorts of directions.

Now, imagine that you take a completely new vector, of equal dimensions and you stick it in this space at some orientation. You can then estimate how many vectors are aligned with this direction and make a note of it that looks like (Direction, Number of vectors aligned with that direction). Change the direction of your "controlled vector" and do it again, do it again and again and again and again and again until you have a (very) long list of (Direction, Magnitude) measurements.

Obviously, a non random "vector cloud" such as the one we are imagining here will have some sort of orientation. That is, most of the vectors will be pointing mostly (not necessarily EXACTLY) towards some specific direction. This direction will have very high "counts" (of aligned vectors) and as the direction of your "controlled" vector moves away from this your counts will be decreasing rapidly.

Let's get back to the eigendecomposition now. Eigendecomposition is a way of discovering these pairs of (Direction, Magnitude) WITHOUT having to stick a "controlled" vector at different orientations and count alignments. In the end of this process, you get pairs of (Eigenvector (a.k.a Direction), Eigenvalue (a.k.a Magnitude)).

Just as a side note, if you can decompose the "vector cloud" into (Eigenvector, Eigenvalue) pairs then you can recompose it from these pairs too.

Back to our original problem: "So what? I don't have vectors, I have an image, how does playing around with little arrows affects the image?".

Let's say that your "vector cloud" assumes a particular form (a shape) in this multidimensional space. That's not hard to do because an image that is not noise does show something so the pixels do have some sort of order and are not all over the place. BUT, the vectors in this "vector cloud" might be corrupted by some random fluctuations in their directions (a.k.a noise). So, in this case, we can take advantage of the fact that a very large number of vectors are pointing in some general direction which we can recover. All we have to do now is RECOMPOSE the vector cloud keeping only a subset of the stronger (Eigenvector, Eigenvalue) pairs.

This, is the main idea. But it is not always applied directly on a given image.

For example, in the case of eigenfaces, the eigendecomposition is applied on the covariance matrix between images. This provides the general direction of a "vector cloud" that describes covariances between images. In effect, eigenfaces are a small set of elementary images that describe a very large set of faces as linear combinations.

In the first paper that you are linking in your original post, the eigendecomposition is applied on the Discrete Cosine Transform (DCT) field of the actual satelite images, again, the main idea is the same but here authors are discovering a set of elementary spatial frequency components.

For more information, please see this and this link.

Let us now turn briefly to your experiments. It is probably obvious by now that by focusing on the eigenvalues you are obtaining a rough approximation of the scene's brightness. Furthermore, the logarithmic dependence that you see in your graph is perfectly normal and due to the characteristic curve of the imaging device you are using (e.g. a digital camera at different f-number settings.)

If all you are trying to do is...

  study the variation of eigen values with respect to luminance of an
  image

...then there will be some dependence but it will be indirect. As you are modifying the exposure of a given scene, you will pass from a set of obvious points.

Point-1: Image severely underexposued, in this case most of the
image is dark, it has very low contrast and is associated with very
low eigenvalues and very few eigenvectors...There is not much to
describe here.
Point-2: Image perfectly exposed, in this case you are getting the
maximum of contrast that can be obtained from a given scene with a
given imaging device. You will get higher number of eigevectors and
eigenvalues than point-1 because every pixel of your image now
contains some information about the scene being imaged.
Point-3: Image severely over-exposed, in this case you are again
operating at a region of low contrast but this time it is because
the sensor is saturated. Any sort of difference in luminus flux goes
unnoticed by the sensor which is simply recording the maximum of
value it can produce at its output. Consequently, you get a mostly
washed out image which again has a lower number of eigenvectors (but
possibly high eigenvalues).

Hope this helps.


Difference equation when transfer function expressed as poles and zeros

If I understood your question right, you would like to obtain a Linear Constant Coefficient Difference Equation (LCCDE) representation of a given system from its given System Transfer Function, $H(z)$, when it's expressed in a pole-zero product form such as:
$$
H(z) = g \frac {(z - d_0)(z - d_1)} {(z - c_0)(z - c_1)}
\tag{1}
$$

instead of a, more directly apparent, power series form such as:
$$
H(z) = \frac {b_0 + b_1 z^{-1} + b_2 z^{-2}} {1 + a_1 z^{-1} + a_2 z^{-2}}
\tag{2}
$$

As known, the latter form can be directly converted to its corresponding LCCDE representation from the coefficients ${b_k}$ and ${a_k}$ as follows:
$$y[n] + a_1 y[n-1] + a_2 y[n-2] = b_0 x[n] + b_1 x[n-1] + b_2 x[n-2]$$
Where the pattern of conversion speaks for itself, without a need of further description.

Now, we'd like to know whether one can obtain a similar conversion from the pole-zero form into the LCCDE without expanding the product?

Well yes and no! Depending on how you consider the process of computing the coefficents. If your intention is on the fact that you do not want to algebricaly perform the paranthesis expansion operation, but instead use a machanism that computes each coefficient separately by a set of arithmetic operations; then yes you can devise such a method to, for example, compute $a_3$ or $b_5$ without computing any other terms. And No, eventually you must always compute the coefficient $a_k$ that multiplies $z^{-k}$ so as to find the coefficient that multiplies $y[n-k]$, hence you need to always expand the brackets and compute multipliers of $z^{-k}$ in one or the other way.

Then I want to present here a simple method to find those coefficients :
given a polynomial in the product of zeros form: $$p(x) = (x-a)(x-b)(x-c)$$ we can perform the following operation to get its coefficients in the expanded form $$p(x) = d_0 x^3 + d_1 x^2 + d_2 x + d_3$$

$$d[k] = (1 -a)\star(1 -b)\star(1-c)$$ where $k =0,1,2,3$ in this particular case.

As clearly seen this is a discrete convolution operation, where each multiplier bracket is turned into a respective convolution operand.  For which a simple MATLAB line to compute $d[k]$ would be:

 d = conv([1 -a],conv([1 -b],[1 -c]));

where a, b, and c would be replaced by their numeric values or if you're competent in the symbolic math capability, you could instead obtain a pure symbolic result for the coeffients $d_k$, which would therefore give you a way to compute LCCDE coefficients from pole-zero products of $H(z)$ in a rather roundabout way.

Finally, the computations for a power of $z^{-1}$ products is basicly the same if you describe the product as $$P(z) = (1 - az^{-1})(1-bz^{-1})(1-cz^{-1})$$ and,
$$P(z) = d_0 + d_1 z^{-1} + d_2 z^{-2} + d_3 z^{-3} $$

where again $$d[k] = (1 -a)\star(1 -b)\star(1-c)$$


Why is there no FPU on (most) DSP chips?

An 8086 used less than 30k transistors. The 8087, which is the FPU for the 8086, is reported to use 45k transistors.  Faster FPUs can be even larger in terms of gate count.  So the cost in silicon die area of an FPU can be significant (over 2X?).  Power and thus heat is proportional on the order of the number of transistors toggling outputs at similar rates.

For real-time DSP, there is also the issue of deterministic latency. A simple integer multiply-accumulate unit has a fixed latency.  A really simple floating point adder does not have a fixed latency (or is either not fast or not simple) due to the potential need to normalize up to twice and/or handle NaNs and denorms.


How to compute number of samples in 1 cycle of signal?

Although it is hardly related to DSP but here is a quick calculation:

$R$ Rotations per minute is $\frac{R}{60}$ Rotations per second which gives $12K$  samples per second. Hence, One rotation is proportional to $\frac{1}{\frac{R}{60}}\times 12000 = \frac{720,000}{R}$ samples. For instance, for  $2000$ RPM there are $360$ samples per rotation.


Difference between Analog IQ sampling and creating IQ digitally

You're right.  Whenever you can do things digitally you should because of all the advantages that you've listed.  But a long time ago it was expensive to do those low pass filters digitally.  If the LPF required is narrow (and requires a steep roll off), it would require many taps (multipliers) and they might require many bits so as to not lose precision to numerical accuracy issues. You also have to have the ADC twice as fast in the all digital case.
The disadvantages of doing the things in the analog realm are many though.  You would have to carefully tune LPFs to have the same group delay.  There could be many calibration issues.
Here's a paper describing a Hilbert Transform design


The distribution of the log magnitude output of a log amplifier would be according to the log of the magnitude of the input. Although in implementation the log amp is not simply taking a log of the input signal, the result is indeed the log with reasonable accuracy (typically linear to the log scale within a fraction of a dB). Below shows how this should appear for both Rayleigh and Ricean distributions. Given the OP is not seeing the same skewness as is the case for the Rayleigh at the log magnitude output of the log amp, I would conclude that the input signal (which he does not have access to) was not Rayleigh distributed.

Notice the skewness decreases significantly as the dominant path in the Ricean distribution (A) becomes more prominent, as would be expected.


Integer upsampling - leave original samples unaltered

You need to be careful with the delay compensation. If your symmetric filter has an even number of taps, then you get a non-integer delay. If you use a symmetric filter with an odd number of coefficients then the delay of $\frac{N-1}{2}$ samples is an integer number ($N$ is the number of filter coefficients).


What is the rate of change of amplitude? Frequency?

Fourier transforms and differentiation are indeed very different, but your intuition is still not entirely wrong. Here's how to make sense of it.

Let's assume you have a signal that is essentially sine-shaped locally, i.e. we can write it as $s(t) = a \sin(2 \pi f t)$ for a times $t$ close to some time $t_0$. In other words: The signal oscillates like a sine and has a sufficiently slowly changing envelope and instantaneous frequency.

The time derivative of this signal is $s'(t) = 2 a \pi f \cos(2 \pi f t)$. This signal is also oscillating with the same frequency and its amplitude is scaled by $2\pi f$ compared to the original signal. So the frequency is encoded in the derivative.

There are a few different methods for decoding the frequency from this derivative, and I'll describe the one that is probably closest to your intuition. We have to get rid of the oscillating part. So let's take first the magnitude of both signals and afterwards apply a low pass filter $H_l$ that removes the residual oscillation and leaves you with just the average value:

$$H_l(|s(t)|) = a \frac{2}{\pi}$$
$$H_l(|s'(t)|) = 2 a \pi f \frac{2}{\pi}$$

( The 2/$\pi$  comes from averaging over $|\sin(t)|$. ) With this, we can find the quotient of the two results and have

$$\frac{H_l(|s'(t)|)}{H_l(|s(t)|)}=2\pi f$$

which demonstrates that the derivative of the signal in fact yields a very accurate approximation of the frequency under the assumptions we made.

For more general signals this method calculates not the instantaneous frequency but the so called spectral centroid frequency.


Figured out how to correctly change from the zero-poles coefficients to the numerator/denominator coefficients.

Essentially this simply requires a convolution of the negative of each value with 1 being added as the first value in the results vector.

The implementation of this for the application of convolution in the instance of filter design is:

var result = new List<Complex>
{
    1
};

for (var i = 0; i < roots.Length; i++)
{
    var value = -roots[i];
    result.Add(result[result.Count - 1] * value);
    for (var j = result.Count - 2; j >= 1; j--)
    {
        result[j] = result[j] + value * result[j - 1];
    }
 }

To visualise this we start with the roots as being:
[0.954697396423354, 0.8855894564017569, 0.86046511627906985, 0.8855894564017569, 0.954697396423354]

Before iteration:

result = [1]

First iteration:

result = [1, -0.9546974]
result[length-1] = 1 * -0.954697396423354

Second iteration:

result = [1, -1.84028685  0.84546995]
result[length-1] = -0.9546974  * -0.8855894564017569
result[length-2] = -0.9546974 + -0.8855894564017569 * 1

Third iteration:

result = [1, -2.70075197  2.42897259 -0.7274974]
result[length-1] = 0.84546995 * 0.86046511627906985
result[length-2] = 0.84546995 + 0.86046511627906985 * -1.84028685
result[length-3] = -1.84028685 + 0.86046511627906985 * 1

And so continues until you reach the end of the input vector.

The main thing to note is that the initial value (1) must never change in the result vector so it is required in the inner loop to make sure you skip the last value, as you've just calculated and added it, and also the first value.


in my opinion and limited experience, you can't just compare spectra of male and female voices to identify which is which.
2 decades ago i worked on two products called PurePitch and PitchDoctor that was able to shift the pitch of a voice independently of shifting the formants.  it also had the ability to increase or reduce the amount of pitch inflection of the voice (all the way to monotone).  the products had factory-supplied "presets", one was called "testosterone" and the other was called "estrogen", that could convincingly make a female voice sound masculine or a male voice sound feminine.
anyway, i don't remember the numerical parameters that were researched literally by my boss in usage of the alpha version of the product.  but it was something like 300 or 400 cents difference between the male and female voices and less than that for the difference between formants.
i think you might need a pitch detector and get lotsa words and do some statistics on the pitch to get a good machine idea of whether the pitch is or was more likely male or female.


Why Does 2D FFT of Gaussian Looks More Sharper than Gaussian Itself?

Have a look at the Fourier Transfrom of a Gaussian Signal.

$$ \mathcal{F}_{x} \left\{ {e}^{-a {x}^{2} } \right\} \left( \omega \right) = \sqrt{\frac{\pi}{a}} {e}^{- {\pi}^{2} \frac{ {\omega}^{2} }{a} } $$

First, Gaussian Signal stays Gaussian under Fourier Transform.

As you can see, the parameter which multiplies the variable is inverted.
Let's say $ a = 5 $, then it means that in time we will have very sharp and thin Gaussian while in frequency we will have very smooth and wide Gaussian.

This is related to other property of Forier Transform.
In simple words, what's thin on Time / Spatial Domain is Wide on Frequency Doman and vise versa.
This is basically the Uncertainty Principle of Fourier Transform.


Digital implementation of type I compensation network

The integrator is unstable, so it will probably just blow up. There's such a thing as a leaky integrator that does an OK job of getting pretty close. Have a look at the algorithm listed here on how you would go about implementing it using code. The math is pretty straightforward.

The differentiator uses b1=1, b2=-1. The integrator is just the inverse, but what you'll find is you need to modify (decrease) the filter coefficient by a subtle amount to make it 'leak'. For the integrator use something like a1=1, a2=-0.99. The closer a2 gets to 1.0 the greater the risk you run on it going unstable. Decreasing a2 to something like -0.95 will result in something that looks less and less like an integrator. Play around with it.


How to filter out noise from non-stationary signal

Break your signal into short overlapping frames over which it is stationary. Say, the length of the frame is $L$ samples. You can window it to smooth out discontinuities. For example, you can use a Hann window with an overlap of $L/2$ samples.
On each windowed frame, use the Wiener filter to find the optimal inverse filter coefficients to cancel your noise. See scipy.signal.wiener.
Overlap add your filtered signal frames to recover the original non-stationary signal.


To increase the speed, I think the dimensionality of the feature descriptors should be reduced. SIFT feature descriptor is of 128 dimensions.

So you can try PCA-SIFT which reduces the descriptor number to 20. It uses Principal Component Analysis to the SIFT descriptors.

You can also try using Speeded Up Robust Features (SURF). The feature descriptors obtained from SURF are invariant to rotation and scale, and it can reduce the dimensions of feature descriptors significantly. it enables faster feature detection as well as template matching. Also, MATLAB has a specific command for SURF i.e. detectSURFFeatures, so if you are using MATLAB for your job you are in luck. I think OpenCV also has the SURF feature command so that you can use it on python too. I have used the SURF features which I think is pretty fast.