why is the z transform transfer function 1/(z-1) called an integrator?

There are a couple reasons.  One is that $(1-z^{-1})$ represents $x[n]-x[n-1]$ which is a finite difference over a very small period of time.  and that is an approximation to a differentiator.  The reciprocal is

$$ \frac{1}{1-z^{-1}} = \frac{z}{z-1} $$

which is the inverse operator.  We normally call the inverse operation of differentiation, we call that "integration".

Another reason is simply to implement that term as a transfer function of a tiny little LTI system:

$$ \frac{Y(z)}{X(z)} = \frac{1}{z-1} = \frac{z^{-1}}{1-z^{-1}} $$

or

$$ Y(z)(1 - z^{-1}) = Y(z) - Y(z) z^{-1} = X(z) z^{-1} $$

that translates to

$$ y[n] - y[n-1] = x[n-1] $$

or

$$ y[n] = y[n-1] + x[n-1] $$

so the current output sample is the previous output added to the (slightly delayed) input.  the output is an accumulation of the input.  similarly an integral of an input $x(t)$ is an accumulation of that input until the present time $t$.

i would say a better representation of a discrete-time integrator is $\frac{z}{z-1}$.  that would correspond to

$$ y[n] = y[n-1] + x[n] $$


When designing a filter, you really care about its behavior in two regions:

Passband: You want little attenuation in this region, and maybe other properties as well, like linear phase, depending upon your application.
Stopband: You want as much attenuation as needed in this region.

Between these two is the transition region. This is treated as somewhat of a "don't-care" band. You don't typically constrain the response too tightly in this area, so you can't really count it being usable. In your example, the passband lies below $0.45 f_s$; after the interpolator, you only plan on using frequency content below this threshold.

This means that you can allow some aliasing in order to simplify your filter design. Your transition region starts at $0.45 f_s$; everything above that frequency in your filter's output will either have a response that is unpredictable given your filter specs (if it lies in the transition band) or one that is highly attenuated (if it lies in the stopband). The takeaway:

  You can't rely upon any frequency content past the passband edge in your filter output anyway. So if it makes the filter design cheaper, why not allow the frequencies above the passband edge in the filter output to contain aliased garbage?

This technique is commonly used in multirate filters as you've noticed, as it allows savings in the required filter order in order to meet a given set of passband/stopband specifications.


Why do we follow Shannon's Theorem in audio sampling?

What your professor probably meant is that in most audio converters a small amount of aliasing is permitted (for example if you put 24.1khz into a 48khz converter you will see an output component at 23.9kHz that is not very heavily attenuated). The reason for this is that manufacturers like to use so- called half-band filters in the decimation chain because this saves a factor of 2 in the number of multiplications. Halfband filters impose a response symmetry around FS_out/2 so the attenuation at Nyquist is only 6db.  Since this is above the audio band it is not audible. However if you attempt to use this shortcut at lower sample rates like 8khz, it can be quite audible so it is not a recommended technique for lower sample-rates.

Bob


Yes, they can be different, depending what your maximum input value is.


The decimation operation, under some acceptations, consists in low-pass filtering the data before downsampling it. The standard in Python scipy.signal.decimate, in Matlab decimate: Filtering Before Downsampling, or even in R [decimate] (but it is borrowed from Octave), is to use a Chebyshev type I (of order 8, 10 and 8  respectively).
In Matlab cheby1, or at Comparison of Classical IIR Filter Types
, one can find comparative plots, such as:

