best window for close frequency components

I'm presuming that you have some or all 9 of the input frequencies available at one time, and you just want to find out which of them is present (you didn't specify in your problem statement).

Two methods come to mind: 1) compute DFTs at the frequencies of interest, or 2) use the Gaussian ratio technique, as described in reference_1 and reference _2.

Regarding 1) above: When Moses came down from the mountain, he did not carry tablets of stone on which it was commanded that in the equation: F = k*(df) = k*sample_rate/N, that k has to be an integer.  In other words, you can compute a DFT bin centered at ANY frequency, not just those that correspond to the usual k = integer of your FFT.  For example, given your sample rate and N, a frequency of F = 16351 would correspond to a DFT bin computed at k = (N*F)/sample_rate = 1024*16351/44100 = 379.669 (ie: between bins 379 and 380 of your FFT).

So for the DFT, X(F) = (sum over n){x[n]*e[-j*twopi*n*k/N]}, where X[F] is the complex result, N = 1024, F = 16351, k = 379.669, and n = 0 to N-1.

Some of the other values for your problem are:

Frequency   and value of   k

16,566  ==> 384.662

16,781  ==> 389.659

16,996  ==> 394.646

I'll leave it to you to figure out the rest.  And, as you can see, you've got about 5 bins of separation between the input frequencies.  That may, or may not be enough based on your SNR and method used.  And you might want to use N = 1025 instead – it makes your input frequencies much more 'bin-centered.'

Method 2) above uses an FFT, which, given the right circumstances, can be extraordinarily accurate when computing frequencies.  Note that in the first reference above, by using a Gaussian window, and then taking the log ratio of amplitudes from adjacent bins, you are actually solving for frequency, subject to the requirements of the technique.  Some code (DevC++ compiler) for your specific problem, followed by the programs' output is shown below (explanation to follow):

#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <cmath>
using namespace std;
void fft_recur(long, double* r, double* i);
int main (int nNumberofArgs, char* pszArgs[ ] ) { // arguments needed for Dev C++ I/O
const long N = 1024;   double sample_rate = 44100., r[N] = {0.}, i[N] = {0.}, w[N] = {0.} ;
long n;      double  t, twopi = 2*3.141592653589793238;
double amp_sq[N], F, fbin, amp, power, c = 2.3;

for (n = 0; n < N; n++)  {  // generate test data
    t = twopi*n/sample_rate;
    r[n] = 1.*cos(16351.*t)+1.*cos(16996.*t)+1.*cos(17211.*t)+1.*cos(17426.*t) ;
} // end for

//Gaussian window as computed in [MCEACH94] program
double z;
for (n = 0; n < N/2; n++) {
    z = (.5+(N/2)-1.-n)*(twopi/2.)/N;
    w[n] = exp(-c*z*z);
    w[N-1-n] = w[n];
}

//multiply input data by window points
for (n = 0; n < N; n++)  {
    r[n] = w[n]*r[n];    // we're overwriting r[n] with the windowed r[n]
} //end for

fft_recur(N, r, i);       // note: fft is not scaled

// compute amplitude squared
for (n = 0; n < N/2+1; n++)  {
    amp_sq[n] = r[n]*r[n]+i[n]*i[n];
} //end for

// compute natural log of FFT outputs
for (n = 0; n < N/2+1; n++)  {
    r[n] = .5*log( amp_sq[n] ) ; // r[n] is now the natural log of the amplitude of FFT bin 'n'
} // end for                     // .5*log(A squared) = log(A)

cout<<"\n\nfrequency and amplitude estimation results\n\n n        magnitude      frequency       amplitude\n\n";
for (n = 370; n < 415; n++)  {
    fbin = n+.5+.5*c*(r[n+1]-r[n]);   // r[n+1] and r[n] are ln(b) and ln(a)
    F = (sample_rate/N)*fbin;
    amp = sqrt(c)*exp(r[n]+((fbin-n)*(fbin-n))/c)*(2.0*sqrt(twopi/2.)/N) ;
    printf("%2d\t%9.5f\t%9.5f\t%9.5f\n",n,amp_sq[n],F,amp);
} //end for

system ("PAUSE");
return 0;
} // end main
//******************** fft_recur ***********************
void fft_recur(long N, double *r, double *i)  {
long h, i1, j = 0, k, i2 = N/2, l1, l2 = 1;
double c = -1.0, s = 0.0, t1, t2, u1, u2;

for (h = 0; h < N-2; h++) {    // ***** bit reverse starts here ***
    if (h < j) {
       t1 = r[h]; r[h] = r[j]; r[j] = t1;
       t2 = i[h]; i[h] = i[j]; i[j] = t2;
    }
    k = i2;
    while (k <= j) {
          j = j - k;      k = k/2;
    }
    j = j + k;
}   //****** bit reverse done ******

for (k = 1; k < N; k = k*2) {
    l1 = l2;     l2 = l2*2;
    u1 = 1.0;    u2 = 0.0;
    for (j = 0; j < l1; j++) {
        for (h = j; h < N; h = h + l2) {
            i1 = h + l1;
            t2 = (r[i1] - i[i1])*u2 ;
            t1 = t2 + r[i1]*(u1 - u2) ;
            t2 = t2 + i[i1]*(u1 + u2) ;
            r[i1] = r[h] - t1;
            i[i1] = i[h] - t2;
            r[h]  = r[h] + t1;
            i[h]  = i[h] + t2;
        } // end for over h
        t1 = u1 * c - u2 * s;
        u2 = u1 * s + u2 * c;
        u1 = t1;  //x = u1 - u2;    y = u1 + u2;
    } // end for over j
    s = - sqrt((1.0 - c) / 2.0);
    c =   sqrt((1.0 + c) / 2.0);
} // end for over k

} // end fft

frequency and amplitude estimation results

 n        magnitude      frequency       amplitude

370       0.01156       15959.53142       0.00065
371       0.01328       16002.86947       0.00070
372       0.01542       16046.20441       0.00076
373       0.01810       16089.42371       0.00082
374       0.02138       16136.92108       0.00095
375       0.03020       16209.62311       0.00211
376       0.14116       16370.55303       3.20360
377      77.00512       16349.81449       0.95808
378     3193.91932      16351.21160       1.00396
379     24622.71009     16350.90700       0.99982
380     32939.04014     16351.10559       0.99850
381     7803.22610      16350.67285       1.01020
382     319.11953       16353.81879       0.87280
383       2.60322       16324.16039       7.57917
384       0.00113       16605.10228       0.00051
385       0.00724       16606.09917       0.00052
386       0.00851       16650.89175       0.00058
387       0.01072       16694.51462       0.00065
388       0.01381       16737.66957       0.00074
389       0.01787       16785.61197       0.00090
390       0.02814       16862.78585       0.00251
391       0.17568       17013.74935       2.83597
392      85.58946       16994.84074       0.95935
393     3413.62331      16996.22071       1.00408
394     25288.20406     16995.88948       0.99974
395     32472.56144     16995.79357       1.00043
396     7296.89833      16957.17339       4.09340
397      60.55471       17218.13259       1.24775
398     3330.74781      17212.33274       1.01810
399     25463.82887     17210.95691       0.99978
400     32350.54927     17210.69761       1.00169
401     7144.91966      17170.91626       4.35710
402      55.60789       17436.25628       1.40301
403     3423.65695      17427.15709       1.01474
404     25642.98140     17426.03425       1.00002
405     32244.80235     17425.93874       1.00073
406     7095.46971      17426.35995       0.98918
407     278.99401       17421.20777       1.27242
408       1.56515       17494.25360       0.02622
409       0.02946       17619.31213       0.00091
410       0.01521       17677.62446       0.00071
411       0.01452       17719.83432       0.00069
412       0.01340       17762.64014       0.00066
413       0.01224       17805.61294       0.00063
414       0.01113       17848.67454       0.00060

First, 4 tones of the input are generated (I use twopi when computing time because it's easier than writing it 4 times into the cosines).  Then, the Gaussian window is generated, after which it is applied to the input (by the way, the window shown here, which is from reference_1, is slightly wrong – I don't want to use the 'corrected' version of it because I don't want to write a page describing how I derived it).  Then we take a 1024 point FFT, followed by the (non-normalized) calculation of amplitude squared.  Then comes the calculation of the natural log of each FFT bin.  Finally, I compute the frequencies and amplitudes via the equations in reference 1 and reference 2above.

As you can see from the results following the code, we have an isolated peak between bin 379 and 380 from the first input frequency, which is estimated as 16350.90700 and an amplitude of .99982.  We can also see that the other 3 inputs exhibit a little 'peak merging' from bins 392 to 407, but the peaks and associated estimated frequencies and amplitudes are quite accurate.

The results could be even more accurate if I had: 1) adjusted parameters, 2) used the 'correct' window, 2) used a non-recursive FFT, or 3) used one of the variations of the technique (which I am not going to show here).

You might also notice that, as your value of the parameter 'c' gets smaller ('c' is related to the 'beamwidth' of the Gaussian window), you can actually improve some of your estimates, but it can cause some results to blow up due to numerical sensitivities.  The effective 'Q' or 'gain' of the Gaussian ratio technique is on the order of 5,000 – 10,000, and it can be much higher if you know how to adjust the parameters properly.

Of course, if your SNR is too low, you may have to increase N to get more integration time, subject to the limits of signal stability.


Positive and negative frequencies: how to prove they are independent?

This is a bit waffly, but gives the general idea.

The usual reason for this is because of expectations used when we talk about Power Spectral Densities:
$$
S_{xx}(\omega) = \int_{-\infty}^{+\infty} \gamma(\tau) e^{-i\omega \tau} d\tau
= \int_{-\infty}^{+\infty} E\left[ x(t) x(t+\tau) \right] e^{-i\omega \tau} d\tau
$$

So what you have to look at is how the $E[\cdot]$ works on the cross-products between the positive and negative frequencies. Assuming $x$ is just:
$$
x(t) = \sum_{k=1}^{N_+} a_k e^{j\omega_{+,k} t} +  \sum_{k=1}^{N_-} b_k e^{j\omega_{-,k}t}
$$
then each of the cross terms will look like:
$$
\lim_{T\rightarrow \infty} \frac{1}{T} \int_{-T}^{+T}e^{j\omega_+t}e^{-j\omega_-t}dt = \delta(\omega_+ - \omega_-)
$$
which will be zero if $\omega_+ \not= \omega_-$. That makes all the cross-terms zero, so you're just left with the $x_+(t)x_+(t+\tau)$ and $x_-(t)x_-(t+\tau)$ terms, i.e. your $A^2_{pos,k}$ and $A^2_{neg,k}$ terms.


In addition to @hotpaw2 explanation, a graphic. There are two analog square waves (red and green), with different lengths. They are depicted with a fine sampling, denoted by crosses. Their actual sampling is denoted by circles. The  red one is shorter than the green one, as can be seen in the interval $]0.7\;0.8[$. Yet, the  sample points are the same. Thus, from the circled samples, you cannot tell one from the other. This  would happen for any square wave whose end falls inside the  blue segment. You can see here the loss of information.

With the  most simple schemes  (using cardinal sines), these different sampled squares would be reconstructed as a single analog signal (as their  discrete samples are the same.) Depending on the  precision on the sampling, several analog-digital-analog signals can be theoretically recovered,  possessing a certain amount of ripples, or overshoot, or Gibbs phenomena.

You can reduce the  ripples in the interior of the square wave (as seen going from blue to red to green), but the amplitude of the overshoot will remain constant  at the  vicinity of the  edges.

As commented by @robert bristow-johnson, there are a number of alternative techniques. The document Alias-Free Digital Synthesis of Classic Analog Waveforms (Stilson and Smith) provides some of them. I add a recent technique, based on Finite-Rate-Innovation: Feasibility of {FRI}-based square-wave reconstruction with quantization error and integrator noise, B. He et al., ICASSP 2015.

Finally, for your last question "Why does superposition hold on such a simple construct?" This can be regarded as a problem of approximating a non-continuous function on a basis on continuous ones: you can get pointwise convergence (all samples converge), but no uniform convergence (not all converge at the same rate).


The PCM values are proportional to the sound pressure. 8-bit PCM is stored with an offset of 128, so PCM coded values from 0 to 255 code for 'pressure values' from -128 to 127.
These values are often normalized to the range -1.0..+1.0.
The PCM values are linear representation of the sound pressure, and a log operation is only required if you want to express the sound pressure in decibels.
To get actual sound pressure levels in Pascal or dB(A) you would need a calibrated .wav signal.


I have run an FFT now want to get the mean value in the 8hz - 13hz range?

Did you ever tried to follow example from MATLAB help? Result of the FFT is complex, containing information both for magnitude and phase. Thus in order to extract the magnitude you must take the absolute value of your result - abs function. What you are doing currently is plotting complex values.

Please refer to basics of DFT if you have any doubts (i.e. to understand why DFT spectrum is mirrored) and mentioned example for MATLAB fft. Here is one of ways to create frequency vector for magnitude response.

% fs --> sampling frequency to be defined.
N = length(data);
X = abs(fft(data))/N;
f = (0:N-1).*(fs/N);

If you wish to show your results in decibel scale, then convert the values with usage of 20*log10(X/ref); where ref is your reference value - can be max(abs(X)).

I am sure that having the frequency vector and magnitude values you can take it from here in order to calculate the mean value in appropriate range.

Extra things you might want to consider:

Apply the windowing to your signal if it is necessary (probably it is). When windowing is applied, scale not by N, but use the sum of window samples.
Use the Power Spectral Density estimation instead of FFT. Second figure you provided is the example of PSD - in fact this is copied from MATLAB help regarding PSD estimates using FFT. I believe you have all necessary info in there regarding periodogram.


What is the interpretation of Fourier Transform containing only imaginary part?

The integral doesn't exist in the conventional sense. How did you evaluate the limit $\lim_{t\to\infty}e^{-j\omega t}$?

If one allows distributions (generalized functions), such as the Dirac impulse (and its derivatives), as Fourier transforms then the class of functions that can be transformed becomes much wider, and the unit step is one of the functions that has a Fourier transform in that sense. However, its transform is not given by the expression in your question. The Fourier transform of $u(t)$ is

$$\mathcal{F}\{u(t)\}=\pi\delta(\omega)+\frac{1}{j\omega}\tag{1}$$

Take a look at this answer for a derivation of $(1)$.

A simple sanity check also shows that the expression in your question cannot be the transform of $u(t)$. A purely imaginary Fourier transform corresponds to an odd time domain function, which $u(t)$ is not. The expression in your question is actually a valid Fourier transform, but of a different function. Its inverse transform is $x(t)=-\frac12\textrm{sgn}(t)$, which is obviously an odd function.


How to deal with the mirrored components of a FFT? And another question

It depends on your spectrum analyzer.

If you want to analyze a one channel real baseband signal, the mirrored frequencies are redundant.  you should keep the first N/2+1 points.

For communications signals, there can be an I and Q (two reals as a single complex), you want to keep them all.  the frequencies will not be (in general) symmetric.

If you want to calculate the cross spectrum of two reals, there will be mirrored response, but if you want to do something like taking the inverse of the cross spectrum, you should keep both, but only display half.

It really depends on the actual application. There really isn't a right or wrong way. There are conventions and exceptions. There are niche applications. Spectrum analyzers that are commercial standalone test equipment, usually have lots of buttons, lots of modes.

for your second question, this is what spectral aliasing looks like if you don't apply  an antialiasing filter prior to the DFT.


Below summarizes efficient phase estimators for this application updated to include both a phase range of +/- 30 degrees and +/- 60 degrees. This is given in two parts, estimators for a real IF (intermediate frequency) signal, and estimators for a baseband complex signal. At the end are additional considerations related to acquisition.

Efficient Phase Estimators for Real IF Signals

Product detector: For real signals, a common phase estimator (detector) is a multiplier followed by a low pass. For this application where sensitivity is desired over $\pm 30$ degrees, the signals are nominally 90 degrees in phase resulting is an estimate that is proportional to the sine of the phase between two signals:

$$y(t,\phi) =A_1\cos(\omega_ct)A_2\sin(\omega_ct+\phi) = \frac{A_1 A_2\sin(\phi) + A_1A_2\sin(\omega_ct+\phi)}{2} $$

Where when followed with a low pass filter removes the time varying component resulting in:

$$y(\phi) = \text{LPF}[y(t,\phi)] =\frac{A_1A_2}{2}\sin(\phi) $$

Where $\text{LPF}[\cdot]$ is the time average provided by a low pass filter.

As a demodulator, this can be implemented either in a coherent receiver where $A_1\cos(\omega_ct)$ is estimated during acquisition and provided as an NCO in a digital implementation (or VCO in an analog), or in a non-coherent receiver (where the interest is in the phase difference between two successive symbols) where the demodulation is done by multiplying the received signal with a time delayed copy of itself, delayed by one symbol duration plus the duration of a quarter cycle of the carrier to convert $\cos$ to $\sin$ (when the IF carrier is sufficiently larger than the symbol rate):

$$y(\phi_2-\phi_1) = \text{LPF}\bigg[A\cos(\omega_ct+\phi_1)A\cos(\omega_c(t-T_s-T_c)+\phi_2)\bigg]$$

Where $T_s$ is the symbol duration in seconds and $T_c = 1/(4f_c)$ is a quarter cycle of the IF carrier in seconds with the IF carrier frequency as $f_c$ in Hz. Resulting in:

$$y(\phi_2-\phi_1) = \text{LPF}\bigg[A\cos(\omega_ct+\phi_1)A\sin(\omega_c(t-T_s)+\phi_2)\bigg]$$

$$ y(\phi_2-\phi_1) = \frac{A^2}{2}\sin(\phi_2-\phi_1)$$
$$ y(\Delta\phi) = \frac{A^2}{2}\sin(\Delta\phi)$$

A very efficient way to implementing either of the above approaches digitally is by hard limiting the input signal which reduces the above to a simple XOR of the most significant bit of the waveform. For accuracy, this requires ensuring the inputs to the XOR operation are a 50% duty cycle, but the result is linearly proportional to phase! It is usable over a $\pm 90°$ range with a linear phase result. Further, hard limiting a phase modulated waveform provides a 3 dB SNR improvement in positive SNR conditions (since all AM noise is removed), but can be more susceptible to jamming an interference (3 dB loss in negative SNR conditions). This is an approach to be considered due to the simplicity and high phase linearity.

As above, the X-or phase detector could be used in a coherent receiver where the NCO is also simplified to a 1 bit output (basically the MSB of a counter and you increment the count rate to adjust the frequency as part of a carrier tracking loop), or non-coherently where the MSB of the received signal is XOR'd with a delayed copy. As with the multiplier the two input signals would be in quadrature to center the detector over its unambiguous range.

Efficient Phase Estimators for Complex Baseband Signals

Given a generic signal as $Ae^{j\phi}= I + jQ$, the actual phase is given by $\phi =\tan^{-1}(Q/I)$ or $\phi =\sin^{-1}(Q/A)$.
The following is a summary of efficient phase demodulation approximations for a variation over $\pm30$ degrees and $\pm60$ degrees, assuming carrier recovery and timing is established during the $0 / 180$ acquisition period. Initial thoughts on approaches to efficient acquisition are also included at the bottom of this post.

Summary of Results

Below is a table summarizing the peak and rms phase error for various estimators. Estimators that were included in earlier versions of this post that offer no advantage to those listed below have been removed. As Ben suggests in the comments, the Q/A estimators are attractive for FPGA implementation since A is assumed to be constant over the duration of the packet.

Plots showing the relative performance are included below:

Detailed Descriptions

The estimators that are scaled by the envelope magnitude $A$ (Q/A, Q/A Juha and Q/Est(A)) are preferred since $A$ can be readily determined during acquisition of the 0/180 signal, and only needs to be determined once for relatively short packets, or is a parameter from the AGC otherwise. In a constant envelope phase modulated signal such as this, the received signal can be simply hard limited if there isn't a concern with potential 3dB loss with stronger out of band interference (or the complete loss from hard-limiting in the presence of a coherent jammer). Further, there is no need to actually divide by $A$, assuming $A$ is maintained to be constant over the packet duration the result will be linearly proportional to the phase and the decision thresholds can be set accordingly.

Q/A

$$\phi =\sin^{-1}\bigg(\frac{Q}{A}\bigg)$$
$$\frac{Q}{A} = sin(\phi)$$

for small $\phi$,  $sin(\phi) \approx \phi$ for $\phi$ in radians:
$$\phi \approx \frac{Q}{A}$$

Q/A Juha

Similar to  @JuhaP's suggestion in the comments of removing the linear slope error for the $Q/I$ estimator, here applied to the Q/A Estimator. The coefficient is found from the linear portion of the remaining terms in the Taylor Series expansion that weren't used, minimizing the error:

For ±30° Operation:
$$\phi \approx 1.0475\frac{Q}{A}$$

For ±60° Operation:
$$\phi \approx 1.150\frac{Q}{A}$$

Q/Est(A)

A fast and very efficient approach for estimating magnitude is the $\alpha$ max plus $\beta$ min algorithm where the maximum between $|I|$ and $|Q|$ scaled by coefficient $\alpha$ is added to the minimum scaled by coefficient $\beta$. At 30° range, $Q$ would always be the minimum and  $I$ always positive so this would simplify to $\alpha I + \beta|Q|$. A common choice for FPGA implementation is $\alpha = 1$ and $\beta =1/2$ since this minimized the error over all phases with bit shift divisions, but in this case $\alpha = 1$ and $\beta =1/4$ is a better choice given the narrowed phase range of $±30°$. If  multipliers were acceptable, the optimized coefficients are $\alpha = 0.961$ and $\beta =0.239$.  The plot below summarizes the two choices:

$$\phi \approx \frac{Q}{\alpha I + \beta |Q|}$$

option 1:  $\alpha =1$, $\beta = 0.25$
option 2:  $\alpha =0.961$, $\beta = 0.239$

Also not plotted below but shown above is the option optimized for use over ±60°:
 $\alpha =0.85$, $\beta = 0.45$

Note these are not optimized for the estimate of $A$, but to minimize the phase estimation error.

Q/I Phase Approximation

$$\phi =\tan^{-1}\bigg(\frac{Q}{I}\bigg)$$
$$\frac{Q}{I} = tan(\phi)$$

for small $\phi$,  $tan(\phi) \approx \phi$ for $\phi$ in radians:
$$\phi \approx \frac{Q}{I}$$

As @JuhaP mentioned in the comment, the linear slope component of the error could be removed by multiplying by 0.9289 resulting in (This one is labeld Q/A JuhaP in the plot). The coefficient below is slightly different than his suggestion but minimizes the error as it was found from the linear portion of the remaining terms in the Taylor Series expansion that weren't used rather than his approach of a first order terms of a polynomial fit to arctan:

$$\phi \approx 0.9289\frac{Q}{I}$$

Taylor Series Phase Approximations

The first term is the Q/A and Q/I approximations covered above for $\sin^{-1}$ and $\tan^{-1}$ respectively. Going beyond that is NOT recommended if efficiency is paramount but included for accuracy comparison.

arcsin

$$sin^{-1}(n) = \sum_{n=0}^\infty \frac{2n!}{2^{2n}(n!)^2}\frac{x^{2n+1}}{2n+1} \text{ for } |n|\le1$$
$$sin^{-1}\bigg(\frac{Q}{A}\bigg)= \frac{Q}{A} +\frac{1}{6} \bigg(\frac{Q}{A}\bigg)^3 +\frac{3}{40}\bigg(\frac{Q}{A}\bigg)^5 ... \text{ for } |Q/A|\le1$$

arctan

$$tan^{-1}(n) = \sum_{n=0}^\infty (-1)^n\frac{x^{2n+1}}{2n+1} \text{ for } |n|\le1$$

$$tan^{-1}\bigg(\frac{Q}{I}\bigg) = \frac{Q}{I} -\frac{1}{3} \bigg(\frac{Q}{I}\bigg)^3 +\frac{1}{5}\bigg(\frac{Q}{I}\bigg)^5 ... \text{ for } |Q/I|\le1$$

Using the first two terms for each results in:
$$\phi \approx \frac{Q}{A} +\frac{1}{6} \bigg(\frac{Q}{A}\bigg)^3$$.

$$\phi \approx \frac{Q}{I} -\frac{1}{3} \bigg(\frac{Q}{I}\bigg)^3$$.

A linear slope can also be removed from either of these with a gain constant multiplication as was done with the Q/A and Q/I estimators to further minimize the error.

Other Estimators

Juha Squared

@JuhaP offered this interesting estimator in the comments. Not very efficient but highly accurate with square terms:

$$\phi \approx \frac{3QI}{Q^2 + 3I^2}$$.

Acquisition

Efficient Acquisition for 0/180 preamble:

One idea that comes to mind for acquisition during the $0/180$ transitions is to use $\text{sign}(I_2)Q_1-\text{sign}(I_1)Q_2$ to get the change in phase between two symbols, which can the be corrected in a fast coverging and simple loop by derotating the incoming signal.  This appraoch would work well if the frequency offset is such that the phase does not rotate more than $\pm \pi/2$ between successive signals, otherwise a course FLL can be used first to get the offset within this acquisition range.
For a coherent receiver approach a PLL would be used to lock/track an NCO or PLL to the carrier and my squaring the received signal a reference tone at twice the carrier can be tracked for all modulations presented here (both the bi-phase acquisition interval and the 30 degree phase modulation when doubled will produce a distinct tone at 2x the carrier). Similarly a Costas Loop would track both signals while providing the reference signal that is nominally 90 degrees in phase with the carrier, thus providing both carrier recovery and phase demodulation.

Sources:
Taylor Series Expansions for arcsin and arctan:
https://proofwiki.org/wiki/Book:Murray_R._Spiegel/Mathematical_Handbook_of_Formulas_and_Tables


Detecting fundamental zero crossing point

I am not confident I have the best algorithm, but short of other suggestions I could offer you this idea:

Use an NCO (Numerically Controlled Oscillator) clocked to a known and accurate time reference to create a 50 Hz or 60 Hz reference signal, and lock the NCO to the power line signal using a Phase-Locked-Loop (PLL). The NCO will lock to the power line signal in quadrature phase with optimum noise suppression (based on the loop design and desired time averaging versus dynamic requirements). By reading the phase accumulator word the instantaneous phase offset between the power line signal and the time reference can be determined.

This idea is developed as follows:

If you have a known time reference and can generate a reference 50 Hz or 60 Hz sinusoidal signal, then you can use a phase detector (by multiplying the two signals and low pass filtering) to determine the phase difference between the two, and from the determined phase can determine the zero crossing relative to the known time reference. There will be calibration involved to offset the time delay of the low pass filtering. The advantage of this approach over single sample decision metrics (threshold zero crossing detectors) is its relatively high noise immunity based on the frequency cut-off (averaging time) of the low pass filter used.   The phase detector is based on the mathematical relationship between the product of two sinusoids at the same frequency; the result is an output magnitude that is proportional to the cosine of the phase of the two signals. This provides an unambiguous mapping of phase to magnitude over a 0° to 180° range.  Most significant is the result is averaged (via the low pass filter) over all samples taken, so benefits from significant noise averaging- EVERY sample in the waveform is used and contributes to the estimate of the waveform's phase position!

This is not yet complete as the primary drawback with this approach is that if the two signals approach 0° and 180° phase, the slope of the phase detector output approaches 0; meaning the result is not very sensitive (accurate) to phase shifts. The ideal operating point is when the two signals are in quadrature as the output would have the highest slope (sensitivity) to phase changes.
This leads to the improved approach of using a reference signal where you can precisely control the phase, and adjust the phase to the quadrature relationship between the signals by monitoring when the output of the phase detector is zero- this is basically how many phase locked loops are designed to operate, and you could use such a phase locked loop (PLL) approach to do this adjustment automatically; as long as you have accurate knowledge of the relationship between the control to your reference waveform and its phase. Further this gives you well designed control of your averaging time via the loop bandwidth of the PLL design. A numerically controlled oscillator is such a device where you have precise knowledge of its phase at any given time, based on its phase accumulator.

For more details on the NCO implementation, see this link: Numerically Controlled Oscillator (NCO) for phasor implementation?


Removing the 'red' artifact from my JPEG

Did you try to change the white balance of your image ? Apparently there's a function (balanceWhite) in OpenCV that tries to correct the white balance, you might want to give it a go.


The Guided Image Filter is a filter which assumes Piece Wise Linear model of the image.
Basically, it tries to estimate the best Linear Estimator per patch / pixel of the image.

Since it is Linear its coefficients are based on the first and second moments of the local patch.

You can see the code in my Q42415 GitHub Repository which also includes an implementation in Kernel Form in ApplyLocalLinearKernel().
Here are some results (Output of Q42415):


As Florent told: you must take fourier transform from your signal, by fft, then you will get your spectrum in normalized frequency, actual frequency is correspond to your sampling rate, then you can see amplitude and phase of each frequency(or harmonic) by these command: abs and angle.
For example, assume x is signal and T is time interval between your samples and

x=[10 -10 10 -10 10 -10 -10 10 10 -10 -10]
X=fft(x);%you can set frequency resolution by setting fft(x,resolution)
Fs=1000;
f=linspace(1,Fs,length(X));
plot(f,angle(X))
figure
plot(f,abs(X))


Impulse response rest conditions of LTI systems

Consider the second order LCCDE related to an LTI system with initial rest is given as:

$$a_0 y^{''}(t) + a_1 y^{'}(t) + a_2 y(t) = x(t) $$

The solution of the given equation for the input $x(t) = \delta(t)$ is defined as the impulse response of the LTI system and denoted as $h(t)$. The formal solution of the associated second order differential equation in time domain requires that a homogeneous and a particular solution be found and are added to yield the total solution; $$h(t) = y_h(t) + y_p(t)$$

Since the particular input is an impulse, then the particular solution is zero and therefore $$h(t) = y_h(t)$$ which is computed based on the initial conditions. For a second order equation you need two initial conditions and these are derived based on initial rest assumption so that:
$$y_h(0) = 0$$ and $$y_h^{'}(0) = 1/a_0$$.

Since most typically $a_0 = 1$ then $y_h^{'}(0) = 1$. This result comes from integrating the differential equation from $t=0^{-}$ to $t=0^{+}$ while recognising the fact that the output, $y_h(t)$, cannot involve any impulses in its defintion. In fact neither $y_h(t)$ nor any derivate of it upto (but not including) second order can include any impulses in their definition, because otherwise, then on the right side of the equation derivatives of impulse should have appeared, which is not the case, by definition, however. Hence

$$\int_{0^{-}}^{0^{+}} \left( a_0 y_h^{''}(t) + a_1 y_h^{'}(t) + a_2 y_h(t) \right)dt = \int_{0^{-}}^{0^{+}} \delta(t)$$

$$a_0 \int_{0^{-}}^{0^{+}} y_h^{''}(t) dt + a_1 \int_{0^{-}}^{0^{+}} y_h^{'}(t)dt + a_2 \int_{0^{-}}^{0^{+}} y_h(t)dt = 1
$$

Now because of the intuition based on the previous paragraph, the second and the third integrals turn out to be zero (since they don't include impulses, their integrals at a point from $0-$ to $0+$ iz zero) but the first integral includes an impulse and it will be nonzero as:

$$a_0 \int_{0^{-}}^{0^{+}} y^{''}(t) dt = a_0 ( y_h^{'}(0^{+}) - y_h^{'}(0^{-}) ) =  a_0 y_h^{'}(0^{+}) = 1
$$ and it follows that

$$y_h^{'}(0^{+}) =  \frac{1}{a_0} $$

since $y_h(t) = h(t)$ and since we can replace $0+$ with $0$, and taking $a_0 = 1$ as usual, then it follows that:
$$ \boxed{ h'(0) = \frac{1}{a_0} = 1 } $$

As the required initial condition on th impulse response $h(t)$.
Similar argumentation follows for systems of higher orders.


What is the difference between image restoration and image reconstruction?

The introduction of this paper explains the difference and gives an example. In short:

  Image restoration techniques presume that data are acquired in the image space; that is, the raw data represent a corrupted version of the image scene. In contrast, images are not directly observed in reconstruction problems. Instead, projections of an image are obtained from two or more angles, and the image scene is pieced together from these projections.


A voice signal is band-limited to 3.3 kHz. What is the Nyquist frequency?

If an analog signal is limited in frequency, and $W$ is its highest frequency component, then the minimum sampling frequency to avoid time-aliasing is $f_s = 2W$.

The Nyquist frequency is half the sampling rate of the system.

In your case, and to avoid time-aliasing, the Nyquist frequency would be $f_{Nyquist} = 3 300\text{ Hz}$. But, since an additional guard band is considered, then the Nyquist frequency is:

$$f_{Nyquist} = 3 300\text{ Hz} + 1 400\text{ Hz} = 4 700\text{ Hz}$$


Fourier Transform with both Time Delay and Frequency Shift

If you are ever unsure, just go back to the definition and work out the Fourier Transform property for the specific situation:

$$\begin{align*}\mathscr{F}\left\{x\left(t-t_0\right)e^{j2\pi f_0\left(t-t_0\right)}\right\} &= \int_{-\infty}^\infty x\left(t-t_0\right)e^{j2\pi f_0\left(t-t_0\right)} e^{-j2\pi f t}dt\\
\\
&= \int_{-\infty}^\infty x\left(\tau\right)e^{j2\pi f_0\tau} e^{-j2\pi f \left(\tau+t_0\right)}d\tau \\
\\
&= e^{-j2\pi ft_0}\int_{-\infty}^\infty x\left(\tau\right)e^{j2\pi f_0\tau} e^{-j2\pi f \tau}d\tau \\
\\
&= e^{-j2\pi ft_0}\int_{-\infty}^\infty x\left(\tau\right) e^{-j2\pi (f-f_0) \tau}d\tau \\
\\
&= e^{-j2\pi ft_0} X(f-f_0)\\
\\
\end{align*}$$

$$\begin{align*}\mathscr{F}\left\{x\left(t-t_0\right)e^{j2\pi f_0 t}\right\} &= \int_{-\infty}^\infty x\left(t-t_0\right)e^{j2\pi f_0t} e^{-j2\pi f t}dt\\
\\
&= \int_{-\infty}^\infty x\left(\tau\right)e^{j2\pi f_0(\tau+t_0)} e^{-j2\pi f \left(\tau+t_0\right)}d\tau \\
\\
&= e^{-j2\pi (f-f_0)t_0}\int_{-\infty}^\infty x\left(\tau\right)e^{j2\pi f_0\tau} e^{-j2\pi f \tau}d\tau \\
\\
&= e^{-j2\pi (f-f_0)t_0}\int_{-\infty}^\infty x\left(\tau\right) e^{-j2\pi (f-f_0) \tau}d\tau \\
\\
&= e^{-j2\pi (f-f_0)t_0} X(f-f_0)\\
\\
\end{align*}$$


How to correctly classify voice from other sound events?

It is very possible that you don't have enough data to train your classifier.  Although several thousand seems like a lot, it means that you only have a fraction of the actual thing you are trying to identify.  It may be helpful to expand your training data by slightly tweaking your samples to improve the robustness of your machine learning algorithm.  You could take the existing human voice samples you have and try slowing each one down slightly, increasing the speed slightly, introducing pink/white noise, etc...  This technique is commonly used with images in CNNs (by introducing small rotations and stretching different ways) to increase the size of the training set when training image classifiers.

Also, don't forget to reduce your learning rate as learning proceeds.  With too large of a step size your gradient decent could actually overshoot the loss function minima and result in poor results or saturate learning prematurely.

Hope that is helpful.


Phase difference between two sinuisoidal signal

You can find phase between 2 sines, if you write it as a complex signal where 1st sine is a real part and second one is imaginary. But to do so it's necessary to rescale one component to fit another in its amplitude variation. Or another approach is to normalize both sines in amplitude value. There are some explanations:

We have signals $s_1(t)$ and $s_2(t)$, which are assumed to be sine waves. Before the measurement of their phase differences, we have to normalize there signals. The normalization process is:

$s^{norm}_i(t) = \frac {s_i(t)} {1.414 \cdot \mathrm{E}[\mid s_i \mid]}$

where $\mathrm{E}[\mid s \mid ]$ is mean value of module of $s(t)$. To perform normalization in real time it's necessary to track $\mathrm{E}[\mid s \mid ]$ constantly since the signal $s(t)$ can be non-stationary. Recurrent moving average filter is the case:

$s^{mean}_i(k) = a \cdot \mid s_i(k) \mid + (1-a) \cdot s^{mean}_i(k-1)$

where $a$ is inverse of integration time (or equivalently filter's length $K$) $a = \frac {1} {K}$ and $k$ sample index in discrete time representation.

Once we normalize both signals it is possible to estimate phase shift between them as

$s(t) = s^{norm}_1(t) + j \cdot s^{norm}_2(t)$

$\phi(t) = arg(s(t))$

$\phi(t) = tan^{-1}(\frac {\Im{s(t)}}{\Re{s(t)}})$

If it's real time implementation good practice is to use CORDIC to calculate the phase instead of $atan$ function.


Doppler frequency when down-converting a receiving signal

When there is relative motion between the transmitter and receiver the frequency seen by the receiver is not $F_c$, so let's say you incur a Doppler shift of $\Delta_f$, then the shifted signal frqeuency is $F_c - \Delta _f$ now when you downconvert with $F_c$ the frequency offset of $\Delta_f$ still remains in the signal. That means the baseband spectrum "centred around freqeuency" of $0$ will be shifted by this value

In the receivers, this offset is estimated by the baseband signal processing (by performing correlation of time domain samples, linear shift in freqeuency domain means phase shift in time domain, so these correlations of samples simply leave the phase corresponding to the frequency offset due to doppler), there are other methods as well. Once these are estimated the clock at the receiver is further corrected with the estimation to tune the receiver to the "seen" signal freqeuency. So that the baseband spectrum is again centred around $0$ frequency


Do you want to remove the Yellow backgorund or make the Yellow light white?

If you want to color correct the image so that Yellow appears white then simply measure the R,G,B values of 'yellow', calculate the ratio between them. Then scale all the pixels in your image by that ratio.

Typically you would assume that green = brightness and just calculate red/green and blue/green then scale the red/blue pixels and just scale green to leave the overall brightness unchanged.

See http://en.wikipedia.org/wiki/Color_balance


How can I use a trained model to generate related images?

There are multiple aspects which render your approach infeasible:

SVM is a learning technique where high-dimensional data is mapped to lower dimension (a couple of hyperplanes, decision boundaries, which pass through support vectors). It doesn't have an inverse mapping: It is not super trivial to go from a point within a decision space back to the space of actual data. In other words, it's not generative. If your learning method supports going back from low-dimensional embedding to the higher dimensions space (such as GPLVM), your task could be easier.
HOG is not exactly reversible. Please check Hoggles, where the authors reveal what HOG features tend to understand, via formulating the visualization as in inverse problem. Therefore, it might not be able to synthesize consistent colors, or even detailed features. The synthesized image might end up looking just like a Hoggle image.

I would suggest the Generative Adverserial Networks for this task, which have been quite successful in image synthesis:

Here are some papers about it:

https://arxiv.org/pdf/1612.03242v1.pdf
https://arxiv.org/pdf/1406.2661v1.pdf

Here are some open source software:

https://github.com/junyanz/iGAN
https://github.com/carpedm20/DCGAN-tensorflow
https://github.com/jacobgil/keras-dcgan


Why there is no $dxdy$ and $\Delta k_x \Delta k_y$ in MATLAB's $\tt fft2()$ & $\tt ifft2()$ functions?

The Matlab functions fft,fft2,fftn and their inverses ifft,ifft2,ifftN implement what is called as the Discrete Fourier Transform (DFT) in $1$,$2$ and $N$ dimensions. The mathematical expression for those transforms is:

$$ X[k] = \sum_{n=0}^{N-1} x[n] e^{-j \frac{2\pi}{N} n k} ~~~ \text{, for} ~~~ k=0,1,2,...,N-1$$

$$ x[n] =  \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{j \frac{2\pi}{N} k n} ~~~ \text{, for} ~~~ n=0,1,2,...,N-1$$

Which is the 1D forward and inverse DFT in particular. Here x[n] is a discrete-time sequence of length $N$ samples and $X[k]$ is a complex valued discrete sequence of length $N$.

Note that the DFT can be considered as the samples of a continuous frequency function (the Discrete-Time Fourier Transform, DTFT) $X(e^{j\omega})$ which is given by:
$$ X(e^{j\omega}) = \sum_{n=-\infty}^{\infty} x[n] e^{-j \omega n } $$

for which the inverse is:
$$ x[n] = \frac{1}{2\pi} \int_{-\pi}^{\pi} X(e^{j\omega}) e^{j\omega n} d\omega $$

Note that DFT $X[k]$ is considered to be as the samples of DTFT $X(e^{j\omega })$ and both are operating on the same sequence $x[n]$ (assuming $x[n]$finite length) but produce two different outputs; one of them is a discrete sequence $X[k]$ while the other is a continuous function of $\omega$, $X(e^{j\omega})$.

Furthermore their inverses are the same; $x[n] = \mathcal{IDFT} \{X[k]\}$, $x[n] = \mathcal{IDTFT} \{X(e^{j\omega n})\}$. Hence this might have tempted you to think that the matlab function ifft2 was actually implementing the inverse DTFT by an approximation; which is not the case as it computes the inverse of the DFT which is a finite summation exactly computed in a computer.


Optimization of execution time of audio capture call in c/c++ (ALSA)

Capture audio into a FIFO, and process the samples in that FIFO in a separate thread.  That way, timing variations in your processing thread won’t block the audio input.

Very commonly the audio input FIFO is implemented using a lock-free ring buffer so it never blocks.  A separate processing thread can poll the buffer at some suitable fixed rate (video frame rate or network data rate for a given MTU size, etc.)


Derive Bode Diagram from difference equation

You can evaluate the frequency response by evaluating $H(z)$ in the unit circle: $H(e^{j\theta})$, where $\theta$ is the discrete-time frequency.

If you prefer a frequency in Hz, then use $\theta = 2\pi/f_s$ (i.e. the sampling frequency corresponds to $\theta = 2\pi$).

The above expression will give you the frequency response as a complex number, signifying both amplitude and phase response.

If you don't have access to complex number algebra on your platform, use either GNU Octave or Matlab to plot the frequency response (using the function freqz, which takes the a and b coefficients you have already).

Note that Bode diagram refers to an asymptotic approximation of the frequency response, valid only in continuous time systems.


Calculating values of frequency bins in Python

This is simply how Discrete Fourier Transform (i.e. Fourier Transform theory applied on sampled signal) works.
You get an output of length N if your input has length N, and after removal of symmetric part, what you get are $\frac{N}{2}$ points that span frequencies 0 (DC component) to Nyquist frequency ($\frac{F_s}{2}$).

For this same reason, the more points your signal have (let's say for a fixed sampling rate $F_s$) the more resolved your FFT output is since the spacing between each frequency bin in the output is:
$$
\frac{F_s/2}{N/2} = \frac{F_s}{N}
$$

UPDATE:

Regardless of the sampling frequency the FFT returns $N$ data points for an input with $N$ samples. What we say is that all the points of the output spread the circle: $[0;2\pi]$ in radian/sample or sometimes denoted $[0;1]$ in cycles/sample (called the normalized frequency). All these bins correspond to how many cycles, or oscillations you might have between each samples. You can measure at best half a cycle between each sample (for real signals). This is usually converted back into frequency simply by knowing the sampling rate.

So indeed if your signal is longer, you will have a longer output. Now what I was pointing out is that if your signal is longer but recorded with same $F_s$, i.e. it actually last longer in time, then you will resolve smaller frequency (since you have more time in total, this gives finer information in the frequency domain). That is, the bins are closer to each other in frequency domain ($\Delta f = F_s/N = 1/T$ get smaller). So basically, the longer you record, for whatever sampling rate, the finer resolution you get in the frequency domain.

In summary, if you do not know $F_s$ the bins are located (after removal of symmetrical part):

span $[0,0.5]$, and spaced by $1/N$, and the unit is cycle/sample, OR
span $[0,\pi]$, spaced by $2\pi/N$, unit is radian/sample

Then if you know $F_s$, simply transform back into Hertz, knowing that it must span $[0, F_s/2]$.


Computing narrowband spectrograms using MATLAB

No one is going to give you a full working code - the point is that you should learn and understand how to do it. Following guidelines for you:

As given, use the wavread function to load the recording (please mind that in the newest 2014a version this function is deprecated, you should use audioread instead). MATLAB has wonderful help with very good examples, so simply use it.
Having your signal loaded, together with its sampling frequency, you need to split it into chunks (many short-time frames). This can be done in two ways, either you are using a loop construct (for, while - whichever you prefer), or (my favourite) use the MATLAB buffer function (please consider using 'nodelay' option in it, as first frame will have zeros in it up to overlap length). Although you require the Signal Processing Toolbox license (I'm sure your University has it). This function allows you to specify the length of frame, overlapping and other things. In result you will obtain the matrix which columns are your signal frames. Having the sampling frequency from wavread, it is very easy to calculate corresponding length of your frame in samples for given number of miliseconds
Now you can iterate over the columns of your matrix and calculate the FFT for each of them. You might want to use the fft function to get it. In extra parameter you can specify the length of your signal (to apply zero padding eventually). So let's say that now you have FFT for first frame. Now take the magnitude of it (abs), and choose first 512 samples. Store this result into fist column of your super matrix $S$. Optional step: before calculating the FFT, I suggest you to multiply your frame by the window function, i.e. hamming window.
Now you have the matrix $S$, so plotting can be done by using various functions (pcolor, imagesc, and so on - pick whichever you prefer).

I suggest you to compare the final result with MATLAB built-in spectrogram function.

Good luck!


At the start of the inner loop, you assign

uchar* t_p= (uchar*) temp->imageData;

later on, you use it as

t_p[x+cx];

Aren't you missing the cy here?


I don't have experience with continuous time kalman filters.  However from your description it sounds like you are making measurements of position and acceleration over time.  If these measurements are sampled, meaning you have individual measurements associated with time, you should be using the discrete time kalman filter.  I feel comfortable broadly stating if you are trying to implement the filter on a pc (matlab), microcontroller, fpga, or dsp based on a set of measurements you are implementing the discrete time filter.

EDIT

Based on your edits, I more clearly understand your question.  I don't have experience with Matlab's built-in kalman filter functions but a quick read of the comments in kalmd seem to indicate to me you want to use kalman and not kalmd.  Below is a snippet pasted from kalmd that should make it clear.

%   The LTI system SYS specifies the plant data (A,[B G],C,[D 0]).
%   The continuous plant and covariance matrices (Q,R) are first
%   discretized using the sample time Ts and zero-order hold
%   approximation, and the discrete Kalman estimator for the
%   resulting discrete plant is then calculated with KALMAN.

This text tells me the function is expecting a continuous function which it then evaluates/discretizes, and then calls the kalman function to create a discrete filter.  From the sounds of it you already have discrete data (measurements over time).  Hence kalman is what you want to use.


Does a closed loop system with a reference of its previous iteration always reach a point in which it adventually loops, even if that loop is stable?

As is usual with these sort of questions, the answer is "it depends".

For example, suppose
$$
F(x) = x
$$
then the loop is immediate.

For example,  suppose
$$
F(x) = x/2
$$
then the loop never occurs unless $x = 0$.

The same will be true for more complex $F$s.

For example
$$
F(x) = \mathbf{A} x
$$
where $\mathbf{A}$ is a matrix. If $\mathbf{A}$ is a rotation matrix of a certain type (so that it rotates by $360/n$ degrees where $n$ is an integer) then the period will be $n$. If $\mathbf{A}$ is something else, then all bets are off.


Obtaining Correct (Least Squares Sense) Affine Transform Parameters Between Two Images

As a complement to the comments above, I am adding the slides for Image Registration by Leow Wee Kheng, where one can find:

  Given two images, how to register one with the other?

  Determine the corresponding points between the images. Manually mark corresponding points, or
  Detect and match features between views (see lecture on feature detection and matching).
  Determine the transformation between corresponding points. Assume that all pairs of corresponding points are related by the same
  transformation. Compute parameters of transformation given
  corresponding points.

Importance is laid on corresponding points. Indeed, even with an over-determined system, if the chosen pixels an weakly informative (on flat or noisy parts of the image), the matching might not be efficient. One can instead detect features (edges, corners, textural patches), and compute the registration based on them.

Sometimes, one does not have the same points, or the same number of features between images. I recently have tested a efficient algorithm (CDP, or Coherent Point Drift) for non-standard chromatographic images (made of peaks or blobs), that can be tuned from rigid to smooth deformations:

Point-Set Registration: Coherent Point Drift, 2010
Our application, for the record: BARCHAN: Blob Alignment for Robust CHromatographic ANalysis, 2017


Random 2D noise generated in spectrum domain always have maximum at (1,1)

If you use random real valued frequencies (e.g. just magnitudes), they all correspond to cosine functions, which all peak at zero (cos(0 * f) == 0).  To get more random peak locations, use random complex frequencies as the input to your ifft2d, which randomizes the phase so all the peaks don't line up at zero.


Information to be transmitted being a stochastic process?

The excerpt is intimating that, for communication applications, the information source is stochastic in the sense that its output is not fully known a priori. If you did know its output, then you wouldn't need a communication link at all! With that said, there are levels of gray to this effect, and there is a large area of study on its own known as information theory. Here's a very basic high-level summary:

The amount of information conveyed by a transmission is related to the amount of entropy contained in the transmission. Its entropy measures how unpredictable the transmission is. In your example, if a receiver observes a signal that is fully deterministic (i.e. it is completely predictable), then that transmission contains zero information. So, in a sense, the less we know what to expect about the next observed transmission, the more information that it can carry.

So, in terms of maximizing information rate, you want to maximize the entropy in your transmitted signal. However, there are some cases where reducing your transmission's entropy could be advantageous. Specifically, error-correcting codes are encoding methods that allow a receiver to recover from symbol errors by exploiting intentionally-inserted redundancy into a transmitted signal (which reduces its entropy). Thus, the information rate is lower, but the probability of error is reduced as well.


Traditional (analysis) filter banks are banks of linear filters followed by downsampling. This  is not so efficient, because you convolve and then subsample: possibly a waste of operations.  People have found ways to invert the operations, downsample and convolve, with other filters involve the polyphase components (even and odd ones in the two-channel case) of the signal (polyphase matrix). Such matrix can be factored into a product simpler matrices, like diagonal, upper and lower diagonal ones, acting more separately on the polyphase components.
A reinterpretation consists in using the upper matrix on some polyphase component (say, odds) to predict and combine outcome to other components (say, even) to get what is called an update, possibly with scaling factors.
Properties of polyphase matrices tell that this is always possible, yet not unique. And one can save some operations, compute inverses more easily (upper and lower matrices have simpler inverses). So all invertible linear filter banks can be converted in lifting.
But more importantly, one can insert non-linear blocks and still be invertible and efficient. So the lifting is more generic that "linear" filter banks (as long as they have finite support). An example of the above is known as the S-transform: given a sequence of integers  $c[n]$, you can decompose it into ($\lfloor \cdot\rfloor$ denotes rounding down):

$l[n] =\Large\lfloor \frac{c[2n]+c[2n+1]}{2}\Large\rfloor$
$h[n] = c[2n]-c[2n+1]$

How to design a lifting scheme from scratch requires to:

understand how to split the data into several subsets
design proper predict and update operators
choose how to iterate them on different levels.

I do like a couple of papers from people at KU Leuven:

Integer wavelet transfroms using the lifting scheme, Geert Uytterhoeven, Dirk Roose, Adhemar Bultheel
Wavelet Transforms Using the Lifting Scheme, Uytterhoeven, Geert  and Roose, Dirk  and Bultheel, Adhemar
WAILI: Wavelets with integer lifting, G. Uytterhoeven, F. Van Wulpen, M. Jansen, D. Roose and A. Bultheel


Have a look at these two python libraries that provide a number of audio features easily from WAV files, including MFCC.

Librosa: MFCC docs, github
Madmom: MFCC docs, github

Good luck!


How to 'interpret' the Fourier Transform (specifically, of a convolution kernel)

Convolving a signal with a kernel is equivalent to applying a filter to it. The kernel is the impulse response of the filter, and the Fourier transform of the kernel is thus the frequency response of the filter.

In your case, the filter's impulse response is a rectangular function of width 2 and centered at 0. You can interpret this filter this way: for any $t$, the filter output $y(t)$ will be the average of the input signal $x(t)$ over the interval $[t - 1, t + 1]$. This is an averaging, low-pass filter, and indeed, the Fourier transform of the kernel is a decreasing function, showing that the higher frequencies are attenuated.


The difference is pretty quickly explained: the CTFT is for continuous-time signals, i.e., for functions $x(t)$ with a continuous variable $t\in\mathbb{R}$, whereas the DTFT is for discrete-time signals, i.e., for sequences $x[n]$ with $n\in\mathbb{Z}$.

That's why the CTFT is defined by an integral and the DTFT is defined by a sum:

$$X(j\Omega)=\int_{-\infty}^{\infty}x(t)e^{-j\Omega t}dt\tag{1}$$

$$X(e^{j\omega})=\sum_{n=-\infty}^{\infty}x[n]e^{-j\omega n}\tag{2}$$

Note that in both cases, the frequency variables $\Omega$ and $\omega$ are continuous variables.

There is also the discrete Fourier transform (DFT), which can be applied to finite length (or periodic) sequences:

$$X[k]=\sum_{n=0}^{N-1}x[n]e^{-j2\pi kn/N}\tag{3}$$

It is easily shown that DTFT and DFT are related as:

$$ X\left(e^{j\omega}\right)\Bigg|_{\omega = 2 \pi\frac{k}{N}} = X[k] $$

when $x[n]=0$ for all $n<0$ or $n\ge N$.

Note that the DFT of the length $N$ sequence $x[n]$ is also a length $N$ sequence $X[k]$. The DFT is important because there exist efficient algorithms to compute $(3)$. This class of algorithms is called Fast Fourier Transform (FFT).

The DFT can be used to approximately compute samples of the DTFT and of the CTFT. For finite length sequences, the DFT computes (exact) samples of the DTFT. For longer (or infinitely long) sequences, there is always a truncation error caused by considering only a finite length window of the data. When approximating the CTFT, there are two error sources: the truncation error, and the aliasing error caused by sampling the continuous-time signal. Since a signal cannot be time-limited and band-limited at the same time, at least one of those two errors is non-zero.


This is just about obtaining a symmetric signal after zero-padding. Take a symmetric signal (w.r.t. to time index $n=0$) and append zeros. Due to the implicit periodicity of the time signal used as input to the DFT, it does not matter if you append all zeros to the right, to the left, or on both ends of the symmetric signal. In any case, since the DFT/FFT assumes the first sample in the buffer to be the time origin, the input vector to the FFT needs to be arranged as

[right half of the signal | left half of the signal]

which is basically one period of the periodically continued zero-padded signal between $n=0$ and $n=N-1$, where $N$ is the FFT length, i.e. the length of the signal after zero-padding. Note that now all the zeros end up in the middle of the vector to be transformed. If done correctly, the result of the FFT will be purely real-valued (apart from numerical errors). If the signal is asymmetrical, its FFT will be purely imaginary.


Overlap and windowing a signal for signal processing

Don't overlap.  Don't use a Von Hann window.  Zero-pad by at least the length of the impulse response of your desired equalization filter.  Then use an overlap add or overlap save fast convolution method or algorithm to combine your zero-padded and filtered FFT/IFFT blocks.

Otherwise circular convolution artifacts will corrupt (severely distort) your attempts at filtering your signal.

If you want to analyze or visualize the audio, use separate overlapped and windowed FFT blocks, not the ones used for overlap add/save fast convolution filtering.


What is the relationship of S-parameters to frequency response?

Let's use an example of a Wilkinson power divider that is used to take an input and split the power in half between two branches.
The ideal S-parameter matrix then looks like this
$$S = \begin{bmatrix}
s_{11} & s_{12} & s_{13} \\
s_{21} & s_{22} & s_{23} \\
s_{31} & s_{32} & s_{33} \\
\end{bmatrix} =
\begin{bmatrix}
0 & \frac{-j}{\sqrt{2}} & \frac{-j}{\sqrt{2}} \\
\frac{-j}{\sqrt{2}} & 0 & 0 \\
\frac{-j}{\sqrt{2}} & 0 & 0 \\
\end{bmatrix}
$$
Practically speaking, you are interested if the device is doing its job. In this case, you want to see if the power divider is indeed splitting the input power in half for the frequencies of interest. In a real divider, there will be deviations on how it performs over frequency with the hope that the best performance is achieved for the bandwidth of interest. To see this, you can plot the S-parameters which are a type of frequency response.
Take this example of a Wilkinson power divider designed to operate at 2.5 GHz

"Frequency response" is a general term for a particular output you want to see over some frequency range. As you suggest, you can plot the S-parameter values over their frequencies, which is a type of "frequency response".


$\omega$ is angular frequency in radians, $\omega =2\pi fT$, where $f$ is frequency (for example in Hz) and $T$ is sampling period (for example in seconds). This is revealed by Eq. 3:

  $$Y(\omega, k) = \textstyle\sum_{n=-\infty}^\infty y(n)w(k-n)e^{-j\omega n}\tag{3},$$

which represents discrete-time Fourier transform (DTFT) of windowed signal $y(n)w(k-n)$.

It's not common to call $\omega$ "discrete angular frequency index", which gives just 3 google hits.


Determining LTI system response to a scale change

I think it's instructive to look at this in the time domain. For LTI systems, the output $y(t)$ is given by the convolution of the input $x(t)$ with the impulse response $h(t)$:

$$y(t)=(x\star h)(t)=\int_{-\infty}^{\infty}x(t-\tau)h(\tau)d\tau\tag{1}$$

From $(1)$, the response to $x_a(t)=x(t-a)$, $a\in\mathbb{R}$, is given by $y(t-a)$:

$$\begin{align}(x_a\star h)(t)&=\int_{-\infty}^{\infty}x_a(t-\tau)h(\tau)d\tau\\&=\int_{-\infty}^{\infty}x(t-a-\tau)h(\tau)d\tau\\&=y(t-a)\end{align}\tag{2}$$

which proves time-invariance.

However, the response to $x(at)$, $x\in\mathbb{R}$, is

$$\int_{-\infty}^{\infty}x(a(t-\tau))h(\tau)d\tau\neq y(at)\tag{3}$$

Consequently, in general we cannot deduce from the response $y(t)$ to an input $x(t)$ the response to a time-scaled input $x(at)$.


How to scale Phase Transition Diagram for Compressed Sensing?

I recommend first fixing $N$. By the way, remember here that compressed sensing works asymptotically well; for $N \rightarrow \infty$ your phase transition is going to move up/left (better) in the $\delta, \rho$ diagram, see for example Donoho & Tannner 2010.

Your size of $N$ essentially determines the resolution you can achieve in your phase transition diagram. For example we can look at Maleki & Donoho 2010, Section IV: they describe how they select $N = 800$ and then calculate $n$ and $k$ from a grid of values $(\delta, \rho)$ between 0.05 and 1 in 30 steps. For each value $\delta$ or $\rho$, the corresponding $n$ or $k$ are calculated as $$n = \lceil \delta N\rceil,$$ respectively $$k = \lceil \rho n \rceil.$$

As you can see there is generally some imprecision involved during to the rounding. If you select a too fine grid of values $(\delta, \rho)$ for a given $N$, several of the resulting $n$ or $k$ values are going to round to the same number. You will have to either select a larger $N$ or choose a coarser grid for your $(\delta, \rho)$ to avoid this.

Finally, you should (if you want the same type of plot as Donoho et al.) plot your success ratio in each of your points in a $(\delta, \rho)$ coordinate system. Your figure above looks as if it was perhaps plotted with $k$ on the first axis?

Now, if you wish to estimate the actual phase transition boundary, you can also see in Maleki & Donoho 2010, Section IV, how they fit a logistic curve to the success ratios and can then read a specific percentage reconstruction boundary $\rho$ for each value $\delta$ from the fitted logistic function.


Connection between SNR and the dynamic range of the human ear

Where does this quote come from? It does NOT look like a reliable source of information to me.

The dynamic range of the 16-bit quantizer is 6B = 6·16 = 96 dB.

Wrong. The quantization noise for a 16 bit bipolar rounding quantizer is about -101dB dFS. That makes the signal to noise ratio for a sine wave 98dB and for the "average" music file about 8dB-90dB.

Note that the dynamic range of the human ear is about 100 dB.

Grossly oversimplified and misleading. The dynamic range of human auditory systems depends highly on frequency (and unfortunately also and age and how much rock& roll you have played, sigh). Max dynamic range is about 120 dB at the midrange.

This is the reason why “CD quality” digital audio requires at least 16-bit quantization.

16-bit quantization was chosen as a complicated trade off between audibility of quantization noise and cost & technical complexity at the time the standard was created. 16-bit per se is in many cases not sufficient unless proper "dithering" and "noise shaping" is applied as well.
The explanation really misses the main point. SNR and audibility of quantization noise depends A LOT on frequency and the specific signal. 16-bit PCM quantization is a "brute force" approach with most bits being wasted and a few important ones missing. Perceptual codecs such as AAC, MP3, Ogg, etc. do this much better and generate way more "sound quality per bit".
Finally, the question:

How can we obtain(calculate) the fact that the dynamic range of the human ear is about 100 dB?

Through listening test performed over a wide range of diverse subjects and test signals. The hearing threshold is determined as the sound pressure where the subject says "I can't hear this anymore" and the pain threshold corresponds to "ow, that hurts. Stop it!". The pain threshold is particularly difficult to test since it's fairly subjective and the testing can cause physical damage to the test subject. The results of the hearing test are then compiled using statistical analysis. You end up with a set of curves like this: https://en.wikipedia.org/wiki/Equal-loudness_contour .


How to determine SNR from FFT of measured data?

Thanks for the suggestions. I ended up just picking out the biggest nearby peaks, and if the difference between the frequency of interest and the biggest peak is less than 10dB, save a plot for manual inspection. I found it to be too difficult to distinguish between all the different scenarios that might be "good" or "bad" data.


How to analyze audio input/spectrum correctly

Let's see. I have no idea about Java classes or what type of support it has for signal processing, but i will give you some guidelines. The Particular details of the implementation in the language, you have to figure out.

To do any kind of processing on an audio file you need the "raw" data, meaning an audio file wich has uncompressed audio samples. For exmaple, the WAV format, which normally has samples in 16bit signed 2's complement.

If you are going to use a computer and do offline processing, might be a good idea to cast that value to a double (Which in most processors today means 64bit floating point number).

Therefore you will have a continuos stream of double numbers coming from the file. Then you need to define a way to handle the continous data. A standard, very widely used way is to use a circular buffer (even if it is offline, i assume you want your code to be efficient, and loading the entire file to memory as an array is not the best solution). Either that or just a normal buffer, is your choice. The buffer length should be a power of 2 (Do to efficiency of the Cooley-Tukey radix-2 algorithm).

Now you need to do the actual FFT. This is merely a multiplicacion of the buffer (Which mathematically is a vector) by the FFT matrix. How this operation is actually performed in Java, no idea. In C it would be just to pass the pointer to the array and the length to an FFT routine that, either returns a pointer to dinamically allocated memory, or leaves the result in an array that you pass to it.

finally, you arrive at an array of M complex numbers (Assuming the length of the array/buffer/vector with the signal is M). And then you do whatever you want with it.

For instance, you could take the magnitude of each of the complex numbers, and find the maximum, to detect where the fundamental frequency could be (very approximate though).

Extras:
More advanced techniques to process would include a pre windowing to avoid leakeage, doing zero-padding to obtain more resolution on the windowed spectrum, etc.

Hope it Helps.


How can I change parameters into a text string in MATLAB?

You can use this code

text = 'cmor1-';
for i=1:10
     s{i}=[text num2str(i)];
end

where s{i} contains the string cmor1-i


Concepts: Ideal and Non Ideal Channel Resulting in Inter Symbol Interference

The book doesn't say that the impulse response must be zero for an ideal channel. It says that an ideal channel has exactly one, and not more than one, non-zero component, i.e. the ideal channel's impulse response is an impulse, which means that the signal is only delayed but not distorted.


How Could One Imitate Photoshop's Vibrance Filter with OpenCV?

The Vibrance in Photoshop is a mask guided saturation.
What does it mean?
It means we add saturation only to pixels marked by the mask while leaving others untouched.
The mask evaluates the original level of the saturation of each pixel.
If the pixel is highly saturated, then when Vibrance is applied this pixel will be excluded.
Namely vibrance works only on less saturated pixels.

How could we create such a mask?
We transform each pixel from RGB Color Space to HSL (Hue, Saturation, Lightness) Color Space.
Then we take only the Saturation values for all pixels and normalize into the [0, 1] range where 0 is not saturated pixel and 1 is fully saturated pixel.

Then we add saturation only to pixels with low value of saturation.

I'm not an expert on OpenCV, but the above should be pretty easy to integrate in OpenCV (Given it has RGB to HSL Conversion, You can use HSV or HLS)


The problem is a bit more complicated that mere denoising, since a highly non-linear operator affects the original signal. Restoration or recovery would be more appropriate.

The floor function acts as an instance of data quantization (and possibly saturation or clipping). And adding noise is also known as dithering, a method to  reduce sensitive quantization effects (dequantization).
Those are keywords you can use along with: restoration, recovery, inference... Without more information on the signal's dynamics, models, and the noise, a couple of pointers that I am aware of (to be updated):

Wavelet denoising of coarsely quantized signals
quantized compressing sensing
piecewise regular approximations


How to convert filter into frequency domain to do filtering?

Your filter is an FIR filter, therefore its coefficients are simply the impulse response $h[n]$. Your signal $x[n]$ can be filtered in time domain by convolving with impulse response: $$y[n]=x[n]\star h[n]$$
Knowing that convolution in time domain is equal to multiplication in frequency domain $$Y[f]=X[f]H[f]$$
you only need to calculate two DFT's.

One thing to remember - output will be of the same length as input signal.

Speaking Matlab'ish:

filter_2 = firceqrip(2,0.6,[0.05 0.03]);

freqs = [0.08, 0.2, 0.32, 0.4];
periods = 1./ freqs;
t_max = 4 * periods(1);
t = linspace(0, t_max, 50);

signal = sin(2*pi*0.08*t) + sin(2*pi*0.2*t) + sin(2*pi*0.32*t) + sin(2*pi*0.4*t);

%% NEW CODE HERE %%
% Output is same length as input signal
M = length(signal);
% Find next power of two for sake of calculations
M2 = 2^nextpow2(M);
% Calculate the FFT's
Hf = fft(filter_2, M2);
Xf = fft(signal, M2);
% Multiply in frequency domain
Yf = Hf .* Xf;
% Perform the IFFT
y = ifft(Yf);
% Take only M samples
y = y(1:M);
% Remove this nasty imaginary parts
% comming from round-off errors
y_r = real(y);

% Call the in-built function for sake of comparison
y_l = filter(filter_2, 1, signal);

% Plot the original signal
subplot(311)
plot(signal)
grid on
title('Original signal')

% Plot the results
subplot(312)
plot(y_l);
hold on
plot(y_r, 'r--')
grid on
legend({'in-build function', 'user function'})
title('Filtered signal')

% Plot the relative error
subplot(313)
plot(abs(y_l - y_r))
grid on
title('Relative error')

Comparison of results below. As you can see, it matches perfectly MATLAB filter function.


Decimator effect on wide sense stationary input

The decimator by integer $M$ can be shown to be the following block:

$$ x[n] \longrightarrow \boxed{ \downarrow M } \longrightarrow y[n] = x[Mn] $$

Assuming that the input $x[n]$ is WSS it has the ACS as
$$ r_x[k] = E\{ x[n] x^*[n+k] \} $$

Then the output autocorrelation can be defined as:
$$ r_y[n,n+k] = E\{ y[n] y^*[n+k] \} = E\{ x[Mn] x^*[Mn+Mk] \} = r_x[Mk]$$

as can be seen, the output auto-correlation $r_y[n,n+k] = r_y[k] = r_x[Mk]$ also depends on the lag $k$ and therefore we can conclude that $y[n]$ is also WSS.

(of course the mean of $y[n]$ is also independent of time $n$.)


Difference between real and complex impulse responses?

Your question might be the result of a misunderstanding. There is only one type of impulse response, and this is a system's response to a Dirac impulse $\delta(t)$ (in continuous time), or - in discrete time - to a discrete impulse $\delta(n)$. For a real-valued system, i.e. a system with real-valued parameters, the impulse response is always real-valued. For a complex-valued system, i.e. a system with complex-valued parameters (e.g. complex-valued filter coefficients), the impulse response is complex-valued. Note, however, that complex-valued systems are just a convenient way to describe certain real systems, just like complex numbers are often more convenient to establish certain mathematical results.


Question #1: Are these 2 methods equivalent?

No, method #1 is taking an FFT of the entire sequence resulting in frequency bin widths of fs/N = 1024/4096 = 0.25 Hz. Method #2 is summing shorter FFTs at a frequency resolution (bin width) of 1024/512 = 2.0 Hz. For #2, if you divide by the number of FFTs after you sum them you are then averaging them. This is often times used to reduce the variance of the result.

  Question #2: Is the frequency resolution for both methods equal to 0.25Hz per bin?

No, they are 0.25 & 2.0 Hz. The resolution, or bin width, of the FFT is fs/N.


This was my first google hit: https://medium.com/analytics-vidhya/understanding-the-mel-spectrogram-fca2afa2ce53
It seems to explain both how and why?
Mel frequency binning is afaik only used when human perception is involved. Either directly, or for signals that have probably been adapted to human hearing (such as human speech). It is a way to represent signal information in a way that mimics «low level auditory perception», or could be seen as removing perceptually redundant information.
I have not understood why it is that we generally have the same temporal reolution in narrow low frequency bands, as in wide high frequency bands. Historically, it may have been due to limited processing power and the efficiency of a single STFT windowing process. And/or it could reflect the limited neuron firing rate after the acoustic-mechanical filtering carried out by our ear.
But from a dsp perspective, I would expect a variable time/frequency resolution, where high & wide bands gets more temporal resolution.


I'm not comfortable with either of your two equations, in part because you don't define any of the three variables involved, but also because I think there's an easier way to understand $E_b/N_0$ and simulate it.

I think the most important thing to know about $E_b/N_0$ is that it is measured at the matched filter's output. The first consequence of this is that the measurement is done at baseband.

Let's assume that your signal $s(t)$ is linearly modulated, so $$s(t) = \sum_k a_k p(t-kT_p),$$ where $R_p = 1/T_p$ is the baud rate and the information is transmitted by the symbols $a_k$ which are elements of an $M$-PAM constellation with average energy $E_p$. Furthermore, the Nyquist pulse $p(t)$ has unit energy.

Let's assume the system has zero noise. In this case, the output of the matched filter, sampled at rate $R_p$, is equal to $a_k$ and its average energy is $E_p$. Then, the average energy per bit at the matched filter's output is $$E_b = \frac{E_p}{\log_2(M)}.$$

Now let us consider the noise. Assume that the input to the matched filter is $n(t)$, which has power spectral density $N_0/2$. Then the filter's output sampled at rate $R_p$ is a Gaussian random variable with variance $\sigma_n^2 = N_0/2$.

The signal-to-noise ratio is defined as the ratio of the average energy per bit to the average energy of a noise sample: $$\text{SNR} = \frac{E_b}{\sigma_n^2} = \frac{2E_b}{N_0}.$$ You often see the SNR defined instead as $E_b/N_0$; the reason is that, in a quadrature system, the total noise energy in the receiver is $2\sigma_n^2 = N_0$. As long as the definition is made clear, one could use either one.

Let's see how to simulate this. Say you want to simulate 4-PAM at $E_b/N_0 = 10$. First, generate 4-PAM symbols with $E_b=1$, or $E_p = 2$. Then, generate noise with variance $\sigma^2 = N_0/2 = E_b/(2 \times 10) = 1/20$. Assuming Matlab:

% generate 100 random a_k
a = 2*(randi([0,3],1,100)-1.5);
% scale them so that Ep=2, Eb=1
a = a*sqrt(2/5);
% Noise has N0 = 1/10, so N0/2 = sigma2 = 1/20
n = sqrt(1/20)*randn(1,100);
% matched filter's output:
y = a + n;

Now y has $E_b/N_0 = 10$.

Note that extending this to quadrature modulation is trivial: just think of the quadrature system as two parallel PAM systems, each with noise variance $N_0/2$.

I realize that I have only indirectly answered your question, but I hope this is still useful.


Predict auditory masking with a crude estimation?

By looking at the spectrum's probability law you are completly missing the temporal information in the two signal you wish to compare thus it would be a bad predictor.
Imagine you analyse two 2048 pts signal, one is null during the last 1024 points whereas the other one is null during the first 1024 pts, they could have the exact same spectrum but they won't interfere each other at all.
This is the whole purpose of time-frequency analysis, unless you have some a priori knowledge, you can't avoid using a short time fft.


can $\frac{1}{H(z)}$ be causal and stable?

For 1/H(z) to be causal, it must not have more zeroes than poles. Therefore, H(z) must not have more poles than zeroes.

If H(z) has linear phase, then it must have zeroes either in the unitary circle or simultaneously inside and outside of it. Therefore, 1/H(z) can't be stable if H(z) has linear phase.


What does nonnegative zero-phase response mean?

What they mean here is that the real-valued amplitude function of the linear phase FIR filter that you provide to the function must be non-negative, because it is interpreted as the desired squared magnitude of the minimum-phase filter. This amplitude function is referred to as zero-phase response, which is the frequency response of the filter when the coefficients are shifted such that they are centered around the origin, i.e., $b[n]=b[-n]$ is satisfied. Note that the filter length must be odd.
If $N$ is the (odd) filter length, the zero-phase response is given by
$$B_{zp}(e^{j\omega})=B(e^{j\omega})e^{j\omega(N-1)/2}$$
where $B(e^{j\omega})$ is the frequency response of the causal filter.


Why Coherence is not a valid metric when performing impact excitation?

a good reference is:

  Carter GC. Coherence and time delay estimation. Proceedings of the IEEE. 1987 Feb;75(2):236-55

The coherence is actually complex valued and bound by a magnitude of one.  I believe you are referring to the magnitude squared coherence as used in Carter's paper.  The term coherence by itself is used in optics.

The problem with the MSC is the probability density of the estimate of MSC, it takes a lot of independent samples to get a reasonably accurate estimate. When applied to a single sample, the confidence region is [0 1] at any type two error, which essentially implies that a single sample has to considered a deterministic signal and not the random model.  The following paper covers how to calculate the pdf of the sample MSC.

  T. Barnard, "Legendre polynomial expressions for the probability
  density function of magnitude-squared coherence estimates," in IEEE
  Transactions on Acoustics, Speech, and Signal Processing, vol. 29, no.
  1, pp. 107-108, February 1981. doi: 10.1109/TASSP.1981.1163516
  Abstract: Carter, Knapp, and Nuttall [1] presented a discrete Fourier
  transform (DFT) method for estimating the magnitude-squared coherence
  between two zero-mean wide-sense-stationary random processes. This
  paper shows that a simple Legendre polynomial expression can replace
  the hypergeometric function used in the probability density function
  expression of Carter et al. [1] (see also Fisher [2] and Goodman [3])
  for the coherence estimate obtained from nonoverlapped DFT's. A
  standard Legendre polynomial recursion provides a recursive method for
  evaluating the density expression. keywords: {Polynomials;Probability
  density function;Density functional theory;Frequency;Discrete Fourier
  transforms;Signal processing algorithms;Random
  processes;Instruments;Detectors;Distributed computing}, URL:
  http://ieeexplore.ieee.org.mutex.gmu.edu/stamp/stamp.jsp?tp=&arnumber=1163516&isnumber=26151

If you have the symbolic toolbox, you can use the following matlab code to calculate the MSC pdf using the expression from the Carter paper.

clear all
%nd=128
n=[ 1 32 512 768 ]
c=linspace(0,1,512);
p=zeros(length(n),length(c));
for k=1:length(n)
nd=n(k)
%syms  C Ch
C=.5;
p(k,:)=Magsquarecoherepdf(c,C,nd);
end
figure(1)
plot(c,p,'Linewidth',2)
h1=legend('1', '32','512', '768');
title(h1,'Number Trials');
title(['Sample Density given $\mid \gamma \mid^2 = $',num2str(C)],'FontWeight','bold','Interpreter','latex')
h=xlabel('Sample $\mid \tilde{\gamma} \mid^2$','FontWeight','bold','Interpreter','latex');

function [p]=Magsquarecoherepdf(c,C,nd)
%
p=zeros(size(c));
for i=1:length(c)
Ch=c(i);
t1=log10(hypergeom([1-nd 1-nd],1 , Ch*C));
assert(~isnan(t1),'t1 nan')
t2=log10((nd-1).*((1-Ch.*C)./(1-C).^2));
assert(~isnan(t2),'t2 nan');
t3=log10(((1-C).*(1-Ch)./(1-Ch.*C).^2))*nd;
assert(~isnan(t3),'t3 nan');
p(i)=10.^(t1+t2+t3);
%assert(~isnan(p(i)),['p  assert ',num2str(t1),' ', num2str(t2),' ', num2str(t3)])
end
return

If you repeat your impacts, you have an ensemble and I believe you could apply the MSC to your problem.

edit:

If you might permit an analogy. lets say you are asked to determine if a coin is fair in a toss. You can't say very much from a single flip of a coin. you have to repeat the toss , and keep on repeat tossing it. the number  of tosses required for a specific accuracy is determined from a probability law. a single toss tells you nothing.

the same is true for the MSC. a single trial says nothing. it needs to be repeated and the probability density in the code above specifies the number of trials required for  a specified accuracy


continuous chain of ones or zeros for transmit

A dominant concern with long runs of ones or zeros is the ability of the receiver to establish symbol timing. The receiver clock is independent from the transmitter clock, so without transitions in the data the receiver will not be able to determine the clock used in the transmitter for each symbol, and will ultimately drift off. Whitening the data has other benefits for equalization and overall efficiency.


Why are there only a few high frequency components in the Discrete Cosine Transform?

There are two parts to the answer.  One has to do with the nature of sharp edges, the second has to do with the how imaging systems (lenses and the human eye) work.

As @PeterK's example demonstrates artificial (human created) images tend to have large regions of one color with sharp transitions between the regions.  This is similar to a square wave.  But if you look at the fourier transform of a square wave it will be a decreasing sum of the odd frequencies:
$$
\sum_{n=1, 3, 5, ...}\frac{1}{n} \sin{ n \omega t}
$$

Images of the natural world also often have regions of slowly changing color with sharp transitions where one object overlaps another.  There can actually be a lot of high frequency "texture" in an image (grass, leaves, a pebbly beach, waves on water that is relatively far away, whispy clouds) but the sharp edges of the occlusions (overlaps) obey the same square-wave fall-off in frequencies beyond the fundamental frequency (the average distance between tree trunks, blades of grass, etc.)

The second part of the answer is that imaging systems (camera lenses/eyes) can not make out features smaller than a certain minimum angular resolution.  Images outside of the imaging systems depth of field will be blurred because lenses have focal points.  But even objects inside the depth of field will be blurred.  This is because the aperture through which the light must pass must be reasonably large (several mm) in order to let in enough light.  But this means that point-sources of light actually create a circle of confusion that essentially acts as a low-pass filter.

Because of this the human eye, for example, is more sensitive at some frequencies than others.  In the below image the contrast differences across horizontal scans in the image are the same, but your eye is more sensitive to mid-frequency contrast differences than high-frequency contrast differences.  (You may have to adjust the distance from the image to your eye to see the effect line up with the black line drawn on the image.)


Implication of Cramer Rao Lower bound and relation with SNR : non-technical overview

This means that FI is calculated based on the true, known value of the estimate.

No! The Fisher estimation is based on the expectation of the likelihood function. The likelihood is just a normal probability – just that the role of the "original" variable and observation have been reversed.

I don't understand your example – spheres don't have a CRLB,  the problem of estimating a certain parameter has.

Maybe you're confusing the "set of all estimators for $\theta$" (which has a CRLB) with "a specific estimator" (which simply has a variance) or a with "a specific estimate" (which has an error).
That is, the Fisher info is a function of the "original" $\theta$.

  So, Can Fisher Information be high and does it vary with the noise?

So imagine a binary source: there's only two possible $\theta$, and they are sent through a memoryless AWGN channel – that is, you observe $x$, and it's either caused by $\theta=0$ or $\theta=1$, and figuring that out is your estimator's job (which gives you $\hat\theta$).

The likeliness function $f(x;\theta)$ hence gives you an idea of how likely an outcome $x$ is explained by the original effect $\theta$ was.

Obviously, if there's very little noise, that function will have a very "sharp" shape – there's little chance a 0 sent will be interpreted as 1 in the receiver.

Increasing the noise power makes it possible that this actually happens – the likelihood changes, and it becomes much more probable that the other $\theta$ value is the reason for what you're seeing.


position of useful filtered data inside DFT based filtered output

Lets put forward the intuition behind the concept of the group delay before further discussing how to find the delay of FIR filters.

Consider an input signal $x[n]$ of length $L_x$ which is nonzero between $n=0$ and $n=L_x-1$. And let a simplistic filter impulse response to be $h[n] = \delta[n-d]$. The output is immediately shown to be $y[n] = x[n-d]$ which is input shifted by $d$ samples to the right (for $d$ positive).

Let's look at the frequency response of the filter; $H(e^{j\omega}) = e^{-j \omega d}$ which is nothing but a phase term $\phi(\omega) = - \omega d$.

Consider a slightly complex filter $h[n] = 0.1 \delta[n] + 0.8\delta[n-1] + 0.1 \delta[n-2]$. Its frequency response can be written as $H(e^{j\omega}) = e^{-j 1 \omega} B(\omega)$ where $B(\omega) = 0.2 \cos(\omega n) + 0.8$ is a real valued function and has no phase-term. The output of this filter can be written in terms of delayed inputs due to each distinct impulse such as :
$$y[n] = 0.1 x[n] + 0.8 x[n-1] + 0.1 x[n-2]$$

It can be proposed that $y[n]$ is dominated by the middle term $x[n-1]$ of a shift $1$ sample. (a more detailed analysis would say that low frequencies of $x[n]$ are dominated by the middle term and high frequencies are tending to cancel etc. which makes it a lowpass filter eventually)

In this latter example the phase function of the filter was $\phi(\omega) = -1\omega$ whereas the dominating component $x[n-1]$ at the output also had a delay of one sample to the left.

This observation can be used as an intuitive basis for the claim that the delay (or the group delay) associated with an FIR filter is given by the negative of the derivative of its phase response wrt frequency $\omega$; i.e.,
  $$\boxed{ \tau = - \frac{d \phi(\omega)}{d\omega} }$$

Lets restrict the discussion to LTI, FIR symmetric filters which have linear phase responses. An LTI filter is described by its impulse response $h[n]$ which enables us to compute its response to an arbitrary input $x[n]$ via the discrete convolution sum:

$$ y[n] = x[n] \star h[n] = \sum_{k=-\infty}^{\infty} x[k] h[n-k] $$

where the dummy index $k$ would actually run over those valid range of nonzero multiplicands.

The filter also has a frequency response which is the DTFT of its impulse response:
$$ H(e^{j\omega})=  \mathcal{F} \{ h[n] \} =\sum_{n=-\infty}^{\infty} h[n]e^{-j\omega} = |H(e^{j\omega})| e^{j \phi_h(\omega)}$$

Now we define the group delay (in samples) of filter as:
$$\tau = - \frac{d\phi_h(\omega)}{d\omega}$$

In general group delay will be a function of frequency $\omega$, however linear phase FIR filters have the property that their phase function is a linear function of frequency such as $$\phi_h(\omega) = K - d \omega$$

From which we deduce that the group delay will be $\tau = d$ which is independent of frequency; i.e. all components of an applied signal $x[n]$ will shift the same amount when filtered by a linear phase filter.

Given such a linear phase filter whose delay is $\tau = d$ and assume that you apply a very narrowband signal $x[n]$ such as:
$$ x[n] = \cos(\omega_0 n) \frac{ \sin(\omega_1 n)}{\pi n} $$
where $\omega_1 \ll 1$  and $0 < \omega_0 < \pi$.  Then the associated output for such a signal will be of the form
$$ y[n]=A(\omega_0) x[n-d] = A(\omega_0) \cos(\omega_0 (n-d)) \frac{ \sin(\omega_1 (n-d))}{\pi (n-d)}$$

Finally, for example,  it can be shown that a causal and symmetric FIR filter of length $M = 2L + 1$ from $n=0$ to $n=2L$ such that $h[n] = h[2L-n]$ will have a linear-phase term of the form $\phi_h(\omega) = K - \omega L$. Where $K$ can be zero. Hence its group delay will be $L$ samples. Lets show it here:

$$
\begin{align}
H(e^{j\omega}) &= \sum_{n=0}^{2L} h[n] e^{-j \omega n} \\
&= \sum_{n=0}^{L-1} h[n] e^{-j \omega n} + \sum_{n=L+1}^{2L} h[n] e^{-j \omega n} + h[L] e^{-j \omega L} \\
&= \sum_{n=0}^{L-1} h[n] e^{-j \omega n} + \sum_{n=L+1}^{2L} h[2L-n] e^{-j \omega n} + h[L] e^{-j \omega L} \\
\end{align}
$$
Substitute $2L-n = m$ in the second sum:
$$
\begin{align}
H(e^{j\omega}) &= \sum_{n=0}^{L-1} h[n] e^{-j \omega n} + \sum_{m=L-1}^{0} h[m] e^{-j \omega (2L-m)} + h[L] e^{-j \omega L} \\
&= \sum_{n=0}^{L-1} h[n] e^{-j \omega n} + e^{-j\omega 2L} \sum_{m=L-1}^{0} h[m] e^{j \omega m} + h[L] e^{-j \omega L} \\
&= \sum_{n=0}^{L-1} h[n] e^{-j \omega n} + e^{-j\omega 2L} \sum_{m=0}^{L-1} h[m] e^{j \omega m} + h[L] e^{-j \omega L} \\
\end{align}
$$
Recognize the similarity of the first sum and second sum such that denoting the first sum as $A(\omega)$ then the second sum is $A(\omega)^*$ (assuming real $h[n]$), Hence

$$
H(e^{j\omega}) = A(\omega) +  e^{-j\omega 2L} A(\omega)^* + h[L] e^{-j \omega L} \\
$$

Taking all of them into the $e^{-j\omega L}$ parenthesis yields:
$$
H(e^{j\omega}) = e^{-j \omega L} \left( A(\omega)e^{j \omega L} +  A(\omega)^* e^{-j \omega L} + h[L] \right) \\
$$

Since $A(\omega)e^{j \omega L} +  A(\omega)^* e^{-j \omega L} = 2 \mathcal{Re} \{ A(\omega)e^{j \omega L} \}$

We reduce the frequency response of the FIR filter to:
$$H(e^{j\omega}) =  e^{-j \omega L} ( 2 \mathcal{Re} \{ A(\omega)e^{j \omega L} \} + h[L] ) $$

Where the thing in the parenthesis is a real valued function and has zero phase term, therefore we call it $B(\omega)$ and deduce that the DTFT $H(e^{j\omega})$  of the causal, symmetric, real valued impulse response $h[n]$ is

$$ \boxed{ H(e^{j\omega}) = B(\omega) e^{-j\omega L} }$$

from which we find that the associated phase response of the filter is:
$$\phi_h(\omega) = -\omega L $$

and the associated group delay is therefore
$$\boxed{ \tau = - \frac{d \phi_h(\omega)} {d\omega} = L}$$

samples...

% This script demonstrates the use of time-aliasing in an advantageous way
% in obtaining the "central" portion of the circular convolution implemented
% via the DFT approach...
%
% Let x[n] of length Lx and h[n] of length Lh = 2*L+1, are convolved linearly
%
%                     y[n] = x[n] * h[n]
%
% so that "full" length of y[n] is:  Ly = Lx + Lh -1
%
% However in some cases we are only interested in those samples of length Lx
% that relate to input x[n] which are at the "central" portion of y[n], whereas
% the initial and final L samples carry transient information. Hence we take
%
%                    yc[n] = y[L:L+Lx-1]
%
% When a linear convolution is implemented with DFT, the length of DFT must
% satisfy the following; if we are to avoid time-aliasing: M >= Lz=Lx+Lh-1
% However this method computes ALL the full samples of y[n] in addition to the
% central ones. Hence if we wish to compute only the central portion of
% linear convolution y[n], then we can allow some aliasing in the DFT method
% The amount of which is indicated in the length of DFT as:
%
%                    M = Lz - L
%
% Then this will produce a result ya[n] whose last Lx samples carry the
% required information such that
%
%                    yc[n] = y[L:L+Lx-1] = ya[Lz-L-Lx:Lz-L-1]
%
% in matlab notation the indices will be added by 1 as:
%
%                    yc(:) = y(L+1:L+Lx) = ya(Lz-L-Lx+1:Lx-L);
%
% ==============================================================================
clc; clear all; close all;

% S0 - Define the parameters:
% ---------------------------
Lx = 32;
L = 15;
Lh = 2*L+1;
Lz = Lx+Lh-1;

% S1 - Generate the signals:
% --------------------------
x = randn(1,Lx);
h = randn(1,Lh);

% S2 - implement time-domain convolution first:
% ---------------------------------------------
yt = conv(x,h);
y1 = yt(L+1:L+Lx);           % get the central portion of the convolution

% S3 - Implement frequency domain DFT based aliased-circular conv:
% ----------------------------------------------------------------
M = Lz-L;       % DFT length (with aliasing)

yt = real( ifft( fft(x,M).*fft(h,M) , M) );     % compute aliased result
y2 = yt(L+1:end);

% S4 - compare the results
% ------------------------
figure,stem(y1,'b')
hold on
stem(y2,'g')

figure,stem(y1-y2);title('the error');


Your sequence is

$$x[n]=a^n\sin(n\omega_0)u[n]\tag{1}$$

First of all, note that the Discrete-Time Fourier Transform (DTFT) of $(1)$ only exists if $|a|<1$. (The case $|a|=1$ can be handled by using Delta impulses). Anyway, for $|a|>1$ the DTFT of $(1)$ does not exist. Second, you can write $x[n]$ as the multiplication of $x_1[n]=a^nu[n]$ and $x_2[n]=\sin(n\omega_0)$. The DTFTs of these two sequences are given by

$$\begin{align}X_1(e^{j\omega})&=\frac{1}{1-ae^{-j\omega}},\quad |a|<1\\
X_2(e^{j\omega})&=\frac{\pi}{j}\left[\delta(\omega-\omega_0)-\delta(\omega+\omega_0)\right],\quad -\pi <\omega <\pi\end{align}\tag{2}$$

Of course, $X_2(e^{j\omega})$ is also $2\pi$-periodic, but in $(2)$ only the interval $-\pi<\omega<\pi$ is considered (assuming $0<\omega_0<\pi$).
Note that in your expression for $X_2(e^{j\omega})$ you have a sign error. Furthermore, note that $X_2(e^{j\omega})$ is the DTFT of $x_2[n]=\sin(n\omega_0)$, and not of $x_2[n]=\sin(n\omega_0)u[n]$, as claimed in your question.

Now we have

$$X(e^{j\omega})=\frac{1}{2\pi}X_1(e^{j\omega})\star X_2(e^{j\omega})=\frac{1}{2\pi}\int_{-\pi}^{\pi}X_1(e^{j\theta})X_2(e^{j(\omega-\theta)})d\theta\tag{3}$$

The only thing you need to know to compute $(3)$ is the property

$$F(\omega)\star \delta(\omega-\omega_0)=F(\omega-\omega_0)\tag{4}$$

for any function $F(\omega)$. With $(4)$, evaluating $(3)$ with the DTFTs in $(2)$ results in

$$\begin{align}X(e^{j\omega})&=\frac{1}{2\pi}\frac{\pi}{j}\left[X_1\left(e^{j(\omega-\omega_0)}\right)-X_1\left(e^{j(\omega+\omega_0)}\right)\right]\\&=\frac{1}{2j}\left[\frac{1}{1-ae^{-j(\omega-\omega_0)}}-\frac{1}{1-ae^{-j(\omega+\omega_0)}}\right]\tag{5}\end{align}$$

The two terms in $(5)$ can be combined resulting in

$$X(e^{j\omega})=\frac{ae^{-j\omega}\sin(\omega_0)}{1-2a\cos(\omega_0)e^{-j\omega}+a^2e^{-2j\omega}}\tag{6}$$


This is easier to see graphically, on an s-plane, knowing that the phase is from the frequency response of the system and the frequency response is determined when we restrict s to be the $j\omega$ axis.

Consider the simple ratio of a single pole and a single zero expressed as

$$H(s) = \frac{s-z_1}{s-p_1}$$

Given s as any point on the s-plane, and $z_1$ and $p_1$ as specific points, we see that the above expression can be viewed as the ratio of two phasors, alternatively expressed as $$H(s) = \frac{A\angle\phi_1}{B\angle\phi_2}$$

To see this clearly, and how it will ultimately answer your question; first consider two arbitrary points on the complex plane, for example m and n, as shown in the graphic below:

The equation $m-n$ is a phasor with a magnitude and phase $K\angle \phi$. By subtracting n from m, we are shifting the origin to be at position n with relation to m. And thus we get a phasor with length K which is the cartesian distance between m and n, and the angle $\phi$ which is the angle of m-n with the origin at position n:

As stated initially, the phase contribution of each pole and zero is given by the frequency response, and the frequency response is given by restricting s to be the $j\omega$ (or imaginary) axis. So as we sweep s along the $j\omega$ axis, each zero in the numerator and each pole in the denominator contribute to the magnitude and phase of the overall response. The phase of the numerator will add, and the phase of the denominator will subtract (for example the simple case of 1/j = -j wich is $-\pi/2$ (-90°): As phasors, this expression is magnitude 1 at angle zero divided by magnitude 1 at angle $\pi/2$).

Thus the angle contribution for each pole and zero is as described in the graphic below. We see that as s approaches infinity, the angle contribution from the zero will approach +$\pi/2$, and the angle contribution from the pole (given it is in the denominator) will approach -$\pi/2$. So in this case at infinity the resulting angle will be angle 0, and at any given frequency we can solve for the angle analytically by using $s = j\omega$ for frequency $\omega$:

You can apply the above to your specific case, except $0.9e^{j\pi/4}$ is a pole in the right half plane, which would be unstable (will not converge). If you corrected that for a case where both poles are in the left half plane, then each pole will contribute phase as described above and the overall phase response will be the sum of the two. Thus for two poles in the left half plane, the phase will approach -$\pi$ (-180°) as s approaches infinity.