Basically, Chebyshev filters aim at improving lowpass performance by allowing ripples in either the lowpass-band (Type I) or the highpass-band (Type II), whereas the behavior is monotonic in the complementary band.
With Type I, you are ensured that, if two frequency components in the highpass-band have the same amplitude, the highest frequency amplitude will be lower after filtering, which can be comforting when downsampling afterward.
With Type II, you have the converse in the lowpass-band.
Personally (a kind of coming-out, and I'd be glad to be corrected by my colleagues):

if my data has precise high-frequency limit (esp. for lectures), and I'm pretty sure of the cutoff, I tend to use type I
if my data is quite imprecise, noisy, and I care more for relative amplitudes in the low-pass band, I tend to use type II (empirically) to reject the noise, even if I lose sharpness.

Edit from r b-j:
Note in this modified plot, the cheby2 filter had its stopband frequency adjusted so that the -3.01 dB cutoff frequency coincides with the other three filters (and I normalized it to 1).  Note that the 5th-order Type I and Type II Tchebyshev filters both transition from -3.01 dB (the edge of the passband) to -30 dB (the edge of the stopband) in exactly the same transition band from a normalized frequency of 1 (passband) to 1.364 (stopband).
So the Type I and Type II Tchebyshev filters have the same performance in terms of sharpness of cutoff, if you compare apples to apples.  The main tradeoff between Type I and Type II is where the ripple is.  In Type I, the ripple is in the passband (where I usually would rather not see it) and in Type II the ripple is in the stop band (where I usually don't care).

MATLAB code:
clear;
close all;

N = 5;
f0 = 1;
w0 = 2*pi*f0;

passband_gain = 10*log10(0.5);          % -3.01 dB (half power)
stopband_gain = -30;                    % dB

% Design a 5th-order analog Butterworth lowpass filter with a cutoff
% frequency of 1 Hz. Multiply by  to convert the frequency to radians per
% second. Compute the frequency response of the filter at 4096 points.

[zb,pb,kb] = butter(N, w0, 's');
[bb,ab] = zp2tf(zb,pb,kb);
[hb,wb] = freqs(bb,ab,4096);

% Design a 5th-order Chebyshev Type I filter with the same edge frequency
% and 3.01 dB of passband ripple.

[z1,p1,k1] = cheby1(N, -passband_gain, w0, 's');
[b1,a1] = zp2tf(z1,p1,k1);
[h1,w1] = freqs(b1,a1,4096);

% Design a 5th-order Chebyshev Type II filter with the passband edge
% frequency and 30 dB of stopband attenuation.

w_stopband = w0*cosh((1/N)*acosh(sqrt(10^(-stopband_gain/10)-1)));
[z2,p2,k2] = cheby2(N, -stopband_gain, w_stopband, 's');
[b2,a2] = zp2tf(z2,p2,k2);
[h2,w2] = freqs(b2,a2,4096);

% Design a 5th-order elliptic filter with the same edge frequency, 3.01 dB
% of passband ripple, and 30 dB of stopband attenuation.

[ze,pe,ke] = ellip(N, -passband_gain, -stopband_gain, w0,'s');
[be,ae] = zp2tf(ze,pe,ke);
[he,we] = freqs(be,ae,4096);

% Plot the gain in decibels. Compare the filters.

plot(wb/(2*pi),mag2db(abs(hb)))
hold on
plot(w1/(2*pi),mag2db(abs(h1)))
plot(w2/(2*pi),mag2db(abs(h2)))
plot(we/(2*pi),mag2db(abs(he)))
hold off
axis([0 3 -40 5])
grid
xlabel('Frequency')
ylabel('Gain (dB)')
legend('butter','cheby1','cheby2','ellip')


Why the time step series affect the noise in FFT results?

Just use MATLAB's interp1 function before FFT


Why must the order of a Band-Pass and Notch filter always be even?

In principle there is no reason why the filter order of a general bandpass or bandstop filter must be even. Such a restriction is a consequence of a specific design procedure. In classic IIR filter design (Butterworth, Chebyshev, Cauer) you start with an analog prototype lowpass filter. Bandpass or bandstop filters are then obtained by a frequency transformation. And it is this frequency transformation that doubles the order of the prototype lowpass filter, hence the even filter order for that specific design method. Note that the bilinear transform has nothing to do with that restriction on the filter order.
Of course, there are specific designs for which an odd filter order doesn't make much sense. E.g., for a notch filter (with a notch frequency greater than zero) you want exactly one zero at the notch frequency plus its mirror image at the negative notch frequency. So for each notch frequency you get two zeros, and, consequently, two poles.


Estimating an Image Filter from Examples

I will talk about the Linear Spatially Invariant (LSI) Model.

Since the Filter (System) is LSI it can be described in the domain of Fourier Transform:

$$ Y \left( \omega \right) = H \left( \omega \right) X \left( \omega \right) $$

Where $ Y \left( \omega \right) $ is the Blurred (Output) Image, $ X \left( \omega \right) $ is the Input Image and $ H \left( \omega \right) $ is the filter model.

Now, assuming the filter doesn't vanish in the Fourier Domain, then:

$$ H \left( \omega \right) = \frac{Y \left( \omega \right)}{X \left( \omega \right)} $$

Now, to practical issues with this model:

Noise.
Vanishing Frequencies.

There are many methods to deal with those issues (Pseudo Inverse, Wiener Filter, Lucy Richardson) all are having different assumptions (White Noise, Gaussian Noise, Linear Model, Poisson Noise, etc...).

Till now I described the classic "Signal Processing" Model.

What you can do is build this model in the Spatial Domain.
The Filtering Operator will be a very large Toeplitz Matrix which you can estimate using LS Model or any solver for any cost function you'd want.

Practically in the world of Image Processing, the quality is in the regularization yo enforce.

P.S.
There is a new community proposal in the SE network:

http://area51.stackexchange.com/proposals/86832

It would be great if you could follow, up vote question with less than 10 points and propose your own questions.

Thank You.


Constraints for Type 3 Linear Phase Filter

There is no such restriction for a type III linear phase filter. Neither is there such a restriction for any other type of linear phase filter. Type III linear phase filters have an odd filter length $N$ and odd symmetry. I.e., as long as the filter's impulse response $h[n]$ satisfies

$$h\left[\frac{N-1}{2}+n\right]=-h\left[\frac{N-1}{2}-n\right]$$

you have a proper type III linear phase filter.

If you look at the impulse responses of type III differentiators or Hilbert transformers then you'll notice that in both cases the filter taps are decreasing in magnitude as you move away from the center tap. So these common applications of type III filters are counterexamples to your assumption.


Confusion regarding calculation of fundamental period?

This seems like more of a semantics problem.

A signal is periodic with time $T$ if

$$x(t+n\cdot T) = x(t), n \in \mathbb{Z}$$

So the signal is periodic in $0.5$ since the for $T = 0.5 \cdot n$ the argument of the cosine is an integer multiple of $2 \pi$. Since it's periodic in $0.5$ it's also periodic in all integer multiples of $0.5$, i.e $1$, $1.5$, $2$ etc.

In this case it's also periodic in $0.25$ since
$$ \cos^2(4 \cdot \pi \cdot t ) = 0.5 \cdot (1+\cos(8 \cdot \pi \cdot t))$$

So any periodic signal has an infinite number of periods, the fundamental one is the smallest one and all the others are integer multiples of the fundamental.


Why removing noise increases my audio file size?

This is has nothing to do with the noise removal process but with your encoding. In order to do any type of processing with an encoded audio file you need to

Decode the input file to PCM
Process PCM data
Encode the PCM into the desired output format

Please note that steps and 1 and 3 are completely independent. The decode is determined by the input file format and the encode by the desired output file format. If you want those to be the same, you need to capture all encoding parameters from the input file and pass them to the encoding process.
In your case, you are probably using a different bitrate for input and output.


Before applying the algorithm below, you should first apply a DC blocking filter. That should remove the DC offset (and partially the step at the beginning), which will allow more usual periodicity techniques to be used.

You can try using the average square difference function. There is a nice write up on this answer and a way to implement it in an efficient manner here.

Here's an example plot. The $\color{red}{\bf red}$ circle shows the selected minimum. Note how the "octave" problem shows up: if the signal is periodic with period $P$, it's also periodic with period $N P$ where $N$ is an integer.

R Code Below

#30511

pulse <- 16 - seq(-4,4)^2
len <- length(pulse)

T <- 1000
P <- 49

ix <- seq(1,T - len, P)

data <- rep(0,T)

for (this_one in ix)
{
  data[this_one + seq(0,len-1)] <- pulse
}

data <- data #+ rnorm(T,0,1)

plot(data, type="l")

asdf <- function(data, Pmax)
{
  N <- length(data)
  asdf <- rep(0,Pmax)
  for (p in seq(1,Pmax))
  {
    ix <- seq(p+1,N)
    asdf[p] <- sum((data[1:length(ix)] - data[ix])^2)
  }

  return(asdf)
}

calc <- asdf(data, 5*P)

par(mfrow=c(2,1))
plot(data, type="l")
title('Repeating pulses')

plot(1:(5*P), calc, type="l")
period <- which.min(calc)
points(P, calc[P], col="grey", lwd = 12)
points(period, calc[period], col="red", lwd = 5)
title('ASDF')


How to explain IDFT of square root of a DFT?

Since $V(\omega)V(\omega)=W(\omega)$, in the time domain we expect $v[n]*v[n]=w[n]$. That is, the convolution of $v[n]$ with itself should become $w[n]$.

$w[n]$ is in the form of a triangular signal. A triangular signal of odd length $N$ can be constructed by convolving a half-rectangular signal of length $(N+1)/2$ with itself. So $v[n]$ is a rectangular pulse with $3$ elements equal to one and the remaining $2$ elements equal to zero.


Linear systems: Square root of input product


STFT - Plotting is not a true representation

Regarding your first 2 plots: for FFT results, the X axis is usually labeled with the FFT_bin_index * sample_rate / FFT_length, and the Y axis with either the magnitude (absolute value of the complex bin value), or the log magnitude of the FFT result bin.


Calculation of cross-correlation of a system with noise

\begin{align}
R_{W,Y}[n] &= E\big[W[m]Y[m+n]\big]\\
&= E\big[W[m](Z[m+n] + W[m+n])\big]\\
&= R_{W,Z}[n] + R_{W,W}[n]\\
&= R_{W,Z}[n]+ P_w\delta[n].
\end{align}
So the question is: can you figure out what
$R_{W,Z}[n] = R_{W,h\star X}[n]$ is where you know that $\{X[n]\}$ and $\{W[n]\}$ are uncorrelated processes, that is, $R_{W,X}[n] = 0$ for all $n$.
Hint: $Z[n]$ is a weighted linear combination of the $X[i]$, and, as Marcus Muller's comment points out, the curious fact about linear operations is that they are linear!


From IQ signal to FM modulated carrier, how its done?

The solution to this is realizing the relationship between frequency and phase:
$$f(t) = \frac{d\phi(t)}{dt}$$
Where $f(t)$ is the instantaneous frequency as a function of time, and $\phi(t)$ is the instantaneous phase as a function of time. So by such definition, a linear change in phase represents a step change in frequency. This is exactly how we can use vector modulators as frequency translators (and equivalently single-sideband modulators). A rotating phasor on the IQ plane, if rotating at a constant rate represents either a positive or negative frequency depending on the direction of rotation. So with FSK using two frequency symbols we get positive and negative rotating phasors. The frequency realized is the rate of that rotation (cycles/sec).
Starting with MSK, which is a subset of FSK, and specifically the minimum frequency shift between the two FSK tones such that the two signals are orthogonal, for this we rotate the phase 90° over one symbol duration, such that the direction of rotation determines if a one or zero was transmitted. The reason for 90° is this results in the minimum spacing of the two frequency tones as first described (hence MSK). If we rotated further within the same time of a symbol period this would result in a higher frequency and thus increase the spacing between the tones (as in generalized FSK). Below shows the phase diagram of rotation versus symbol "1" or "0". If we continued to send a long series of either symbol, it is quite clear that we are transmitting either a positive or negative frequency (relative to a central carrier frequency, so basically one of two possible tones: FSK).

When we transmit a series of symbols the frequency shift keying that is occurring is clear from the relationship between phase and frequency:

The sharp transition in frequency as a step is what creates a high spectral occupancy (sharp edges create high frequencies, consider all the spectral components of a square wave). By rounding these edges we minimize the spectral occupancy, and Gaussian Minimum Shift Keying specifically replaces the straight phase slopes with a Gaussian filtered transition; the phase specifically transitions as an integrated Gaussian pulse (thus the frequency follows a Gaussian transition- see diagram below). So the phase starts out slowly at the beginning of the symbol, reaches maximum rate mid way and then terminates slowly at the end where a possible otherwise abrupt direction change could occur. You can make a very simple (no multipliers so can be done in a micro-controller) GMSK modulator using a Gaussian Filter on the Frequency Control Word input to an NCO which is otherwise connected to the data stream (and similarly with a VCO in the analog world). (This is actually what I was getting at with this question: Simplest All Digital GMSK Modulator which can be combined with this solution for the simple no-multiplier Gaussian filter: Gaussian FIR filter with no multipliers?).  The GMSK phase and frequency diagrams below further and most immediately shows how this is a form of FM modulation.

GMSK is also commonly implemented as partial response signals where we overlap consecutive pulses with intentional inter-symbol interference which can later be completely undone (with Viterbi decoding for example), and thus further improves spectral efficiency since we can send more data in the same amount of time (at the expense of receiver complexity).  The strategic use of the wonderful Laurent Decompositions significantly simplifies those receivers— but that is another topic for another day.
Further FSK, MSK, GMSK are all constant envelope modulations so have the noted advantage that we can run the transmitter PA into saturation which offers significant power efficiency and without regard to immunity against intention jammers can hard limit the input instead of using an AGC (both of particular interest to low cost battery operated devices)


Hardware advice for machine vision project

The things you can control with the camera selection are:

Resolution
Sensor size
Sensitivity
Frame rate

The things you can control with the lens selection are:

Focal length
Minimum depth-of-field (DoF)
Field of view (FoV)
Minimum aperture
Zoom functionality

Each of these interact with each other so it is often difficult to find the right trade off.

The sensor size and resolution are usually related. For example, it is hard to find sensors with resolution above 1024 * 1200 in a 1/3" sensor.

The resolution and thus sensor size will dictate what lenses can be used with the camera, then combined with the FoV of the lens will dictate how many pixels the face covers in the image.

First decide how far off-centre the face can move. That will guide you minimum FoV. The work out how far away the face can be, and how many pixels wide it must cover. Using the average width of a face, the FoV angle, and the depth you can work out the minimum resolution you need.

Using the minimum resolution and DoF, choose a camera and lens combination. Pay close attention to the maximum sensor size the lens can accomodate, and note that the FoV changes with the sensor size.

The DoF is a function of the lens focal length and the aperture. To increase the DoF you have to reduce the aperture. E.g. f/1.8 is wide open, f/22 is almost closed. A smaller aperture will let less light in, and thus the camera exposure and gain will need to be increased. Choosing a camera with high sensitivity will allow you to get better images with less gain increase while keeping the exposure short enough to not be blurry.

Motor driven zoom lenses are available.

Near-IR is possible with many camera manufacturers, as they actually have to put in a filter to block IR, so making a near-IR one is easy - just take the filter out.

But probably my best piece of advice is to contact the camera manufacturer and ask them. They have people that will take your specifications and recommend a camera / lens pair for you. Personally, I have used Point Grey cameras.


How is compression ratio related to number of dct coefficients

Compression ratio is defined as the amount of bits in the original image divided by the amount of bits in the compressed image. In the original image the number of coefficients is 64, while you choose 15 coefficients after DCT in the compressed image with mask A. The ratio is 64*6/15*8. With mask B it is 64*6/41*8


The formulation for MEL frequency scale is logarithm scale based, it clubs(put together) the higher frequencies more closely(so less resolution), because logarithmic functions vary slowly for larger numbers(see the logarihmic functions plot).

So some function of exponential form will definitely work out if one wants to give more resolution to higher frequencies.

so based on particular requirements one can define exponential scale for effect that looks like reverse of MEL scale.

Following is the filter bank I created using exponential function, similar to the formulation of MEL scale.

$fnew$ $=$ $exp(f / 7000)$.

I think this is what you are looking for, more resolution for higher frequencies as compared to lower frequencies.


Modeling a signal through the environment

You can model frequency shift this way:

$s = s \circ exp(j \cdot 2\pi \cdot (f-\hat{f}) \cdot t)$,

where $\circ$ means Hadamard product and $s$ is input signal vector (e.g. frame or part of data stream of interest). You can implement this operation in the baseband or in carrier frequency. The effect will be generally the same.

As you've told your signal time can increas if multipath exists. You can model it as FIR filter of the order of your impulse response length. Refer to the channel modeling techniques for details.

Also your signal time can vary (increase and decrease) if relative motion of receiver and transmitter exists. This effect appears in the case of so called wideband model, e.g. sonic channel in the underwater communications. For such a case signal bandwidth is comparable to carrier frequency so Doppler shift is very different for low frequencies and high ones in the band of interest. It causes varying of signal duration in the receiver end. If you deal with radio signals your signal model is probably narrowband so this effect is neglectable.


Is there a need of Point interpolation before proceeding for gaussian smoothing of an incomplete distribution?

Irregularly sampled data need not necessarily be interpolated prior to smoothing. You can try to smooth your data in the following way:

$$y(t_n)=\frac{\sum_{k}g(t_n-t_k)x(t_k)}{\sum_{k}g(t_n-t_k)}\tag{1}$$

where $t_n$ are your sample times, $x(t_n)$ are your data, and $g(n)$ is your smoothing kernel. Note that $t_n$ need to be integers (as in your example), which is no restriction, just a matter of formulation. Since $t_n$ are irregular, you do not always use all values of the chosen kernel $g(n)$, but this is compensated for by the normalization in (1). The kernel width has to be chosen such that the kernel always covers at least several of your irregularly sampled data. The smoothing in (1) is also used in the Nadaraya-Watson estimator.


Help me understanding this source code of a Phaser

_zm1 is the internal state variable of the allpass. The allpass is an IIR filter and hence it's current output is dependent on past inputs and outputs. This is the way to keep track of this.

In this particular case it's implemented using a set of coupled difference equation in so-called "Transposed Form II" filter topology. So $zm_1[n] = y[n]*a_1+x[n]*b_1 $ . For a first order allpass we have $b_1 = 1$


You might try Vincent Falco's Collection of Useful C++ Classes for Digital Signal Processing.

The StackOverflow Question A C++ library for IIR filter offers a few more suggestions.

Finally: the hard part for Butterworth/Chebychev filters is really the design, not the implementation.  You might consider doing the design in Matlab (or Gnu Octave), or an online site like: http://www-users.cs.york.ac.uk/~fisher/mkfilter/.  Writing the C code should be trivial once you know the coefficients for the recurrence relation.


How to understand the implementation codes in Debevec and Malik's HDR paper

OK, this is incomplete but this is how far I've gotten:

My previous answer, which you can see in the edits was incorrect.

Your question boils down to:

  However, one thing I am confused is that why the first part and the second part of the object function are treated independently.

So let's try having a look at the objective function:
$$
\begin{array}{rcl}
{\cal O} &=& \displaystyle \sum_{i=1}^{N} \sum_{j=1}^P \left\{
w(Z_{ij}) [ g(Z_{ij}) - \ln E_i - \ln \Delta t_j ]
\right\}^2 + \lambda \sum_{z =Z_{\rm min}+1}^{Z_{\rm max}} [w(z) g^"(z)]^2\\
&=& \displaystyle \sum_{i=1}^{N} \sum_{j=1}^P w(Z_{ij})^2
 [ g(Z_{ij})^2 + (\ln E_i)^2 + (\ln \Delta t_j)^2
- 2 g(Z_{ij}) \ln E_i - 2 g(Z_{ij}) \ln \Delta t_j + 2 \ln E_i \ln \Delta t_j
 ]\\
&&+ \lambda \displaystyle \sum_{z =Z_{\rm min}+1}^{Z_{\rm max}} [w(z) \{ g(z-1) - 2g(z) + g(z + 1) \} ]^2
\end{array}
$$
where the second line just expands the first term and substitutes for $g^"$ as the paper suggests.

Now, we need to choose the $\ln E_i$ and $g$ values to minimize this so we must take the derivative of $\cal O$ with respect to $\ln E_i$ and to $g$.

$$
\frac{d \cal O}{d \ln E_i} = \sum_{i=1}^{N} \sum_{j=1}^P w(Z_{ij})^2
 [  2 \ln E_i
- 2 g(Z_{ij}) + 2 \Delta t_j
 ]
$$


QAM Constellation Construction - Basic 4/8 QAM

You can choose any mapping that you desire from your unique data sequences to the symbol locations on the QAM constellation. It is common and usually best to Gray-code the mapping such that only one bit changes between adjacent constellation locations (those with minimum Euclidean distance). This way if there is a symbol error due to noise crossing an immediate decision threshold, there will only be one bit error.

See figure below showing the common mapping for QAM using Gray Coding from Wikipedia under "Gray Coding". Notice that the adjacent constellation points only differ by one bit. Under AWGN conditions at threshold error conditions with equi-probable symbols only one bit error would result with this mapping approach (the optimum error boundaries with equi-probable constellation points are mid way between the points shown).

For Gray-coded QPSK the mapping is quite trivial as illustrated by the figure from this link:

And here is a link to a simulation showing the significance in BER performance when Gray-coding QPSK.

The amplitude and phase for each possible symbol is directly taken from the constellation's I and Q values as shown in the diagrams above, and as traditionally done for any complex I, Q pair:

The amplitude is

$$A= \sqrt{I^2+Q^2}$$

And the phase is:

$$\phi = \mathrm{atan2}(Q,I)$$


Because when a quantity can be complex, and even when it is just real, the absolute squared difference $|f-g|^2$ can be expressed in both domain (complex and real) as:
$$|f-g|^2 = (f-g)^H(f-g)$$
and of course this is correct as well for reals. This setting is often related to Hilbert spaces.


What's the right scale to represent a frequency spectrum?

FFT scaling for audio signals is indeed difficult. Sometimes the energy is distributed over many frequencies, sometimes the energy is only at few or just one frequency. This can result in a very high "crest factor". i.e the ratio of the peak value to the average power (RMS) value.

This mainly depends on the the length of the FFT: the longer the FFT, the larger the crest factor will get.

The best ways to deal with this is:

Use floating point representation and math if any possible
Choose an FFT length that is optimal for your application and not longer than you really need

Other than that, the differences that you see are real. FFT is a metric of "spectral density", and for pure sine waves this approaches infinity, while it is quite low for noise-like signals.


Software PLL tracking of carrier frequency in bandlimited transmission

Here is a link to an example carrier tracking implementation for a QPSK and QAM system, applicable to full-response raised cosine pulse shaping. This post also has links to many other concerns with demodulating the waveform, including timing recovery and setting the parameters in the tracking loop implementation.

Recovering signal for psk

A x4 carrier recovery approach may be the simplest way to see how the pulse shaping effects the carrier recovery. For QPSK without pulse shaping you can recover the carrier by raising the waveform to the fourth power, and then dividing the resulting 4x carrier by four to complete the recovery. This is because phases add when you multiply, so if we have four possible phases in our modulated signal:

$$e^{(j\omega t + 0)}$$
$$e^{(j\omega t + \pi/2)}$$
$$e^{(j\omega t + \pi)}$$
$$e^{(j\omega t + 3\pi/2)}$$

Raising to the 4th power results in the following for all the cases above.
$$e^{(j4\omega t + 0)}$$

If we observe our waveform over time and it instantly goes from one phase location to the other (no pulse shaping, infinite bandwidth), then the 4th power result will be a pure tone at 4x the carrier frequency with no other components.

However if we slowly transition from one phase to the other, we will still get a pure tone at 4x the carrier frequency, but there will be significant energy in sidebands close to this frequency as well. The 4x tone will still be the strongest in vicinity of this 4x signal and is easily tracked with a PLL to complete recovery. (This is evidenced by a positive mean in the baseband equivalent analytic signal).  It is because this signal is the strongest that a PLL can track it (typically if the signal being tracked is +6 dB higher than the sidebands a PLL can lock to it), and it is because of the additional sideband energy that we need to the PLL to complete the recovery, effectively acting as narrow bandpass filter to provide a clean carrier.

I added some graphics below for the simpler case of BPSK with and without pulse-shaping that will hopefully provide an intuitive answer. Personally I find it much simpler to think and work with describing frequency using the complex exponential ($e^{j\omega t})$ rather than with cosines and sines, and for the same reason I avoid including an actual carrier but model at the baseband equivalent analytic signal. However I think in this case showing an actual sinusoidal carrier will more immediately get across the point I am trying to make. For the same reason, this is very illustrative in showing how the pulse shaping affects a tracking loop but it is not how I would typically do carrier recovery (see the link above for a typical software recovery loop I would implement).

This first image shows a sinusoidal carrier with an example (1 0 1) data pattern. Below the data pattern is shown the result of BPSK modulation with no pulse shaping applied (infinite bandwidth available).  Observe in the third plot showing the modulated signal, that the phase is changing from 0 to $\pi$ radians as the data transitions. This abrupt change back and forth of the phase would prohibit a simple phase tracking PLL from tracking the phase: if it could even lock onto one phase state (because the data doesn't transition for some time for example), once the data transition occurs the very significant phase error that would occur would cause a simple tracking loop to break lock (unless the phase detector had 0°/180° ambiguity -- hint, hint at another recovery approach!). Another way of describing this, with a random data pattern that is 50% ones and 50% zeros, the carrier would be completely suppressed and therefore does not exist, so there is no signal at that frequency for any loop to lock onto. Notice too that in the example plot I give, the data transitions are not commensurate with the carrier. I did that purposely just to show that regardless of the data rate and where transitions occur, squaring this signal will result in a pure tone (when no pulse shaping or other bandwidth limiting is used) at exactly twice the rate of the carrier.

The fourth plot is simply the waveform in the third plot squared with no other filtering applied. This waveform is completely in phase with the original carrier, at exactly twice the frequency.  In this case we could simply pass this signal into a frequency divide by two function (quite simple) and we would completely recover the carrier. No PLL is needed, but if there were any other noise sources present, a PLL would help increase the SNR of our recovered carrier in the presence of noise. However our waveform itself and the recovery process, in this case, contributes no additional noise or distortion.

What is occurring by squaring the signal is we are doubling the frequency (and the phase). This is clear from the trigonometric relationship:

$$cos(\alpha)*cos(\beta) = \frac{1}{2}cos(\alpha+\beta)+ \frac{1}{2}cos(\alpha-\beta)$$

$$cos(2\pi f_c t + \phi)cos(2\pi f_c t + \phi) = \frac{1}{2}cos(4\pi f_c t + 2 \phi)+ \frac{1}{2}cos(0)$$

$$ = \frac{1}{2}cos(4\pi f_c t + 2 \phi)+ \frac{1}{2}$$

So if the phase is only in two states: 0 (0°) and $\pi$ (180°), doubling this results in 0 and 2$\pi$ which is also 0.

Now observe what occurs in the same situation with pulse shaping as shown in the plot below. Here the only difference is a pulse shape has been added to each symbol. Notice that the result is to only add an amplitude variation (the pulse shape) to the squared signal but it has no effect on the phase. This is an important observation as to why a normal tracking loop could still track the phase in this situation, and also illustrates why the loop that is used should be able to "fly-wheel" in the absence of updates (as most tracking loops in both carrier and timing recovery). The phase in the doubled signal is in phase with the carrier at exactly twice the frequency but the amplitude is going up and down. This does of course impact the SNR of the recovered signal, as there are fewer phase error updates for tracking but unlike the case of the modulated signal itself, the phase does not abruptly change in a manner that would be disastrous for the loop to track. This also shows the importance of using a phase tracking loop in this case, which will ride through the amplitude variations and create a clean doubled carrier with no amplitude variation.

This is perhaps clearer by comparing the two final cases of the multiplied outputs directly as shown below. Here we clearly see that the pulse shaping only affects the amplitude and not the phase, hence a phase tracking loop should be able to track and recreate a cleaner replica of the doubled carrier signal, which we would then follow with a divide by two to completely recover the carrier.

The plots above should provide a satisfactory intuitive explanation in the effects of pulse shaping on carrier tracking loops, but not to suggest the favored approach to doing carrier recovery. In a digital or software implementation instead of squaring, PLL and divider, I would suggest decision directed and other similar digital Costas-Loop approaches as I have shown in the link above.

Another intuitive explanation is given by observing the above signals in the frequency domain, and noting how the carrier is suppressed prior to squaring and then recreated at twice the carrier frequency after squaring (or raising to the fourth power for QPSK, to the eight power for 8-PSK etc).


The Hann Window Kernel is the FFT of the Hann Window.

Offset is defined strangely, and unnecessarily. calculating InstFreqs and then dividing it by Sampling Frequency Fs is just throwing away the work of attaining it.

Bin Resolution = Fs / N
//InstFreqs will start at 0Hz and peak at Fs Hz (or half that??)
InstFreqs(i) = i * Fs / N

So in the referenced answer:

Offset = abs(InstFreqs(i) / (Fs/N) - i)

but if I'm correct that: InstFreqs(i) = i * Fs / N then:

Offset = abs(InstFreqs(i) / (Fs/N) - i) = abs( i - i ) = 0

So that line didn't make sense. Too bad there was no comment hinting at what it was supposed to mean.

index is put into Wk it must span the length of the Window, and not the length of the FFT. So let's start by changing that. N * Offset will (likely) exceed the length of M (window size)

index = floor(abs(Offset) * N) + 1
InstMag = Mag(i) * Wk(index)

Unless of course I misinterpreted what InstFreqs(i) is, which it seems I have.


What does it mean for audio to be 'sped up'? How can I undo it?

It could mean:

That the metadata in the .wav file wrongly indicates that the signal was recorded at a sample rate 10x higher than the true value (for example it is advertised as having been recorded at 80kHz while it has been recorded at 8kHz). Simply fixing the metadata undoes this.
That the signal has been downsampled by a factor of 10, but is still tagged as having been recorded at the original sample rate. This can be undone by dividing by 10 the playback sample rate (by tweaking the metadata); or by upsampling the signal by a factor of 10 without changing the playback sample rate. Note that in any case, 90% of the bandwidth of the original signal (the high frequency content) has been lost and cannot be restored (or only very approximately through spectral band replication).
That a time-stretching (time-scale modification) algorithm has been applied to it. This is a lossy process, and cannot be undone without severe loss at such a high stretch factor. The "antidote" is simply to use any time-stretching algorithm with a stretch factor of 10x.


Modelling transmitter and receiver impairments in communication system

It makes sense to me. If the ideal transmitted signal is $s(t)$, the signal that is actually transmitted is $s(t)$ plus noise and distortion. This is the signal that suffers attenuation and fading, and to which an extra heaping of noise and distortion is added in the receiver. Modeling transmitter imperfections is not something that is included in most introductory textbooks, but it's very important to implement a system with good performance. You may see for example "RF System Design of Transceivers for Wireless Communications" by Qizheng Gu for more.

Having said that, modeling the transmitter imperfections as just Gaussian noise is a bit optimistic, but it may be a step in the right direction.


How the white and colored noise differ in time domain

White noise looks like the plot below on the left in the time domain. A key point is no matter how much you zoom into it, it will still generally look like this. Why? Because by definition every sample (let's assume a sampled system although this is not limited to that) is uncorrelated from the next in a white noise signal. The auto correlation is therefore an impulse.

If it happens to be a Gaussian White Noise process, then it will also have a magnitude distributed as shown in the histogram on the right, but these are two different characteristics, it can be White and not Gaussian, or Gaussian and not White.

Now for comparison, below is the same result for a "Colored Noise", and in this case, Gaussian distributed process.  Interestingly here the histogram of the distribution, given sufficient samples is identical, but the key point is as you zoom into the time domain plot on the left, we see the result of "memory" due to the bandwidth limiting that was done. Because of low pass filtering (in this case), the signal cannot transition from one sample to any other arbitrary location on the next sample; the location of the next sample depends on previous samples and is therefore constrained. (The more we low pass filter, the longer the effective memory depth). So there is more predictability on the next sample based on previous samples; samples are correlated to adjacent samples, and the auto-correlation function will no longer be an impulse but will be spread out with a width depending on how narrow the filtering is in frequency.

In the frequency domain, the power spectral density for the "Colored Noise" could look like the plot below; in particular this plot is the low pass filtered signal discussed above (but colored noise could be any selective band of frequencies, or even multiple bands). In contrast the white noise would have the same level of power spectral density for all frequencies. So in case this point wasn't obvious to all readers; white similar to white light that contains all frequencies vs colored with again the analogy to light in that individual colors are given by selective frequencies of light.


Is this an entirely stupid idea?

No, but you've just came to the conclusion that instead of sampling complex, with Nyquist rate being the bandwidth, you should do twice as many samples.
That simply means you're not doing IQ sampling, but low-IF or direct-RF sampling.

Mix the signal with a higher intermediate frequency and filter it in such a way that successive samples are equivalent to the I, Q, -I, -Q

so, multiply the real and imaginary parts with 1, -1, 1, -1 … each , then sample them.
1, -1, 1, -1 is is a cosine at exactly half the sampling rate. What you're mathematically doing is multiplying your IQ signal with a cosine + j· sine, each of half the sampling rate. $\cos ft + j\sin ft = e^{jft}$, i.e. that's a simple frequency shift.
That simply shifts up your signal from $f \in [-b/2; +b/2]$ to $f \in [0; b]$. I.e. you've just put your $b$-Bandwidth signal on an intermediate frequency of $b/2$. That's a low-IF heterodyne receiver, which existed before IQ was even discovered.
So, clearly, you can do that. I'd avoid shifting to only b/2, use something slightly higher if your bandwidth and sample rate allows, and forget all you've learned about complex sampling – we're back in real domain, and you simply have to sample with a rate that's higher than twice the bandwidth.

Also: As much hype there is around the RP2040: There's nothing new about having two cortex-m0 cores. It's not an especially fast microcontroller or anything, it really doesn't open any "new doors". The things that are nice about it are its flexible IO logic, but you need none of that.
In fact, these cores are ubiquitous, basically every semiconductor company sells microcontrollers based on that core. Although it's dual-core and clocked at a rather high rate, it's the last microcontroller that makes sense here: You need something that can process samples, not two high-clocked cortex-M0 cores that have no floating point and only rudimentary integer arithmetics. Get something with a Cortex-M4F instead. You might be getting a faster ADC, too.


The transfer function is $H(s) = \frac{16.94s + 579.5}{s^2 + 507.2s + 1224}$
This transfer function has 2 poles, one slow pole at -2.4248 and a fast pole at -504.7752. The function has a slowish zero at -34.2. Good  news, your poles and zero are all in the left-half plane. It is much easier to control a system with zeroes and poles in the left-half plane than in the right-half plane. Your slow pole at -2.4248 will limit your performance. Your slowish zero at -34.2 will probably limit your performance by creating some overshoot, or maybe not. Your fast pole will have a limited impact on the closed-loop performance, so you can simply ignore it for now.
You could design your controller to cancel the slow pole at -2.4248 and the slowish zero at -34.2. This would be a lead lag compensator looking like this
$$H_{LeadLag} = K\frac{s + 2.4248}{s+34.2}$$
Then you could set a PI in series with the lead-lag compensator and tune the P and I gain to get the closed-loop performance you want.
$$H_{c} = H_{LeadLag} \frac{K_ps + K_i}{s}$$
That being said, I'm not sure my solution is particularly robust. You should check the Gang of four stability functions https://www.cds.caltech.edu/~murray/courses/cds101/fa02/caltech/astrom-ch5.pdf.
Edit 2 :
Another solution, simply use a PI controller and set the ratio of the P and I gain to 2.4248 to cancel your slow pole. This controller will not cancel your slowish zero, but maybe you don't need to cancel it to get the performance you want.
$$H_{PI} = K\frac{s+2.4248}{s}$$


Why do we need to flip the kernel in 2D convolution?

At the risk of repeating what has been said elsewhere, this is how convolution in two dimensions is defined:

$$y(n_1,n_2) = \sum_{k}\sum_{l}x(k,l)h(n_1-k,n_2-l))$$

$x(n_1,n_2)$ is your data and $h(n_1,n_2)$ is the kernel. Note that in order to compute the double sum above, you need to flip $h(n_1,n_2)$ with respect to $x(n_1,n_2)$. Just take $n_1=0$ and $n_2=0$: in this case you need to multiply $x(k,l)$ with the flipped kernel $h(-k,-l)$. In order to compute the output for other values of $n_1$ and $n_2$ you need to move the flipped kernel across the input data.


Here is the entire simulation code.
First, to verify the spectrum response of ZOH, I replaced the DAC_output with an impulse signal. So I can directly check if the simulation confirms the theory or not. ( As ZOH is a kind of LTI filter, the spectrum of impulse response represents the frequency response of this system. You can try other input signals like the one you originally provided.)

Input test signal:

Output Spectrum: Note that theory spectrum entirely cover the simulated spectrum.

Ps. DAC is simulated via an up sampling and a convolution with hold filter, h = ones(1, sample_hold)

A very high pseudo analog sample frequency should be given because there is no real 'analog' in computer simulation.

You can alter different Tp or Ts (not lower than analog pseudo freq.). These codes should work.

There are some rooms for improvement about the ADC process you provided in question. You can refer the my blog.

clear all;
close all;

end_time = 2e-2;
x=0.0001:0.0001:end_time; %10KHz ADC sampling

ADC_response = sinc(x); %% Please check the frequency response!!

fb=1e3;  %1Khz frequency

f_noise=4e3;  %1Khz frequency

yb=sin(2*pi*fb*x); %1Khz frequency sine input signal

yf_noise=sin(2*pi*4*fb*x); %4Khz frequency sine input noise

y_in=yb+yf_noise;%total signal going into the ADC

ADC_out=conv(y_in,ADC_response);%signal after ADC

%% Started here.

% ADC_out = ADC_out(1:length(y_in)); % trim for same length

ADC_out = zeros(1,length(y_in));
ADC_out(length(y_in)/2) = 1; % impulse

%% pseudo analog frequency & ZOH simulation
% parameters
analog_sampling_rate = 1e8; % 100 MHz.
digital_sampling_rate = 1e4; % according to question.

%% Important Time Parameters
Ts = 1/digital_sampling_rate; % 1e-4 as given.
Tp = 2e-5; % ZOH time in analog unit (sec), tunnable.
hold_sample = round(Tp*analog_sampling_rate);
%  This indicates how many samples need to be hold in this pesudo analog simulation.

%% DAC Processing
t = 0:1/analog_sampling_rate:end_time-1/analog_sampling_rate; % pseudo analog time.

% upsampling
y_up = upsample(ADC_out,analog_sampling_rate/digital_sampling_rate);

% zero order hold filter
h = ones(1,hold_sample); % zero order hold

% filter your signal
y_zoh = filter(h,1,y_up);
y_zoh(1:floor(mean(grpdelay(h))))=[];
y_zoh = [y_zoh zeros(1,floor(mean(grpdelay(h))))];

%% Time domain plot
figure;
plot(x,ADC_out);
hold on;
plot(t,y_zoh);
legend('yb','y_zoh');

%% Freq response of discrete time signal and freq unit conversion
d_N = length(ADC_out);
d_F = fft(ADC_out,d_N);
d_freq_step = digital_sampling_rate/d_N;
d_freq = -digital_sampling_rate/2:d_freq_step:digital_sampling_rate/2-d_freq_step;
d_freq = d_freq(d_N/2+1:end)';

%% Freq response of psudo analog time signal and freq unit conversion
c_N = length(y_zoh);
c_F = fft(y_zoh,c_N);
c_freq_step = analog_sampling_rate/c_N;
c_freq = -analog_sampling_rate/2:c_freq_step:analog_sampling_rate/2-c_freq_step;
c_freq = c_freq(c_N/2+1:end)';

%% Theory transform of DAC Transfer function.

c_thy_F = abs(sinc(c_freq*Tp)*Tp/Ts); %% For verification.

%% Plot Freq domain comparison.

figure;
plot(d_freq,abs(d_F(1:d_N/2))/d_N);
xlabel('Freq. (Hz)');

figure;
plot(c_freq,abs(c_F(1:c_N/2))*d_N/c_N);
hold on;
plot(c_freq,c_thy_F);
xlabel('Freq. (Hz)');
legend('simulated','theory');


Determining sampling rate in OFDM demodulation

It is incorrect to say that the sample rate is "much less" than the bandwidth of the signal.  It is true to say that it is a "little less".  Even that statement is only true for part of the demodulator.

Demodulators almost always use a sample rate that is higher than the bandwidth for two reasons: to ensure that all of the signal is captured rather than aliased due to Nyquist frequency limitations, and to do some version of early/late gating when syncing with the symbols.  It is very common to sample at twice the symbol rate because that is very convenient for capturing the symbols and doing early/late gating.

OFDM is different because of the FFTs/IFFTs built into the signal.  There is no advantage for the demodulator to do the FFT at higher than the symbol rate, so they don't.  The roll-off in OFDM is very minimal (an example of typical OFDM roll-off is shown below), though, so the sample rate is not much less than the bandwidth.  In some schemes, like 802.11a, the outer frequency bins are not used which reduces the bandwidth of the signal relative to the sample rate, meaning that the sample rate will be higher than the bandwidth even when taking into consideration the signal roll-off.

Also, before the demodulator does the FFT it has to eliminate as much carrier offset as it can.  This is typically done at a sample rate that is higher than the bandwidth for the same Nyquist reasons as a typical demodulator.

EDIT: I see that I forgot to answer some of your questions.  The sample rate for OFDM is exactly equal to the sample rate that the OFDM transmitter used when it inverse FFT'ed the data.  No, the sample rate has nothing to do with the modulation types used inside the OFDM symbols (e.g. BPSK or QPSK).  No, the sample rate has nothing to do with the FFT length.


Is spectral density sometimes normalized by sampling rate rather than bin size?

The use of $rad/\sqrt{\text{Hz}}$ suggests that this is phase noise specifically (a spectral density due to phase fluctuations), and typically in my use this has been described as a power spectral density (units of $rad^2/\text{Hz}$), so this is just the square root of that quantity.
The reason the DFT (of which the FFT computes) is divided by $N$ is to normalize the FFT to be the same units of the time domain signal, specifically using the following normalized form of the DFT:
$$X_1(k) = \frac{1}{N}\sum_{n=0}^{N-1}x[n]W_N^{nk}$$
Versus the typically version that is not normalized which the FFT returns:
$$X(k) = \sum_{n=0}^{N-1}x[n]W_N^{nk}$$
With such a normalization, the magnitude of $x[n]$ at any specific frequency will match the magnitude of $X(k)$ for that frequency. For example, if we had a time domain waveform of a sinusoidal phase error versus time given as:
$$\phi[n] = A\cos(\omega n) = \frac{A}{2}e^{j\omega n} + \frac{A}{2}e^{-j\omega n} \space \text{rad}$$
Then assuming $\pm\omega$ were exactly on a bin centers (for the DFT due to its circular nature $-\omega = N-\omega$), the resulting two bins in $X_1(k)$ would have a magnitude of $\frac{A}{2}$, matching the magnitudes of the time domain waveform.
As a power spectral density (meaning we are interested in the power over a given frequency range) the normalized power of each frequency index in the DFT (aka bin) is then:
$$|X_1(k)|^2 = \frac{|X(k)|^2}{N^2} \space \frac{\text{rad}^2}{\text{bin}}$$
(Where the units of $\text{rad}^2$ for the power quantity $|X_1(k)|^2$ only make sense if x[n] was the phase noise in units of radians).
$\frac{\text{rad}^2}{\text{bin}}$ is a power quantity per bin. To make this the recognized form of the power spectral density in power/Hz we recognize that $Nd = f_s$ where $N$ is the number of samples in the DFT, $f_s$ is the sampling rate, and $d$ is the spacing of each frequency index (bin as the OP used) in Hz resulting in the spectral width of each bin in Hz:
$$d = \frac{f_s}{N} \space \frac{\text{Hz}}{\text{bin}}$$
Thus
$$ \frac{|X(k)|^2}{N^2} \frac{\text{rad}^2}{\text{bin}} \times d^{-1} \frac{\text{bin}}{\text{Hz}} = \frac{|X(k)|^2}{N^2}\frac{N}{f_s} \frac{\text{rad}^2}{\text{Hz}} = \frac{|X(k)|^2}{N f_s} \frac{\text{rad}^2}{\text{Hz}}$$
This result would specifically be what we typically notate as $\scr{L}_{\phi}(f)$ as the two-sided power spectral density due to phase fluctuations (since the DFT contains both sides of the spectrum, in contrast to the one-sided PSD which is $S_\phi(f) = 2\scr{L}_{\phi}(f)$.).
Note we say "due to phase fluctuations" since the units here were phase. It is also interesting how the phase unit in radians when squared is the power unit relative to the carrier (often expressed as dBc/Hz). This is clear for small angles given the small angle approximation $sin(\theta) \approx \theta$, or geometrically the quadrature component being the noise as phase noise relative to the in-phase component being the carrier that has been rotated due to that phase, such that the ratio of the two is the phase unit in radians, for small angles!) This is why when phase noise is dominant this computation will match the actual power measurement we see under test with a spectrum analyzer.

Further update:
The OP clarified in his comments that his question is specific to the peak at 30 Hz offset as shown in this plot:

It isn't specified but assuming this is a two-sided spectral density, the peak of a single tone would have a total power independent of density, so we would typically report its result as $\text{rad}^2$ and not $\text{rad}/\text{Hz}$ (or the magnitude quantity as the square root $\text{rad}$ as used in this plot, meaning this plot is $\sqrt{\scr{L}_{\phi}(f)}$). The paper also incorporates a moving average of 5 and suggests in a foot note that the peak would be $\approx 1.6 \text{nrad}/\sqrt{1\text{kHz}}/5$, and the plot was scaled (moved up or down) such that the level of the tone landed on this expectation.
I suggest that the peak would be at either $\approx 1.6 \text{nrad}/20$ or $\approx 1.6 \text{nrad} \sqrt{2}/20$ depending on if the spectrum is intended to be double-sided or single-sided which should be specified. The sampling rate does not change the value of the tone on the spectral density when the units are already in nrad, so there should also be no $\sqrt{1\text{kHz}}$ in that answer - The sine wave theoretically occupies zero bandwidth, or for practical reasons we can assume we integrated that power over a small bandwidth in order to measure the peak we see. Either way the density becomes a single figure for the tone independent of bandwidth. Any windowing applied in the time domain prior to the FFT (other than the rectangular window) will also shift the value of the tone differently from the values for the noise. Further details below.
To confirm that assumption, here is my prediction of where such a tone would be:
The 1.6 nrad oscillation is specified as the peak to peak value and thus is of the form:
$$\phi(t) = \frac{1.6}{2} \cos(2\pi f t) \space\space \text{nrad}$$
with $f=30e3$
If the spectrum is two-sided (as $\sqrt{\scr{L}_\phi(f)}$ rather than one-sided as $\sqrt{S_{\phi}(f)}$), then the spectrum is only showing the upper half of this two-sided spectrum, with both sides given by:
$$\phi(t) = \frac{1.6}{2} \cos(2\pi f t) = \frac{1.6}{4}e^{j 2\pi f t} + \frac{1.6}{4}e^{-j 2\pi f t} \space\space\text{nrad}$$
Thus prior to the effect of the moving average filter (MAF), I would predict the tone shown on a double-sided spectrum to be at:
$$\frac{(1.6e-9)}{4} = (4e-10) \space \text{rad}$$
Notice the units are $\text{rad}$ and not $\text{rad}/\sqrt{\text{Hz}}$ as the standard deviation of the tone itself is not a density spread across frequency, unlike that of the noise.
I assume the moving average filter that is mentioned was done on the frequency domain samples. If in the time domain there would be an additional loss of 0.963 but I don't see evidence of such a moving average response in the plot, in which case with a moving average of frequency samples, the tone is reduced by a factor of 5 as the author had done, resulting  in $(4e-10)/5 = (8e-11)$.
If the plot was supposed to be a single-sided spectrum $\sqrt{S_{\phi}(f)}$, then the result would be $\sqrt{2}$ larger or $1.13e-10$, which is consistent with the standard deviation of $\phi(t)$ reduced by the MAF.
Neither of these results match the plot, but this is where I would expect a 30 Hz tone after a moving average of 5 samples when sampled at 1 KHz if the units of the spectral density are $\text{nrad}/\sqrt{\text{Hz}}$, for either case of a single-sided or double-sided spectral density. Also note that my computation was independent of the bin size or number of samples since as the author of the paper was intending to do (and perhaps did if I made an error in my prediction) was to predict the expected value of that tone and then scale the plot accordingly. My earlier answer shows how I would scale the result from the DFT directly in which case the bin size and number of samples would be involved.
As a further note since these spectrums are being derived from FFT's and since the OP is interested ultimately in assessing noise: We must also be careful to account for the equivalent noise bandwidth due the effect of windowing especially if we are normalizing the plot based on the power of a tone. (and other effects such as scalloping loss etc which have been minimized by choosing a tone at or near a bin center as was done). Any windowing done on the time domain signal other than the rectangular window will widen the bandwidth of each bin beyond the single bin as given by the rectangular window, which means that the noise measured will be larger than the actual noise! Further the window has a loss reducing the signal from the tone and the noise, but because of the effectively wider noise bandwidth of each bin the noise will go down less than the tone (the tone only occupies one bin)! The effect of the moving average in frequency on SNR is also affected by the window since the adjacent noise bins are no longer uncorrelated. I detail this further in this post: Find the Equivalent Noise Bandwidth


Using white noise to test filter freq. response?

Yes you can also identify a system (more spcifically an LTI system) by inputing white noise to it and measuring the corresponding output. Yet this will be a statistical approach compared to the sine-wave measurement which is deterministic.

First, you have to be able to generate a WSS (wide-sense stationary) white noise of power $\sigma_x^2$ as the input $x[n]$ for testing.

Then you have to measure the output signal $y[n]$ and estimate it's power spectral density (PSD), assuming it's also WSS (which is indeed). So you will use some statistical approaches to obtain those quantities.

Then the relation between the input PSD of the white noise of power $\sigma_x^2$ and output PSD of the (colored) noise through the LTI is:

$$ S_{yy}(e^{j\omega}) = |H(e^{j\omega})|^2 S_{xx}(e^{j\omega}) $$

where $S_{xx}(e^{j\omega})$, and $S_{yy}(e^{j\omega})$ are the discrete-time input and output PSDs respectively. Since $S_{xx}(e^{j\omega}) = \sigma_x^2$ the output PSD will be:
$$ S_{yy}(e^{j\omega}) = |H(e^{j\omega})|^2 \sigma_x^2 $$

However, what is worth noting here is that,  this method yields not the frequency response $H(e^{j\omega})$, but its magnitude square. So the phase information is lost.

A better approach, which yields the frquency response $H(e^{j\omega})$ directly is to use cross-correlation or cross-power spectral densities instead.
The relation between input noise $x[n]$, output noise $y[n]$ can be expressed as

$$ S_{xy}(\omega) = H(\omega) S_{xx}(\omega) $$

Which shows a relation between input PSD and cross-PSD, and reveals $H(\omega)$ directly when input PSD is known as in the case with white-noise input PSD of $S_{xx}(e^{j\omega}) = \sigma_x^2$:

$$ S_{xy}(\omega) = \sigma_x^2 H(\omega). $$

Here the cross power spectral density, $S_{xy}(\omega)$ between the input $x[n]$ and the output $y[n]$ can be obtained through the DTFT of the cross-correlation between $x[n]$ and $y[n]$ from

$$\phi_{xy}[m] = \mathcal{E}\{x[n]y[n+m]^*\}$$

Which is indeed equivalent to the impulse response of the LTI system to be identified :

$$\phi_{xy}[m] = \sigma_x^2 h[m]$$


How to implement a $j$-level $M$-band wavelet transform of an image?

[EDIT: some code made available]  The  notion of $M$-band wavelet transform, $M\ge2$, generalizes the standard $2$-band wavelet. The theory is provides, for instance, in Theory Of Regular M-Band Wavelet Bases, P. Steffen et al., 1993. Remember that with $2$-band wavelets, one cannot obtain wavelets that are real-valued, symmetric, orthogonal and with finite support, except for the Haar wavelet (not very regular). For $M>2$, such designs can be obtained with sufficient regularity, and $M$-band wavelets are sometimes believed to be more efficient for image analysis.

Alas, regular $M$-band filterbank are generally harder to design, and $M$-band wavelets have been used more rarely. Resultingly, they are less implemented in standard software.

On $M$-band implementations: I am not aware of build-in Matlab functions. However, there exist $1$D versions, that you can built upon for your applications:

LTFAT (Large Time-Frequency Analysis Toolbox) : WFILT_CMBAND - Generates $$M$-Band cosine modulated wavelet filters, from On cosine-modulated wavelet orthonormal bases, 1995
M-band dual-tree wavelet: Shannon & Meyer wavelets, and Modulated Lapped Transforms for specific $M$s, and a $4$-band filter from Design of Efficient $M$-Band Coders with Linear-Phase and Perfect-Reconstruction Properties, 1995, Alkin and Caglar. They  come with their Hiilbert transform, and are detailed in Image Analysis Using a Dual-Tree M-Band Wavelet Transform, 2006; Noise Covariance Properties in Dual-Tree Wavelet Decompositions, 2007; A Nonlinear Stein Based Estimator for Multichannel Image Denoising, 2008. Below, 2D oriented $4$-band wavelets from Alkin and Caglar:

A $2$-D version is now embedded in the Matlab toolbox:  $M$-band $2$D dual-tree (Hilbert) wavelet multicomponent image denoising. Other sources to grab coefficients from:

-An algebraic construction of orthonormal $M$-band wavelets with perfect reconstruction

And one day, I shall rewrite my own codes to make them much cleaner.


Interpretation of Clarke's Doppler power spectral density

A simple, "non-technical" way of thinking of it is the fact that the Doppler frequency is proportional to $\cos\theta$. The amplitudes of cosine, however, are not uniformly distributed, but are heavily weighted towards $\pm 1$.

Example plot to demonstrate, using Python/Pylab code:

theta = linspace(0, 2*pi, 1001)
x = cos(theta)
hist(x)

More rigor can be seen by noting that
\begin{align}
f &= f_d \cos\theta\\
\theta &= \cos^{-1}\left(\frac{f}{f_d}\right)
\end{align}
and the power received at any angle is proportional to a small angle increment $d\theta$:

$$
P(\theta) \propto d\theta = \frac{-1}{f_d\sqrt{1-\left(\frac{f}{f_d}\right)^2}} df
$$

And the total power can be determined by integrating the above quantity, which is identically what defines a power spectral density.


How to get the scale bar of an FFT for a 2d image

If you are having trouble visualizing the reciprocal space, perhaps an analogy to 1D temporal FFTs will help (personally I'm more familar with these, so an analogy helps me anyway)

In a 1D temporal FFT, the input is N samples, at a rate of Fs samples/second, for a total time of N/Fs seconds.  The output is N samples, distributed between 0 and Fs (or between -Fs/2 and +Fs/2 depending on your point of view) (units 1/second).  So each output sample has width Fs/N (again units 1/second).

In a 1D spatial FFT, all we do is replace seconds by meters and samples by pixels. The input is N pixels, and a rate of Fx pixels per meter, for a total width of N/Fx meters.  The output is N pixels, distributed between 0 and Fx (or between -Fx/2 to Fx/2) (units of 1 / meter).  So each pixel has width Fx/N (agains units 1/meter).

So an output scale bar of M pixels has a width of M*Fx/N in units of 1/meters.

To go from 1D spatial FFT to 2D spatial FFT, just need a sampling rate in each direction.  Presumably you have the same sample rate in each directin.


Sum of correlated exponential random variables

The pdf $f_Z(z)$ of the sum $Z=X+Y$ of any two jointly continuous random variables $X$ and $Y$ with joint pdf $f_{X,Y}(x,y)$ is as follows:
$$\text{For all } z, -\infty < z < \infty, ~~ f_Z(z) = \int_{-\infty}^\infty f_{X,Y}(x,z-x) \, \mathrm dx.\tag{1}$$

For the special case when $X$ and $Y$ are nonnegative random variables (including as a special case, exponential random variables) and so take on nonnegative values only, $f_{X,Y}(x,y)$ has value $0$ if at least one of $x$ and $y$ is smaller than $0$. Hence, in this case, the integrand $f_{X,Y}(x,z-x)$ in $(1)$ has value $0$ if $x < 0$ or if $z < x$. Consequently, if $z$ is a negative number, then the integrand in $(1)$ is always $0$ regardless of the value of $x$ and therefore so is the integral. All of which is just a long-winded way of saying that $f_Z(z)$ has value $0$ when $z<0$, that is, $Z$ takes on nonegative values only, which any idiot could have deduced from the fact that $Z=X+Y$ and both $X$ and $Y$ are nonnegative. But the approach is useful even for $z>0$ since now we have that the integrand in $(1)$ is zero when $x<0$ or when $x >z$ and so for nonnegative $X$ and $Y$, we can simplify $(1)$ to
$$f_Z(z) = \begin{cases}\displaystyle\int_0^z f_{X,Y}(x,z-x) \, \mathrm dx, & z \geq 0,\\\quad\\
0, & z < 0\end{cases} \tag{2}$$

  No further simplification of $(2)$ is possible in general.

For the special case when $X$ and $Y$ are independent random variables, $f_{X,Y}(x,y)$ factors into $f_X(x)f_Y(y)$ and so $(1)$ becomes the familiar convolution integral and $(2)$ the somewhat-less-familiar convolution integral for causal signals. But no such simplification is possible for nonindependent random variables $X$ and $Y$; we need the joint pdf to calculate $f_Z(z)$ and just knowing that $X$ and $Y$ are correlated random variables (whether exponential or Gaussian or whatever) is not enough.


General questions regarding OFDM modulation

The output time domain samples are also of same length $N$. Your figure shows a part of a OFDM receiver. At the OFDM transmitter, IFFT is used.
One OFDM symbol duration is $NT$ and obviously one sample of the OFDM symbol has a duration $T$.


Recommended Resources / Literature Search Terms for a Solutions to a Specific Kind of Multi Harmonic Signal Structure

If I understand this problem correctly you have access to 2 signals:

Noise Signal - $ w \left[ n \right] $. It is composed of a linear combination of harmonic signals. Something like $ w \left[ n \right] = \sum_{i}^{m} {a}_{i} \sin \left[ 2 \pi \frac{ {f}_{i} }{ {f}_{s} } n + {\phi}_{i} \right] $.
Input Signal - $ y \left[ n \right] $ which is composed of the signal of interest $ x \left[ n \right] $ and the same harmonics as in $ w \left[ n \right] $ just with different weights. Something like $ y \left[ n \right] = x \left[ n \right] + \sum_{i}^{m} {b}_{i} \sin \left[ 2 \pi \frac{ {f}_{i} }{ {f}_{s} } n + {\psi}_{i} \right] $.

The objective is to estimate $ x \left[ n \right] $ given those signals.
Since we can model $ y \left[ n \right] = x \left[ n \right] + \left( h \ast w \right) \left[ n \right] $ we can utilize the adaptive filter framework and solve this in highly efficient manner.
Remark
In real world each observed signal would also have some kind of measurement noise we can model as AWGN. But let's leave it a side at the moment.


My naive approach would be to just place a grid over the IQ map, and find which grid-square my symbol lands in.

And that's (more or less) exactly how it's done: it's called a hard decision, and your  grid is what we'd call the decision boundaries. It's a result of applying the Maximum Likelihood criterion to the decision process.

Of course, in that decision, a sample that hits a constellation point (nearly) exactly is not distinguishable from a sample that was somewhere "in the middle" between two constellation points (and hence close to one of the decision boundaries).

You can't tell "absolutely certain" from "I'm about 51% sure" decisions that way! So, this hard deciding throws away valuable information.

If you find a way to encode the decision information and the certainty of that decsion into numbers, then you get a *soft-output decider". The most common way of finding a single number that describes that is through likelihood ratios, i.e. divide how likely one symbol is by the likelihood of a different symbol. Because ratios are unelegant to deal with, apply a logarithm to it, and get a LLR.

Now, you want bits, not something like "I'm almost certain this was a 0 or 1", but don't forget:

In modern systems, you use channel coding. If you have an algorithm that doesn't only consider single symbols, but sequences of symbols, then you can use these soft-outputs to infer the likeliest sequence that was sent more reliably.

As said in my comment to your question, you really shouldn't start learning about QAM mapping by reading about soft-output decision on QAM. Start with soft-output BPSK, and then learn about the hard- and soft-output Viterbi sequence estimator. (For Viterbi, you'll need working knowledge of simple convolutional channel codes, but that should be comparatively simple when considering you're doing LLRs on QAM...)


Is there a way to use a Kalman Filter to improve these measurements if I have no physical rule...

No.
Kalman filters add value by using a model of a process to refine an estimate of the states of that process.
If you don't have a model, and you don't have enough data to develop a model, then you can't get there from here.


Is sampling at double the desired reproduction frequency accurate?

There is not only one sampling rate that theoretically allows exact reproduction. The sampling theorem states that it is sufficient for perfect reconstruction if the sampling frequency $f_s$ is greater than twice the highest signal frequency:

$$f_s>f_{max}\tag{1}$$

Without further restrictions, Eq. $(1)$ is also necessary for baseband signals occupying the whole frequency range $[0,f_{max}]$. Some people like to discuss whether in $(1)$ there should be a "greater or equal" sign or just a "greater" sign. In practice, this question is totally irrelevant because you always have to sample at a frequency twice the maximum frequency plus some extra tolerance to allow for non-ideal filters.

As you've noticed, sampling a sinusoid at exactly twice its frequency will generally not allow you to reconstruct it. As soon as you sample at a rate slightly larger than twice its frequency you will be able to reconstruct it perfectly.

Note that there is band pass sampling for band pass signals where the sampling rate can be considerably smaller than twice the highest signal frequency for perfect reconstruction to be possible. However, it's generally not true that sampling at twice the signal's bandwidth is sufficient, even thought that may be the case for specific signals.

These are only the basics of (uniform) sampling of 1-dimensional signals. There is much more to be learned. Browse the "sampling" tag on this site, and read the overview article Sampling - 50 Years After Shannon by M. Unser.


Sample Rate: How to determine sample rate

This question isn't particularly on-topic, but the answer is pretty simple. Look at the second paragraph of the General Description on the first page of the ADC datasheet:

  The conversion process and data acquisition are controlled using
  CS and the serial clock signal, allowing the device to easily interface
  with microprocessors or DSPs. The input signal is sampled on the falling edge of CS and conversion is initiated at this point. There are no pipeline delays associated with the part.

CS in this case refers to the chip select pin on the ADC's SPI-compatible 4-wire serial bus. That signal is asserted each time the host processor initiates a transaction to the converter. The mechanics of an SPI read go like this:

The SPI master (i.e. the host processor) asserts the CS signal.
If the SPI master has any data to write to the device, it presents it in serial fashion on the MOSI (master out, slave in) signal. From the ADC's perspective, this would be the DIN pin.
The SPI master begins to apply a clock signal to the ADC on the SCLK pin. This is used to simultaneously shift data in and out of the ADC serially.
For this ADC device, the transaction continues for 16 bits. After that, the master has received 16 serial bits from the ADC, which contain the conversion result.
The master deasserts the CS pin to end the transaction.

In the above timeline, the speed of the conversion is therefore limited by two factors:

The SPI clock rate (as it determines how long it takes to extract 16 bits).
The maximum conversion rate that the ADC supports.

#2 in the above list is very clearly stated to be 1 MHz. #1 is dictated by the capabilities of your host processor, and can be found in the datasheet. Specifically, look at the bottom of page 166:

  Bit 0 – SPI2X: Double SPI Speed Bit

  When this bit is written logic one the SPI speed (SCK Frequency) will be doubled when the SPI is in Master mode (see Table 18-5
  ). This means that the minimum SCK period will be two CPU clock periods.

It is stated elsewhere in the datasheet that a "CPU clock period" as referenced here is equivalent to $\frac{1}{f_{OSC}}$. Therefore, if your processor has a CPU clock frequency of 8 MHz, its maximum SPI clock rate is:

$$
f_{max,SPI} = \frac{f_{OSC}}{2} = 4\text{ MHz}
$$

Remember though that this is a serial bus with 16 bits per conversion transaction. Therefore, the maximum ADC sample rate is:

$$
f_{max,ADC} = \frac{f_{max,SPI}}{16} = 0.25\text{ MHz}
$$

Thus, the processor is the limiting factor here; it can only support a maximum of a 250 kHz sample rate. Furthermore, this is the maximum throughput that is available to the entire ADC device, so if you're using more than one channel on the ADC, this 250 kHz would have to be split between the channels you're using.


Relative performance of RLS and LMS filters

Using less data is helpful when, as you said, the parameters are time varying, and in particular when they change a lot.

The key difference is that LMS is a Markov process.  It has its current state, but other than that it does not remember data from the past.  For time-varying signals this is a feature because past data will give you erroneous information about the current parameters.

The RLS algorithm uses all of the information, past and present, but that can be a problem if the past data is misleading for the current parameters.

If you are looking for a quantitative rule for when to use one or the other, I don't have one.  RLS is more computationally intensive than LMS, so if LMS is good enough then that is the safe one to go with.  RLS converges faster, but is more computationally intensive and has the time-varying weakness, so I would only use it if the parameters don't vary much and you really needed the fast convergence.


How Random Forests can be used for regression?

The key is in choosing the split function of the decision tree. In classification, the data is split according to a certain homogeneity metric. On the other hand, for a regression tree, a regression model is fit to the dependent variable, for each of the independent variables. Then, the data is split at multiple points and the SSD value is evaluated. The best split is chosen for the dimension (variable) giving the minimum error. The node is then split at this point. The splitting is again recursive.

Finally, each leaf stores a distribution for the continuous output variable/s. In the testing phase, the responses are averaged.


Efficient real-time pitch detection algorithm

YINS fft is the best method that worked for me in the past. I have my open source implementation of it written in awesome D language: link . I compared the results with Praat they match.


How to Deduce a Linear System's Impulse Response from a Set of Input and Output Signals?

There seems to be an image which is gone now and hence I might be missing something.

In order to declare if the system is invariant you should see if a delay of the input yields only a delay in the output.
In you case, does the input $x_1[n-m]$ yields $y_1[n-m]$ and etc...

As written by @DilipSarwate, By looking at the response for $ {x}_{2} - {x}_{1} $ one could see the system is neither causal nor time invariant:

consider that linearity requires that the response to $ {x}_{2} − {x}_{1} $ be
$ T \left( {x}_{2} − {x}_{1} \right) = T \left( {x}_{2} \right) − T \left( {x}_{2} \right) = {y}_{2} − {y}_{1} $.
Now, $ {x}_{2} − {x}_{1} $ is a pulse of amplitude $ 2 $ at $ n = 2 $ and the response to it is nonzero at $ n = 0, 1, 2 $ and so the system is not causal, and since $ {x}_{2} − {x}_{1} $ is just $ {x}_{2} $ delayed by 2 units while the response is not $ {y}_{2} $ delayed by 2 time units, the system is time varying, not time invariant.

If the input signals are band limited and their bandwidth is less than your system you won't be able to restore the impulse response.
You will be able only to get the response in the frequencies the input has energy in.
This could be done by frequency analysis of the input and the output.
If your system is indeed LTI the connection between input and output is given by convolution with the impulse response.
Convolution is multiplication in the frequency domain, hence you could easily get the impulse response (Again, only at frequencies the input has energy in).

Yet since this system isn't time invariant it doesn't have a single impulse response. It has a response to an impulse at $ n = 0 $ but it might have a different one for an impulse at $ n = 1 $.
Update
This is a nice case to show the commutative property of convolution.
Since $ y \left[ n \right] = \left( h \ast x \right) \left[ n \right] = \left( x \ast h \right) \left[ n \right] $ you can switch the roles between the "System" and the "Input Signal".
So this question is equivalent to Given known output and system response, how could one estimate the input of the system?
Well, this is a known problem called Deconvolution and in the case of no noise it is simple to solve.
As written above, one way to do so it write the problem in matrix form.


QAM constellation slicer with unknown carrier phase

What you need is carrier phase synchronization. This is a complicated topic with many different approaches. The approach that you'll choose could depend on things like:

Data-aided versus blind: Does the underlying sequence contain any known data (e.g. a training or sync sequence of some kind) that you can use to divine the phase offset? Or, do you have to synchronize with no knowledge of the modulating symbols?

Blind approaches are more general, but you can get better phase estimation performance if you're aided by the data. Also, blind approaches often introduce ambiguities into the recovered phase (i.e. there are multiple phase offset solutions that equally satisfy the blind optimization criterion).

I haven't come across any blind solutions that work equally well across all modulation types; for instance, most blind phase estimators are best applied to PSK signals. They can be made to work, suboptimally, on QAM. Various blind carrier phase estimators for QAM are documented in the literature, but I don't have any great recommendations for QAM-specific estimators. The most commonly-used blind estimator for PSK-type signals is the power-law detector. Here is an example of a paper that talks about its application to QAM.
Coherent versus differential: As you noted, one way to get around having to synchronize to the carrier phase directly is to use differential modulation. In that case, information is carried by the difference in phase between consecutive symbols. Since the carrier phase is likely to be approximately constant across a two-symbol time period, the carrier component cancels out. This makes synchronization easier, but there is a symbol error rate performance loss on the order of 1-2 dB for differentially-coherent modulations versus fully-coherent operation.
Feedforward versus feedback: Are you processing a continuous, indefinitely-long stream of symbols, or do you have a finitely-sized batch? A feedback approach like a phase-locked loop might be appropriate for the former, while feedforward techniques that estimate the bulk phase offset of a block of symbols at a time are best for the latter (feedback techniques have some acquisition period during which you won't get good output; if you only have short blocks of data at a time, this can be a problem).

If you're looking for a book reference, my go-to is Synchronization Techniques for Digital Receivers by Mengali. It's expensive and hard to find a copy, but I find it very thorough.


What is the Reference in Control Theory?

Imagine that you're heating (or cooling) a home with a modern furnace (or air conditioner).

the reference or set point is the temperature that you set your thermostat to be.  the feedback signal is the actual temperature that is measured with some kinda thermometer.

the actual value that you are trying to control, whether it's temperature or the position of a robot arm or the position of a pointer in an asynchronous sample rate converter, that actual value is compared to your reference value.  the difference between those two signals is what drives your "controller".

in the case of heating your home, if the actual temperature is below the reference (what you want your temperature to be), the furnace heat is increased.  if the actual temperature is above the reference, the furnace heat is decreased.


Try NCH Software

From Wikipedia:

  NCH Software is a software development company founded in 1993 in
  Canberra, Australia. The Colorado office was started in April 2008 due
  to the large U.S. customer base. NCH Software primarily sells to
  individuals via their website.

And further

  Browser Toolbar

  During the installation process NCH Software also by default (but you
  can uncheck the option during install) installs Google Chrome. This
  can be unselected during install or removed after install. In older
  versions of their software a browser toolbar powered by Conduit,[11]
  was installed as well as making changes to the default homepage and
  search engine.


Oversampling in OFDM is usually implemented by modulating some subcarriers at the spectrum margins with zero. This zero-padding in frequency domain corresponds to oversampling in time domain. If only 50% percent of the subcarriers are used, then this corresponds to an oversampling factor of 2. Specifically, if the total number of subcarriers is $N$ and if $M$ of these subcarriers are set to zero, then the oversampling factor is $N/(N-M)$. Note, that the $M$ zero subcarriers have to be arranged in a row, cyclic shifting allowed. Normally, they're arranged symmetrically at the lowest and highest frequencies.

An OFDM transmitter might require a low pass filter in order to keep out-of-band radiation low. The cut-off frequency is usually half of the sampling frequency ($f_\mathrm{s}$).

An OFDM receiver requires a low pass filter before A/D conversion (sampling) with cut-off frequency $f_\mathrm{s}/2$. Otherwise aliasing can occur.


Orthogonal Codes for Band Limited Channel

Let $p_k(t) = p(t-kT)$ an orthonormal signal set for integer $k$ and $T>0$. In other words, we require that:

$$
\int_{-\infty}^\infty p(t-\alpha T) p(t - \beta T) dt =
\begin{cases}0, \text{  if $\alpha \neq \beta$}\\1, \text{  if $\alpha = \beta$}\end{cases} $$
The most common example of an orthonormal set is the square-root raised cosine pulse (SRRC). After matched filtering, SRRC pulses become raised-cosine pulses, which have zero ISI.

Let $s_1(t) = \sum_m a_m p_m(t)$ and $s_2(t) = \sum_m b_m p_m(t)$, where $a_m, b_m \in \mathbb{R}$ (the result below can be easily extended to the complex case).

The dot product $s_1(t) \cdot s_2(t)$ is

$$
\begin{eqnarray*}
\int_{-\infty}^\infty s_1(t) s_2(t) dt &=& \int_{-\infty}^\infty \left( \sum_m a_m p_m(t) \right) \left( \sum_m b_m p_m(t) \right) dt \\
&=& \int_{-\infty}^\infty \sum_m a_m b_m p^2_m(t) dt \\
&=& \sum_m a_m b_m.
\end{eqnarray*}
$$
In the second step, I used the fact that $\int p_m(t) p_n(t) dt = 0$ if $m \neq n$. In the third step, I used the fact that $p(t)$ has energy equal to one.

Then, the dot product is zero only when $\sum_m a_m b_m = 0$. Note that the assumptions stated above about $p(t)$ are crucial. If you use a pulse that does not meet these conditions (for example, a raised-cosine pulse instead of square-root RC), then the dot product $s_1(t) \cdot s_2(t)$ will not be zero even if the sequenes $a_m,b_m$ are orthogonal.


[EDIT] In 1991, Nasir Ahmed wrote: "How I Came Up with the Discrete Cosine Transform". Interesting to read, on how he was inspired by  Chebyshev polynomials, and on how he didn't get funding, for a tool at the heart of JPEG and MP3.
Natural images are not very stationary, but locally, their covariance is often modeled by a first- or second-order process. Being able to "concentrate" those matrices is one key to efficient compression (not the only one, however).
In different words: one aims at projecting, or representing, the pixels in another system of vectors (a basis, a frame), such that the information content of the picture is better represented, or sparser. Classically, one preserves energy (hence orthogonal or close to orthogonal transformations. Two main options:

compute or learn the best adaptive representation from each image. PRO: more efficient in compaction; CON: computationally expensive and storage lost because should be sent for decoding
use a known and fixed transform (or a set of known transforms). PRO: fast; CON: less optimal.

The DCT is option 2, with the performance of option 1: a fixed transform with close to adaptive performance, for most standard images. It worked, people still use JPEG and MP3, at the heart of which one finds the DCT.
The goal is thus to find the eigenvectors of such (covariance matrices, with less computational cost than matrix diagonalisation. The original DCT (DCT-II) from Discrete Cosine Transform (1974) aimed at finding basis vectors that could approximately diagonalize the first-order Markov Toeplitz  matrix:
$$ \left[\begin{array}{c}
1&\rho& \cdots &\rho^{N-1}\\
\rho& 1&\cdots &\rho^{N-2}\\
\vdots& \vdots&\ddots &\vdots\\
\rho^{N-1}& \rho^{N-2}&\ddots &1\end{array}\right] $$

The original motivation for defining the DCT was that its basis set
provided a good approximation to the eigenvectors of the class of
Toeplitz matrices that constitutes the autocovariance matrix of a
first-order stationary Markov process, with the result that it had a
better performance than the discrete Fourier transform (DFT) and some
other transforms [l], 3, 4 with respect to such kinds of
processes. In fact, as shown in 4, the DCT is asymptotically
equivalent to the Karhunen-Loeve transform (KLT) of a first-order
stationary Markov process as $\rho$ tends to 1, where $\rho$ is the
correlation coefficient.

Looking at the basis functions for $\rho=0.9$

the match is pretty good, and it turns out that the processing performance is close as well to that of the Karhunen-Loève transform.
In their 1995 paper Diagonalizing Properties of the Discrete Cosine Transforms, the authors provide an extension of the above property.

In this paper, we obtain the eight types of DCT’s as the complete
orthonormal set of eigenvectors generated by a general form of
matrices in the same way as the discrete Fourier transform (DFT) can
be obtained as the eigenvectors of an arbitrary circulant matrix.
These matrices can be decomposed as the sum of a symmetric Toeplitz
matrix plus a Hankel or close to Hankel matrix scaled by some constant
factors.

In 1991, Nasir Ahmed wrote: "How I Came Up with the Discrete Cosine Transform". Interesting to read, some excerpts:

What intrigued me was that the KLT was indeed the optimal transform on
the basis of the meansquare- error criterion and the first-order
Markov process model, and yet there was no efficient algorithm
available to compute it. As such, the focus of my research was to
determine whether it would be possible to come up with a good
approximation to the KLT that could be computed efficiently.

to study a “cosine transform” using Chebyshev polynomials

The motivation for looking into such “cosine functions”
was that they closely resembled KLT basis
functions for a range of values of the correlation coefficient
$\rho$ (in the covariance matrix).


Warning: this is still an active research topic so there is no "canned" solution.

As you have stated, middle channel suppression is not a very robust method. It can help in getting a very rough "Karaoke" version but you can't infer anything reliable from the resulting audio.

A simple approach would consist in extracting a bunch of audio features on each signal frame, and use machine learning techniques (for example support vector machines) to discriminate the "voice present"/"no voice" classes. This works relatively well (reported recognition rates in the literature are in the 70-80% range), but the downside is that you have to manually annotate material to get training data. You might try to use karaoke files (with precise timing of the lyrics) and align them with the audio to get a kind of automatic transcription, but remember that your system will only be as good as your training data. Your training database will have to cover all the possible genres of music on which your classification system will be unleashed. See Rocamora's Master's thesis or Ramona et al's paper. There are dozens of other papers presenting variations of this scheme - with different feature sets, post-processing/temporal-smoothing operations, etc. Follow the references in the papers.

If you can't afford to spend countless hours in Transcriber or Wavesurfer annotating data, there are still unsupervised methods that would reasonably work. I recommend you to look into lead melody separation methods such as:

"REpeating Pattern Extraction Technique (REPET): A
Simple Method for Music/Voice Separation" (Rafii / Pardo). This is a very simple method which works very well on "commercial" pop music (with a very stereotypical structure and modern computer-based production).
"Source/Filter Model for Main Melody Extraction From Polyphonic Audio Signals" (Durrieu et al). Complex algorithm, but the code is available.
If you want something really simple, just extract and track the most salient pitch and assume it is the lead melody. This works to some extent in recordings in which the singer voice is well isolated from the mix (other instruments EQ'ed so that they do not overlap much with it).

You don't need to bother here with resynthesis - just use the energy of the extracted spectrum / f0 track.

A last angle of attack is to use a multi-f0 estimator or a sinusoids tracker, and check if one of the tracks has frequency fluctuations typical of vibrato or glissando. In commercial pop music, most instruments have a very stable pitch contour (keyboards, fretted bass and guitars), and the voice is the only instrument with a very wobbly pitch contour - unless it has been doctored by autotune.


It's best to see the original text, but I assume $f$ here is the normalized frequency (actual frequency divided by the sample rate), so we have $\omega = 2 \pi f$. This means that if you need to evaluate the transfer function from DC ($f = 0$) to Nyquist frequency ($f = \frac{1}{2}$), you need to evaluate the transfer function with $\omega$ varying from 0 to $\pi$. Practically, you decide on a grid size $N$ (say $N=100$ points) and evaluate the transfer function at the points $\omega_k = \frac{k \pi}{N - 1}$ for $k = 0, 1, \ldots, N - 1$.

You make it sound like there are two independent variables but there is only one, $\omega$ - the other one, $f$ is just the same thing with different units.


What is causing my ANC LMS update to diverge?

Turns out it's just a scaling issue. Because the struck.unpack() gives out a signed 16-bit value, the step size is too large for it. Two ways to go around this, you scale down your input and scale back up when outputting, or you just adjust the step size mu.
I used a mu value of 10e-8 and seems like this value works differently in different environments, reasons unknown. With jupyter-notebook is around 10e-7, python idle is 10e-8, and Matlab is 10e-9. Just tweak around that range with the code I posted above and should give you a converged result.
Hope this helps whoever is reading.
edit; for MatLab, I scaled the input by 32768 so the mu value is 10e-9, if max amplitude is 1, use 0.1.


Why am I not getting the intended coefficients in this 2D Fourier demonstration?

Your misconception is to assume that you can compute Fourier coefficients with the fft. You cannot (really). Fourier coefficients are found by integrating a continuous time (periodic) function after multiplying it with a harmonic kernel. The cosine function gives two real-valued Fourier coefficients equal to half its amplitude in this case.

You are considering a matrix though which means that you have sampled your cosine function at a discrete grid (which is not a problem in your case) and truncated it to a finite width (which is exactly the problem you have). Then you are computing the FFT which is a fast way of computing a DFT (Discrete Fourier Transform). The DFT (of discrete vectors) is not the same as a Fourier transform (which is defined on continuous functions).

That said, your DFT can approximate Fourier coefficients, as the Fourier integral can be approximated by a sum, which is precisely the definition of the DFT then. The approximation errors that happen in this case are sampling and truncation.

For functions that are strictly bandlimited (like your cosine), sampling does not cause any error as long as you satisfy the Nyquist criterion. Conversely, for functions that are strictly time-limited, truncation does not cause any error as long as the whole function fits in your window. Unfortunately, it can be shown that no function can be time-limited and band-limited at the same time. For this reason, the DFT can never accurately predict the result of the Fourier transform.

There is hope though, since you're not after a Fourier transform but instead you're after Fourier coefficients. These are defined on periodic continuous-time functions. For these, the DFT can give you accurate predictions of the Fourier coefficients, as long as you make sure that the samples you are considering contain a whole period of the signal. This needs to be accurate, i.e., if you periodify your extracted window it should look like the original function. For instance, if the period of your function is 1 and you want to extract four samples, they need to be at 0, 0.25, 0.5, 0.75 (i.e., not 0, 1/3, 2/3, 1, which is a very common mistake in this setting).

Your Fourier coefficients will then be real-valued if and only if your function is symmetric, which in the sampled case means x[n] = x[N-n] for n=1, 2, ..., N-1.  If you sample your cosines correctly, this should work.


Instantaneous velocity and displacement from acceleration signal using a proper filtering method

There a a number of conceptual problems with you code.  The first has to do with filtering by zeroing out FFT bins.  This is covered in:

Why is it a bad idea to filter by zeroing out FFT bins?

The second has to do with just taking the real part of the fft.  I don't understand why you did this but this is the wrong way to do the wrong thing.

If your near term priority is to process your data.  I suggest you discard most of your code and take an entirely different approach with using the built in Matlab timeseries class and use one of the built in filter methods on the object you create. The plot methods are also easy to use and you can set units on the object you create.  This is the most direct way to a "correct" result with your data.

If you want to understand how to write your own correct code, this will take more time.  There is more than one way to filter data and choosing which takes some time to learn.  You can look at:

What resources are recommended for an introduction to DSP?

Essentially, I'm pointing you at a number of books that are several hundreds of pages of study, which is the long term answer to your problem. You will know what is correct.

If you do some Googling, you can probably find something in the middle that suits your needs but I can't recommend anything because I don't know what you know and don't know.

If you may permit the analogy, If one sees someone drowning, do you throw them a floatation device or try to teach them to swim.

Someone else here may provide a suitable solution.

It is to your credit that you have shown genuine effort.


Copy the spectrum to an array that is long enough to contain the spectrum of the double-length array. Put all the frequencies in their right places so that the zero frequency bin, the negative frequency bins (do not exist if you used real-input-FFT), and the positive frequency bins fall in their correct places. Zero the remaining (high-frequency) bins. If your original array was of even length, depending on the flavor of FFT, you may need to split the Nyquist frequency (highest original frequency, sampling frequency / 2) bin to its identical negative and positive frequency component and divide their value by two. You may further need to scale all bins by a constant (2 in the example) to account for the way the particular FFT and IFFT implementations are normalized. Then at the end use IFFT to get back to time domain (example in Octave):

octave:1> x = [1 2 3 4]
x =

   1   2   3   4

octave:2> X = fft(x)
X =

   10 +  0i   -2 +  2i   -2 +  0i   -2 -  2i

octave:4> Y = [10 -2+2i -1 0 0 0 -1 -2-2i]*2
Y =

 Columns 1 through 6:

   20 +  0i   -4 +  4i   -2 +  0i    0 +  0i    0 +  0i    0 +  0i

 Columns 7 and 8:

   -2 +  0i   -4 -  4i

octave:5> y = ifft(Y)
y =

   1.0000   1.0858   2.0000   2.5000   3.0000   3.9142   4.0000   2.5000

The process is not equivalent to linear interpolation but to ideal sinc interpolation with wraparound of the sinc function tails to account for the assumed periodicity of the sequence.