You still only get one range and one velocity line per radar observation.

You can, however,  "scan" a region (that's why airport towers and ships have these rotating antennas) and get another axis of information.

The same can be achieved using digital beamforming, which, however, requires multiple signal chains, and specific antenna arrays. For a two-dimensional picture, you'd need an at least two-dimensional array of antennas. For the resolution displayed, you'd need both high bandwidth (think in the order of 1GHz) as well as small wavelength (think of beyond 76 GHz). So, I'd say, without a millimeter-wave lab, this is probably not achievable.

I've not seen pictures like that myself (the article says "source: NXP" without elaborating any further, so this might just as well be "artist's interpretation of what NXP wants to achieve in 10 years).

Whilst this resolution might be impossible to achieve with commodity hardware, in principle, SDRs that are MIMO-capable do exist. If you can get your hand on a 4x4 MIMO phase-aligned SDR system, you could, behind/infront ambiguities notwithstanding, build a simple patch antenna array of 2 rows / 2 columns parallel patches and build a digital beamforming radar with that. You'd probably want to avoid the problems of isolating RX from TX by having a separate array for TX and RX.


Discrete Fourier Series of square of a signal

Here is an outline (in case it's homework)

Use a different frequency index for $X_2(l)$
Split the $X_2(l)$ into two sums: from $0:N-1$ and from $N-1:2N$
Use the periodicity of x[n], i.e.  $x(n) = x(n+N)$
Substitute $m = n-N$ in the second sum
Recombine the two sums again
Stare what happens to even $l = 2k$ and odd $l = 2k+1$ terms


How do I go about detecting whistles, pops and other sounds in live audio input?

While I agree with one of the comments that using speech recognition techniques might be a good start, these sounds are different and I'm not aware of anyone having done any research into categorizing them (the paper cited by Nathan appears only to distinguish between speech and noise), so unless someone presents something to contradict that, you are going to have to invent your own technique, and that's going to take a lot of learning and a lot of work. The best I can do is get you started.

First off, don't expect anyone to be able to produce some magic formula. There's no equation to get from sound to what the sound is. Humans and computers must learn by incorporating data to make educated guesses about what sounds are. The reason people say in their answers "use FFT" or "use zero crossing" is because these are some of the basic DSP building blocks used in speech recognition and related algorithms. But the FFT and zero crossing rate are usually just first steps in creating a set of parameters describing the sound. These parameters are then analyzed statistically (not through some magic function) to determine which category they most likely belong to. Note that I said "most likely": even the best speech detection (and the human brain!) has a high error rate at this stage and further processing using knowledge of the language and grammar is used to correct these mistakes.

So, some parameters you could look for include:

zero-crossing rate
frequency centroid
envelope (this is actually a set of parameters, including, eg attack time)
spectral envelope
even/odd harmonic balance
voiciness
fundamental pitch

Once you have a set of parameters that you think will allow you to distinguish your sounds, you will need to use some statistical method to categorize them. Hidden Markov Model is often used in speech. You can also look into logistic regression, K-means, and I'm sure there are other choices, but I think HMM is tried and true.


Mean Square Error and Gibbs oscillations

The first condition mentioned in your question is absolute summability, which is sufficient for the discrete-time Fourier transform (DTFT) to exist. In this case, the sum given by the DTFT of a sequence converges uniformly. The other condition you probably mean is square summability:

$$\sum_{n=-\infty}^{\infty}|x[n]|^2\lt\infty\tag{1}$$

in which case the DTFT exists if we relax the condition of uniform convergence. That other type of convergence is called mean-square convergence. In this case, the result of the sum oscillates around points of discontinuities, and those oscillations do not decrease with increasing number of elements in the sum. You can find a more in-depth explanation of this so-called Gibbs phenomenon in this article.

The importance of the Gibbs phenomenon in digital filter design is that you might want to approximate an ideal frequency response (with discontinuities) by a filter of finite length (finite impulse response, FIR), and this results in large errors close to the discontinuities. This problem can be alleviated by introducing transition bands and/or by (non-rectangular) windowing instead of simple truncation of the sum.


Simple FFT filtering vs. e.g. butterworth filtering

Edge effects due to circular convolution can occur if you don't zero pad the data sufficiently before FFT/IFFT filtering.  Fast convolution filtering using shorter sequential segments or windows from a longer data stream requires both padding of the FFT and overlap add/save processing of the overlapping IFFT results.

Just naively zeroing FFT bins instead of using the transform of a properly designed time domain filter can lead to extreme ripples in the filter's frequency response and/or ringing in the FFT/IFFT filtered result.

No phase distortion should occur if the filter kernel has appropriate symmetry in the time domain.

The computational cost of an FFT (in CPU cycles or latency time) can be higher than with IIR filtering.  But that depends on the length of the data and the filter complexity (in filter order and/or number of passes).


Questions in a Character Degradation Model Paper

All g is doing is saying whether the area is closer to a circle (less flat) or a line (completely flat ellipse).
Yes.
No.  $N(\mu, \sigma^2)$ is a distribution.  You can say
$$P_k \tilde{} N(\mu, \sigma^2)$$
read as $P_k$ is distributed normal, with mean $\mu$ and variance $\sigma^2$.

To find a value for $P_k$ you need to generate a random number with the mean given by equation (1) in the last screenshot:

$$
\mu = \bar{c}_i + (\bar{b}_j - \bar{c}_i) \left( \frac{d_{ik}}{d_{ij}} \right)
$$

In matlab, this would be something like normrnd(mu, sigma), though what sigma value to choose is not clear from your extractions from the paper.


Removing vocal but retaining all background ambient noises

What you are trying to archive is very hard if at all possible. The pinnacle of research today (including my own) is striving for good results on such tasks. DCASE challanges engage with similar tasks with some interesting results, though the make assumption which cannot be generalized to your case. In DCASE2019 task 3, for example, they assume up to 2 timewise parallel sources where in your case there are no such limitations.

With respect to your conditions and constraints, I would advise looking into source separation and\or classification. If your background audio has a repetitive nature, this paper may help.


Anybody with experience with online double integration of acceleration?

Any filter within your chain can remove useful information. The kalman filter way is applicable when you have/use a good model of the error, and you have some auxiliary sensor (like GPS, or images) to confront to the result of the double integration (DI). I dont know how good your acc is, but a military one, for example, can drift so few after the DI, that it doesnt need any auxiliary sensor. However, you cannot expect the same result using this cell-phone acc, though. The limit of goodness is defined by your application and the type of sensor you have.


Does mixing brown noise and white noise create pink noise?

No.

Pink noise has a spectrum of that falls with 3dB/octave (or 10dB/decade). The spectrum of the sum of white and brown noise will be "brown" at low frequencies and "white" at high frequencies. The spectrum will have two slopes: below the transition frequency it will be -6dB/octave and above it, it will just be flat.

Somewhere in the transition region will be a frequency range that could be considered approximately  "pink" but it's not really a constant slope.

That transition frequency is simply a function of the mix ratio (and the initial energies) of the two noise signals.

If you define "mixing" as a multiplication and not a summation, the result will just be white.


Extracted envelope via Hilbert Transform is tracking the signal

Computing the envelope of a signal via the Hilbert transform only works well for bandpass signals.

But apart from that, you're not taking into account the delay of the Hilbert transformer when creating the analytic signal. It would be easier to choose an odd filter length $N$, because then you get an integer delay:

$$d=\frac{N-1}{2}$$

Also note that there's the command hilbert for directly computing the analytic signal. Furthermore, there's a function envelope that lets you choose different methods for computing the envelope of a signal.


A 1D sinusoidal wave, time dependent, is described by the next formula:

$$
y(t) = A\sin(2\pi f t + \phi)
$$

The trick here is just the concept. When you use the formula $A\sin(\theta) $ you aren't specifying if $\theta$ is time-dependent and you are not thinking in the formula as a wave, but as a mere function. But when you use the argument $2\pi f t$ you are considering the angular frequency $2\pi f$ ($rad/s$), where $ f $ is the normal frequency, and the wave's phase $\phi$ (which represents a shift in the wave).

If you are working with a simple sine wave, you can assume a zero-shift, $\phi = 0$, and select a working frequency. Sample the wave and work normally with your DSP.

I hope it helps!

You can check sine wave description in WIKI


The piece I was missing was the distribution of the initial phase values $\varphi_1$ and $\varphi_2$. It is standard to assume that these are uniformly distributed [^]. This leads to:
$$
\mathbf{R_x} = \mathbf{E[xx^H]} = \mathbf{S \Lambda S^H}
$$
where
$\mathbf{\Lambda} = \begin{bmatrix} A_1^2 & 0 \\ 0 & A_2^2
\end{bmatrix}$
and $\mathbf{S} = \begin{bmatrix}\mathbf s_1 & \mathbf s_2 \end{bmatrix}$ so that
$$
\mathbf{R_y} = \mathbf{S \Lambda S^H} + \sigma^2 \mathbf{I}.
$$

Now note that $\mathbf{R_y}$ is Hermitian symmetric and must possess an eigendecomposition with eigenvalues $\lambda_1 \geq \lambda_2 \geq \lambda_3 \geq \cdots \geq \lambda_L$. However, also note that $\mathbf{R_x}$ has two nonzero eigenvalues, say, $\tilde\lambda_1 \geq \tilde\lambda_2$. Since $\mathbf{R_y} = \mathbf{R_x} + \sigma^2 I$ it is easy to verify that $\lambda_1 = \tilde \lambda_1 + \sigma^2$, $\lambda_2 = \tilde \lambda_2 + \sigma^2$ and $\lambda_i = \sigma^2$ for $3\leq i \leq L$.

In the eigendecomposition of $\mathbf{R_y}$, let $\mathbf e_1$ and $\mathbf e_2$ be the eigenvectors corresponding to $\lambda_1$ and $\lambda_2$ and $\mathbf g_1, \ldots, \mathbf g_{L-2}$ be the eigenvectors corresponding to the eigenvalues $\lambda_3, \lambda_4,\ldots, \lambda_L$. One can now show that $\mbox{span}\{\mathbf e_1, \mathbf e_2\} = \mbox{span}\{\mathbf s_1, \mathbf s_2\}$. This is often called the "signal subspace" and the orthogonal subspace $\mbox{span}\{ \mathbf g_1,\ldots \mathbf g_{L-2}\}$ is called the "noise subspace."

Here's the punchline: To estimate $\omega_1$ and $\omega_2$ we need to solve for $\omega$ in $\mathbf s(\omega)^H \mathbf G \mathbf G^H \mathbf s(\omega) = 0$, where $\mathbf G = \begin{bmatrix} \mathbf g_1 & \ldots & \mathbf g_{L-2}\end{bmatrix}$, and $\mathbf s(\omega) = \begin{bmatrix} 1 & e^{j\omega} & e^{j2\omega} & \ldots & e^{j(L-1)\omega} \end{bmatrix}$.

Reference

[^] P. Stoica, R. Moses, Spectral Analysis of Signals 1st Ed. Ch. 4.


When should the sum of all elements of a gaussian kernel be zero?

The sum of a gaussian kernel cannot be zero, because all the elements are going to be positive. The first kernel you have shown, is most likely an edge detection kernel, (which is a type of high pass filter), so the elements add up to zero because you want to completely null out any DC/constant component.

The second kernel you have shown however, is a low pass kernel, indeed used in smoothing. Here the elements are not going to add up to zero because you are not trying to remove the 'constant' part of the image, in fact, you are trying to remove any high frequency noise.

Hope that helps.


Evaluating the inverse $\mathcal Z$-transform on the unit circle

Hint:

With $z = e^{j\omega}$, you have:

\begin{align}
\displaystyle \frac{dz}{d\omega} &= je^{j\omega}\\
\Rightarrow dz &= je^{j\omega} d\omega
\end{align}

With $-\pi \leq \omega \leq \pi$, your contour integral goes from $-\pi$ to $\pi$.

  Can somebody tell me intuitively how to simplify this equation to look
  like the inverse Fourier transform?

Plugging all the above should do the job.

  Is my confusion due to a lack of understanding of the contour
  integral?

It is more the change of variables.

EDIT:

Note that $\omega$ is taken in the fundamental interval $\left[-\pi, \pi\right]$ because $X(\omega)$ is periodic with period $2\pi$. So, $\displaystyle \int_c$ becomes $\displaystyle \int_{2\pi}$


Solving a second order difference equation

You can indeed proceed with the complex roots $z_1$ and $z_2$ as well as with the real ones.

In order for a particular solution you may assume a solution of type:
$$y_p[n] = K ( \frac{1}{2})^n u[n]$$

And plug it into the LCCDE to find the undetermined coefficient $K=1/5$.

Then assuming a homogeneous solution of the form (for causal system) $$y_h[n] = C_1 z_1^n u[n] + C_2 z_2^n u[n]$$

You have the complete solution as:
$$y[n] = y_h[n] + y_p[n] = (C_1 z_1^n + C_2 z_2^n + \frac{1}{5} (\frac{1}{2})^n) u[n] $$

In order to find $C_1$ and $C_2$ you can get $2$ equations from $y[0]$ and $y[1]$ by explicitly solving them from the LCCDE via recursion, which yields
$y[0] = 1/3$ and $y[1] = 5/18$. Then you have the following:

$$
\begin{cases}
y[0] = \frac{1}{3} &= C_1 + C_2 + \frac{1}{5} \\
y[1] = \frac{5}{18} &= C_1 z_1 + C_2 z_2 + \frac{1}{10} \\
\end{cases}
$$

Solve this system for $C_1$ and $C_2$ to complete the procedure...

Alternateley you can also consider the following.
When have a LCCDE with real coefficients, its consequence is the fact that the roots associated with the homogeneous equation will be either real or in complex-conjugate pairs, as you have found.

Furthermore you can deduce from the LCDDE with real coefficients that its solution $y[n]$ must be real valued, for a real input $x[n]$.

Hence you can deduce from these that the constants $C_1$ and $C_2$ of those two complex roots should be constrained like $C = C_1 = C_2$ or $C = C_1 = C_2^*$  then you can see that the complex roots will merge:

$$ C_1 z_1^n + C_2 z_2^n = C (z_1^n + (z_1^n)^* ) = 2C \mathcal{Re} \{z_1^n\}$$

or like

$$ C_1 z_1^n + C_2 z_2^n = C z_1^n + (C_1z_1^n)^*  = 2 \mathcal{Re} \{C z_1^n\}$$

Conversion to polar form helps like:
$$z_1 = |z_1| e^{j\phi_1} \Rightarrow 2C \mathcal{Re}\{z_1^n\} =2C |z_1|^n\cos(\phi_1 n) $$

or like
$$C = |C|e^{j\theta_C}, ~~z_1 = |z_1| e^{j\phi_1} \longrightarrow 2 \mathcal{Re}\{C z_1^n\} =2 |C| |z_1|^n\cos(\phi_1 n + \theta_C) $$


Note that for any impulse response $h[n]$, the step response is given by the cumulative sum of $h[n]$:

$$a[n]=\sum_{k=-\infty}^nh[k]\tag{1}$$

Its $\mathcal{Z}$-transform is then

$$A(z)=\sum_{n=-\infty}^{\infty}a[n]z^{-n}=\sum_{n=-\infty}^{\infty}z^{-n}\sum_{k=-\infty}^nh[k]\tag{2}$$

The region of convergence (ROC) of $(2)$ is the range of values of $|z|$ for which the sum in $(2)$ converges.

For your example, you have three different impulse responses with the same expression for their $\mathcal{Z}$-transform, each with a different ROC. You are right that the ROC of $A(z)$ is the intersection of the ROC of $h[n]$ and the ROC of the unit step $u[n]$.

For the ROC $|z|<\frac12$ you get an empty intersection. This simply means that the sum $(2)$ converges nowhere in the $z$-plane, i.e. the $\mathcal{Z}$-transform of the step response $a[n]$ doesn't exist in this case. The reason is that the corresponding impulse response is anti-causal and exponentially increasing towards negative values of $n$, so from $(1)$, $a[n]$ is infinite for any value of $n$.

For the ROC $\frac12<|z|<2$ you have the intersection $1<|z|<2$, i.e., the step response $a[n]$ is a two-sided sequence with a $\mathcal{Z}$-transform $A(z)$ that converges in that range of values of $|z|$.

For the ROC $|z|>2$ the intersection is $|z|>2$, and, consequently, the step response $a[n]$ is an exponentially increasing right-sided sequence with a $\mathcal{Z}$-transform $A(z)$ that converges for $|z|>2$.


Digital filter coefficients from low-pass to high-pass

You can apply a so-called all-pass transformation to a discrete-time low-pass prototype filter in order to convert it to other standard filters (such as high-pass, band-pass, and band-stop). This is accomplished by transforming the complex variable $z$ in the transfer function of the prototype filter by a function $G(z)$ which satisfies $|G(e^{j\omega})|=1$, i.e., $G(z)$ is an all-pass function. This makes sure that the transformation maps the unit circle onto itself, i.e., the frequency response of the new filter is just a shifted and/or warped version of the prototype frequency response.
The most straightforward way to transform a discrete-time low-pass filter to a high-pass filter is to use the trivial transform $G(z)=-z$, i.e.,
$$H_{HP}(z)=H_{LP}(-z)\tag{1}$$
If $H_{LP}(z)$ is given by
$$H_{LP}(z)=\frac{\displaystyle\sum_{k=0}^Nb[k]z^{-k}}{\displaystyle 1+\sum_{k=1}^{N}a[k]z^{-k}}\tag{2}$$
then $H_{HP}(z)$ becomes
$$H_{HP}(z)=\frac{\displaystyle\sum_{k=0}^N(-1)^kb[k]z^{-k}}{\displaystyle 1+\sum_{k=1}^{N}(-1)^ka[k]z^{-k}}\tag{3}$$
This transformation shifts the low-pass frequency response by $\pi$ (i.e., by half the sampling frequency). Consequently, if $\omega_c$ is the cut-off frequency of the low-pass prototype filter, the cut-off frequency of the resulting high-pass filter is given by $\omega'_c=\pi-\omega_c$.
Other cut-off frequencies can be obtained by applying the more general lowpass-to-highpass transformation
$$G(z)=-\frac{z+\alpha}{1+\alpha z},\qquad |\alpha|<1\tag{4}$$
The simple transformation shown above is obtained from $(4)$ by the choice $\alpha=0$.
This and other frequency transformations applicable to discrete-time filters are treated in some detail in Chapter $7.4$ of the third edition of Oppenheim and Schafer's Discrete-Time Signal Processing.
Also take a look at these related questions and their answers: Q1, Q2.


You don't need to know the phase at all and you only need to know the frequency approximately. It is the job of your digital carrier recovery scheme to estimate and remove the phase and residual frequency offset.


How to calculate the data rate of IEEE 802.15.4 to be 250kbps?

You have $16$-ary symbols, each symbol requires $4$ bits. With a chip rate of $1000$ K, the symbol rate is $1000/16=62.5$ K, hence the bit rate is $4\times 62.5=250$ Kbps. Did you mean this calculation?


Biquad Filter Transfer Function Frequency Response

To get the frequency response, replace $z$ in $H(z)$ by $e^{j2\pi f/f_s}=\cos(2\pi f/f_s)+j\sin(2\pi f/f_s)$ where $f_s$ is the sampling frequency in Hz and $f$ is the desired frequency in Hz. This will give you the value of the frequency response at frequency $f$. Note that the frequency response is a complex function, so you might want to look at its magnitude.


Finding Fourier series coefficients for discrete time signal

You are asking for DFS (discrete Fourier series) coefficients $b_k$ of the periodic sequence $y[n]$, in terms $a_k$ of the periodic sequence $x[n]$. Since DFS and DFT (discrete Fourier transform) are the same things, I will instead write down an answer for the DFTs $Y[k]$ and $X[k]$ of the associated sequences, in the context of which $X[k] = a_k$ , and $Y[k] = b_k$ will be understood.
I would like to relate DFTs to DTFTs (discrete-time Fourier transform) through the sampling relation, which indicates that :

N-point DFT X[k] of the N-point sequence x[n], is the uniform samples of the DTFT X(w) of
x[n].

The DTFT of $x[n]$ is :
$$ X(\omega) = \sum_{n=0}^{N-1} x[n] e^{-j \omega n } \tag{a}$$
and the DFT of $x[n]$ is :
$$ X[k] = X(\frac{2\pi}{N}k) = \sum_{n=0}^{N-1} x[n] e^{-j \frac{2\pi}{N}nk } \tag{b}$$
Let $x[n]$ and $y[n]$ be two sequences related by:
$$y[n]=\begin{cases}{ x[n] ~~~,~~~n:even \\  ~~~0 ~~~~~,~~~ n:odd }\end{cases}$$
You can also express $y[n]$ as :
$$ x[n] \longrightarrow \boxed{ \downarrow 2 } \overset{w[n]}{\longrightarrow} \boxed{ \uparrow 2 }    \longrightarrow y[n] $$
Applying DTFT relations to the sequences :
$$ W(\omega) = 0.5 \left( X\left( \frac{\omega}{2} \right) +  X\left( \frac{\omega-2\pi}{2} \right) \right) \tag{1} $$
and
$$ Y(\omega) = W(2\omega) \tag{2} $$
hence
$$ Y(\omega) = 0.5 \left( X\left( \omega \right) +  X\left( \omega-\pi\right) \right) \tag{3} $$
Eq.3 describes the relation (for any $N$) between the DTFTs $Y(\omega)$ and $X(\omega)$ of the sequences $x[n]$, and $y[n]$. To turn this into the equivalent relation between the N-point DFTs $X[k]$ and $Y[k]$, (and for your case a relation between DFS  $a_k$ and $b_k$) we shall sample the associated DTFTs. Then we have :
$$  \begin{align}
Y[k] &= Y( \frac{2\pi}{N} k ) ~~~,~~~k = 0,1,...,N-1 \tag{4} \\ \\
     &= 0.5 \left( X\left( \frac{2\pi}{N} k \right) +  X\left( \frac{2\pi}{N} k-\pi\right) \right) \tag{5} \\ \\
     &\boxed{ Y[k] = 0.5 \left( X\left[ k \right] +  X\left[k- \frac{N}{2} \right] \right) }\tag{6} \\
\end{align}$$
Eq.6 describes the requested relation between N-point DFTs $X[k]$ and $Y[k]$.
When $N$ is EVEN (say $N = 2 M$ ) the relation becomes :
$$ \boxed{ Y[k] = 0.5 ( X[k] + X[k-M] ) }\tag{7} $$
which is what you have arrived at $b_k = 0.5(a_k + a_{k -\frac{N}{2}})$.
However, when $N$ is ODD, there is a problem due to $N/2$ not being an integer. The sequence $X[k-N/2]$ is interpreted as an interpolation of $X[k]$ at the fractional indices $k-N/2$.
For $N=2M+1$ , we have $X[k-N/2] = X[k-M-0.5]$ , $k=0,1,...,N-1$ , which evaluates $X[k]$ at the points $m = 0.5,1.5,...,N-0.5$. To get those intermediate samples, we need an interpolation of $X[k]$ by 2.
Interpolation of $X[k]$ by 2, is achieved by zero padding $x[n]$ by $N$ samples, and computing its $2N$-point DFT :
$$ x_e[n] = \begin{cases}{ x[n] ~~~,~~~ 0 \leq n < N \\ ~~~0 ~~~~ , ~~~~ N \leq n < 2N }\end{cases} $$
Then we can see, by the DFT - DTFT sampling relation that :
$$ \mathcal{DFT}\{x_e[n]\} =X_2[k] = X(\frac{2\pi}{2N}k )= X(\frac{2\pi}{N} k/2 ) = X[k/2] \tag{8} $$
for $ k = 0,1,.,2N-1$. Note that DTFTs of $x[n]$ and $x_e[n]$ are the same.
Based on Eqs. 8, you can restate Eq.6 when $N$ is odd :
$$ \boxed{ Y[k]= 0.5 \big( X\left[ k \right] +  X_2\left[2k- N\right] \big) } \tag{9} $$
with your notation it becomes:
$$ b_k = 0.5 ( a_k + a^e_{2k-N} ) $$ where $a^e$ is the DFS of the sequence $x_e[n]$.
As you can see, you cannot express $b_k$ (with a simple formula) in terms of sole $a_k$ when $N$ is odd.


Combining time domain and frequency domain data for BPM detection in real time

BPM detectors work by finding the strongest period in the onset detection function - a function of time indicating how likely it is that there is an onset at time t. Which brings the topic of onset detection...

Tracking energy/envelope in the time domain is a weak strategy for detecting rhythm/notes, because:

For some instruments, there are playing techniques which create new notes without a significant change of energy (for example moving the finger on a string of a cello with no change in bowing, or changing the fingering on a flute with continuous blowing). This has the effect of creating changes in fundamental frequency without a noticeable change in the time-domain envelope.
A lot of commercial music is so heavily compressed that the signal is almost "pumped" to maximum amplitude all the time.
While this is yet another corner case, it is worth reminding that some instruments or playing techniques have very slow attacks.

So naively looking for peaks in the time-domain (via energy) will not be very robust. You have to look for better "clues", and the frequency domain is the answer. The simplest thing that can be done in the frequency domain is to compute the difference of energy between adjacent STFT frames (a measure known as the "spectral energy flux").

A more elaborate metric, which is actually as cheap in terms of computing costs, consists in checking the deviation between the actual and expected value of the complex amplitude in a STFT bin (reference here) - this deviation would be null in the stable segment of a note (as the amplitudes are just sustaining or slowly decaying; and the phases are just keeping on rolling), and high at an onset.

All these techniques have very modest computational costs and can be implemented in realtime even on relatively cheap microcontrollers - a complete beat-tracker (onset detection + BPM detection + beat tracking) is in the 10 MIPS range.


If $x(t)$ is a finite-energy signal with Fourier transform $X(f)$, then $x(t)\cos(2\pi f_c t)$ is also a finite-energy signal with Fourier transform $\left.\left.\frac 12 \right[X(f-f_c) + X(f+f_c)\right]$. This is just the modulation theorem of Fourier transform theory.

The energy spectral density of $x(t)$ is $S_x(f) = |X(f)|^2$ while the energy spectral density of $y(t) = x(t)\cos(2\pi f_c t)$ is
\begin{align}S_y(f) &= \bigg|\left.\left.\frac 12 \right[X(f-f_c) + X(f+f_c)\right]\bigg|^2\\
&= \left.\left.\frac 14 \right[ |X(f-f_c)|^2 + |X(f+f_c)|^2 + 2\Re\big(X(f-f_c)X^*(f+f_c)\big)\right]\\
&= \left.\left.\frac 14 \right[S_x(f-f_c)+S_x(f+f_c)+2\Re\big(X(f-f_c)X^*(f+f_c)\big)\right]\tag{1}
\end{align}
which has a vague resemblance to what the OP wrote.  Now, if $x(t)$ is a low-pass signal with bandwidth $W$ Hz and $W < f_c$, then $X(f-f_c)$ occupies the frequency band from $f_c-W > 0$ to $f_c+W$ while $X(f+f_c)$ occupies the frequency band from $-f_c-W$ to $-f_c+W < 0$. We see that there is no value of $f$ for which $X(f-f_c)$ and $X(f+f_c)$ both are nonzero, and thus $X(f-f_c)X^*(f+f_c)$ equals $0$ for all $f$. We are left with
\begin{align}
S_y(f) &= \left.\left.\frac 14 \right[S_x(f-f_c)+S_x(f+f_c)\right]\\
&= \left.\left.\frac 14\right[|X(f-f_c)|^2+|X(f+f_c)|^2\right].\tag{2}\end{align}
But when $f_c$ is smaller than $W$, then, as the OP peevishly has informed us all that he has placed no restrictions on the value of $f_c$, he must learn to live with the more complicated $(1)$.


Linear Time-Invariant system frequency response and impulse response given input and output in time domain

A classical question on signals and systems, for which I provide you one method of solution.

Now given an LTI system, and its input $x[n]$ and output $y[n]$, the fundamental approach to determine the impulse response $h[n]$ involves the utilisation of a frequency domain approach to assert $H(e^{jw}) = Y(e^{jw})/X(e^{jw})$

In this problem however, I would choose a different approach which involves more of an analysis by inspection rather than pure mathematical operations, because of the fact that the input $x[n]$ involves a cosine term whose DTFT is includes impulses and will be problematic in the division process.

The method proceeds by drawing the DTFTs of input and output:
The input $x[n]$ involves two terms one of which is $$5 \sin(0.4 \pi n)/ \pi n$$ whose DTFT is a zero phase ideal lowpass filter of cutoff frequency $W_c = 0.4 \pi$ with a gain of 5.

And the second term of input is $$10 \cos (0.5 \pi n)$$ whose DTFT is modeled as two impulses at $$ 5 \pi \delta(\omega + 0.5 \pi) +  5 \pi \delta(\omega - 0.5 \pi) $$ Because of linearity property of DTFT, these two DTFT are added to represent the DTFT of the input $x[n]$.

Then, look at the DTFT of the output $$y[n] = \frac {10 \sin(0.3 \pi (n-10) )}{\pi (n-10)} $$ which is again an ideal low pass filter of gain = 10 and $W_c = 0.3 \pi$ with a phase term of $e^{-10 \omega}$ hence $Y(e^{j\omega}) =10 e^{-j10 \omega}$ when $|\omega| < 0.3 \pi$ and zero elsewhere.

Now we want to find a frequency response $H(e^{j\omega})$ such that when it is multiplied by the DTFT $X(e^{j\omega})$ of the input $x[n]$, it should produce the DTFT $Y(e^{j\omega})$ of the output $y[n]$.

The observation of the input and output DTFTs yields the conclusion that one such $H(e^{j \omega})$ is a linear phase ideal lowpass filter of gain = 2, $W_c = 0.3\pi$ and phase term of $e^{-j10 \omega}$, whose time domain correspondence is found by inspection as $$h[n] = \frac {2 \sin(0.3 \pi (n-10))}{\pi (n-10)}$$

The justification of this solution can be made by observing the effect of this filter on its input: it will completely supress any signal above the frequency $\omega = 0.3 \pi$ and it will multiply any signal below that frequency by 2, while adding a phase shift of -10 radians which reflects itself as a group delay of 10 samples in (n-10) term.


How to decide filter order in Linear Prediction Coefficients ( LPC ) while calculating formant frequency features?

You have to define the order based on your application, usually for speech you'll use and order of 10, this definws the order of your all pole filter and it's related with the amount of formants tha you will see in the impulse response of the filter.

https://www.dsprelated.com/showthread/speechcoding/965-1.php


Why do we express complex noise as Gaussian distributed while the channel is Rayleigh distributed?

There are too many things named after Rayleigh floating around in this problem.

Consider first additive white Gaussian noise with two-sided power spectral density $\frac{N_0}{2}$.  This is generally taken as the model for channel noise, though in fact the source of the noise actually is thermal noise in the front-end of the receiver. NASA once even experimented with cooling the receiver front end in a liquid helium bath to reduce the noise temperature to about $5^\circ$K. Be that as is may,
the receiver has filters (and amplifiers) so that the noise at the receiver output is definitely band-limited. Because the front-end noise has been amplified greatly, the contribution of the impedances of later receiver stages to the thermal noise is negligible. Indeed, the noise at later stages is a
(stationary) band-pass Gaussian process $\left\{N(t)\colon -\infty < t < \infty\right\}$ where
$$N(t) = X(t)\cos(2\pi f_ct) - Y(t)\sin(2\pi f_ct)\tag{1}$$ with $\left\{X(t)\colon -\infty < t < \infty\right\}$ and $\left\{Y(t)\colon -\infty < t < \infty\right\}$ being independent zero-mean low-pass Gaussian processes with common autocorrelation function $R(\tau)$. And Yes, I do want a $-$ sign in $(1)$.

The autocorrelation function of the $\{N(t)\}$ process is
\begin{align}\require{cancel}R_N(
\tau) &= E[N(t)N(t+\tau)]\\
&= E\big[\big(X(t)\cos(2\pi f_ct) - Y(t)\sin(2\pi f_ct)\big)\\
&= \hspace{0.2in}\cdot
\big(X(t+\tau)\cos(2\pi f_c(t+\tau)) - Y(t+\tau)\sin(2\pi f_c(t+\tau))\big)\big]\\
&= R(\tau)\big[\cos(2\pi f_c t)\cos(2\pi f_c(t+\tau))+\sin(2\pi f_c t)\sin(2\pi f_c(t+\tau))\big]\\
&= \hspace{0.2in}~ - \cancelto{0}{R_{X,Y}(t)}\big[\cos(2\pi f_c t)\sin(2\pi f_c(t+\tau))-\sin(2\pi f_c t)\cos(2\pi f_c(t+\tau))\big]\\
\\&= R(\tau)\cos(2\pi f_c \tau)
\end{align} which is a band-pass function.  Note that we can also express the $\{N(t)\}$ process in terms of its complex baseband equivalent process $\{\mathcal Z(t)\colon -\infty < t < \infty\}$ where
$\mathcal Z(t) = X(t) + jY(t)$ as
$$N(t) = \mathrm{Re}\left(
\mathcal Z(t)e^{j2\pi f_c t}\right)\tag{2}$$ But, if we choose to represent the  complex number $\mathcal Z(t)$ in polar form as $B(t)e^{j\Theta(t)}$ where $B(t) \geq 0$ and $\Theta(t) \in [0,2\pi)$ are real numbers, then
$B(t) = \sqrt{X^2(t)+Y^2(t)}$ is a Rayleigh random variable and $\Theta(t)$ is uniformly distributed on $[0,2\pi)$. (This idea is the basis of the Box-Muller method of simulating Gaussian random variables. The usual probability transform method does not work well for Gaussian random variables because the Gaussian CDF does not have an inverse that is an elementary function whereas a Rayleigh CDF $1-\exp(-t^2/2\sigma^2)$ is readily invertible. So, the Box-Muller method calls a random number generator twice to simulate $B$ and $\Theta$ and then sets $X = B\cos(\Theta)$ and $Y = B\sin(\Theta)$ to get realizations of two independent Gaussian random variables).

The other place where Rayleigh's name comes in is in the transmission of signals over fading channels. Real-life signals are real-valued, not complex-valued.  When a sinusoid (say $\cos(2\pi f_ct), -\infty < t < \infty$) is transmitted over a (slowly fading) channel, the received signal is of form
\begin{align}r(t) &=  b(t)\cos(2\pi f_c t + \theta(t)) \tag{3}\\
&= b(t)\cos(\theta(t))\cos(2\pi f_c t) - b(t)\sin(\theta(t))\sin(2\pi f_c t) \tag{4}\end{align}
where $b(t)$ and $\theta(t)$ are slowly varying signals. For the purposes of digital communications, it is reasonable to treat $b(t)$ and $\theta(t)$ as constants over the duration of transmission of one bit/symbol, so that $(4)$ simplifies to
$$b\cos(\theta)\cos(2\pi f_c t) - b\sin(\theta)\sin(2\pi f_c t).\tag{5}$$
Slow fading occurs in ionospheric channels because the transmitted signal
is reflected back to earth from the ions in a ion cloud and the received signal is a sum of a very large number of very weak reflections with nearly equal delays. As a consequence, with a lot of handwaving and genuflections in the direction of the Central Limit Theorem, it is possible to regard
$b\cos(\theta)$ and $b\sin(\theta)$ in $(5)$ as realizations of two independent Gaussian random variables which allows us to treat $b$ and
$\theta$ as realizations of Rayleigh and uniform random variables respectively.  Hence, the complex baseband equivalent of a Rayleigh fading channel is that complex input $\mathrm{Re}(\mathcal Ae^{j2\pi f_c t})$ is received as output $\mathrm{Re}(\mathcal{AB}e^{j2\pi f_c t})$ where
$\mathcal B = Be^{j\Theta}$ with $B$ and $\Theta$ being independent
Rayleigh and uniform random variables respectively.  Fading is sometimes referred to as multiplicative noise since the fading multiplies  the input amplitude instead of adding to it as additive noise does.

"But, but, but," you say, "In a Rayleigh channel model, the fading parameters $B_{n+1}$ and $\Theta_{n+1}$ for the next bit/symbol interval are independent of the parameters $B_n$ and $\Theta_n$ for the current interval. How does that fit into your model of $b(t)$ and $\theta(t)$ being slowly varying functions? If $b(t)$ and $\theta(t)$ really vary so slowly that you can treat them as constant over the interval, then $B_{n+1}$ and $B_{n}$ should be very highly correlated, as should with $\Theta_{n+1}$ and $\Theta_n$."  Well, the dirty little secret that no one will ever tell you about is that communication systems intended for use over Rayleigh fading channels are very heavily interleaved so that what appears to be the next symbol in the mathematical model is actually a symbol that occurs $N\gg 1$ bit durations later where $N$ is chosen to be large enough that heavily-faded consecutive symbols appear in different frames rather than wiping out four bars of
Bethoven's Fifth Symphony or a whole stripe from an image of Old Glory if interleaving were not used.


Autocorrelation of filtered cyclostationary random process

Your first two equations generalize to the spectral correlation and cyclic autocorrelation. The most general results are (12) and (13) on my blog post on input-output relations for $n$th-order cyclic cumulants and $n$th-order cyclic polyspectra: signal processing operations and cyclostationary signals.

You can find the generalization of your first equation in several places, including Professor Gardner's book "Statistical Spectral Analysis", Eq (90) in Chapter 11:

$$S_y^\alpha(f) = H(f+\alpha/2) H^*(f-\alpha/2) S_x^\alpha(f)$$

And yes, the output of the filter is also cyclostationary with the same period provided that the transfer function $H(f)$ is such that $H(f+\alpha/2)H^*(f-\alpha/2)$ is not zero for at least some $f$ for which $S_x^\alpha(f)$ is also not zero. In other words, the filter could remove so much of the energy of the signal that the remainder is no longer cyclostationary.


resources for learning more about pedestrian tracking

What I have so far:

A Survey of Human-Sensing: Methods for Detecting Presence, Count, Location, Track, and Identity; Teixeira, Dublon and Savvides
I just found the question People Detection and Tracking on SO which also has some good leads.


A few points first

Since you are using digital filters, you should be using residuez() not residue()
Matlab's residuez() is not particularly well implemented: in case of higher order or steep filters, you can get a lot of numerically noise or results that are just totally wrong. Always check the result against your original target.

Back to the question: each residue comes as a conjugate complex pair (at least if you are lucky and don't have to deal with real or multiple poles). Just add these to form a single real second order section.

$$H(z) = \frac{r}{1-p \cdot z^{-1}} + \frac{r^{*}}{1-p^* \cdot z^{-1}} $$
$$= \frac{r \cdot (1-p^* \cdot z^{-1})+r^* \cdot (1-p^* \cdot z^{-1})}{(1-p \cdot z^{-1}) \cdot (1-p^* \cdot z^{-1})}$$

$$= \frac{(r + r^*) - (r \cdot p^* + r^* \cdot p  ) \cdot z^{-1} }{1 -(p+p^*)\cdot z^{-1}+ p \cdot p^* \cdot z^{-2}}$$

The last equation  gives you directly the real valued biquad coefficients (if I did the math correctly,  that is).


I am not giving you a complete algorithm here, but since I worked on a similar project, I can give you some hints and tips.

First of all, changing an image taken from one perspective to a different one relatively easy only for planar surfaces. E.g., you have an image of a tall building with a road, and since buildings are usually build vertical in the air and roads tend to be flat on the ground, the planes of the building and the road are perpendicular. Now, if you want to get an image where you are looking straight in to the building (e.g. building plane parallel to image plane), you have two options:

either transform the image so that the building plane is the way you want it, and disregard the rest of the image (road plane will not be correctly transformed)
in addition to transforming the building plane, calculate and perform a transformation on the image part representing the road plane, and somehow stitch them together.

The transformation that allows you to "change" the perspective is called the planar perspective transformation or homography. A homography represents a (mathematical) connection between two images of a same object taken from different perspective, e.g. allows you to calculate the unknown image coordinates of a point in space in one of the images from known image coordinates of the same point in the image with a different perspective.

To be able to calculate a homography, you need to know four pairs of corresponding points. In case of bad camera, you will also have to take care of the camera intrinsic parameters. A very good post on this can be found here.

Now, if you really want to do it properly, you should understand all the terms I used up until now really well before doing anything else. That said, I have been working on a fairly large project, and some parts required to transform images of a road taken from a moving car in to images taken from the bird's perspective, which is what you want to do now. The project was dealing with various road signalization, traffic signs and road markings. Useful references would be:

MASTIF project site with all information related to the project
publication list for the project containing all published articles related to the project
Petar Palašek, Petra Bosilj, Siniša Šegvić: "Detecting and recognizing centerlines as parabolic sections of the steerable filter response" - a paper where I am a coauthor, shows a practical use for what you're trying to achieve, but is also a good source of further references. Has a brief explanation of homography - Section III should be the most interesting.
Ivan Sikirić, Karla Brkić, Siniša Šegvić: Recovering a comprehensive road appearance mosaic from video - a paper detailing how to calculate a rectifying transform used to change the image perspective.

This should give you all the material you need to get what you want to do. Good luck, it took us a while.


Calculating phase shift between two sines of equal, unknown, amplitude and frequency

Two samples is sufficient and we can determine if the phase difference is indeed $\frac{3\pi}{2}$ as follows:

Start with the hypothesis that the phase difference $\phi_1-\phi_2$ is $\frac{3\pi}{2}$, which is equivalent to $-\frac{\pi}{2}$

If and only if the phase is in such quadrature, then $s_1(t)$ and $s_2(t)$ will be the real and imaginary components of a single vector on the unit circle with radius A which can be given as $-jAe^{j(t+\phi_1)} = As_1(t) + jAs_2(t)$. This  starts at the angle $\phi_1-\pi/2$ when t = 0 and rotates counter-clock-wise with increasing t, which matches the plot shown with the purple line representing the real axis and the red line representing the imaginary.

Starting from a position on the unit circle at time = $t_1$ with position given by the real and imaginary components $s_1(t_1)$ and $s_2(t_1)$, we can move forward in time to $t_2$ over a determined angle which should predict $s_1(t_2)$ and $s_2(t_2)$. If this prediction matches the given result, then this will confirm the hypothesis:

The starting angle using the first sample as given by B and D is given by $tan^{-1}(s_1(t_1)/s_2(t_1))$
$= tan^{-1}(3.614015/1.174265) = 1.256637 $ radians.

The change is angle to the second sample is simply the difference in time since the frequency unit is 1 radian/sec given by $sin(t)$:

$\Delta\phi = t_2-t_1 = 231-230 = 1$ radian    (NOTE! This alone clearly does not match the plot visually, so the numbers given cannot be actually derived from the plot! There must be a frequency factor multiplied by t in the actual formulas for this to match the plot).

Therefore if the phase difference was $-\frac{\pi}{2}$ then the second angle would be at $1.256637 + 1$ radian or $2.256637$ radians and the ratio of $(s_1(t_2)/s_2(t_2))$ would be $\tan(2.256637)= -1.22195$

This does not match the ratio for the second sample given by A and C and therefore the phase difference cannot be $\frac{3\pi}{2}$ for the formulas as given.


Finding preferred pair of $m$-sequence for Gold code

Printed tables of irreducible polynomials of large degrees are usually incomplete because there are so many of them!  What one can do is use various tricks of the trade to figure out a set of $68 = 2\times 34$ consecutive bits of one m-sequence from the other m-sequence (with known generator polynomial as taken from the Peterson and Weldon table) and then run the Berlekamp-Massey algorithm on these consecutive bits to get the shortest shift register that generates the desired m-sequence.
If ${\bf s} = s_0, s_1, s_2, \ldots $ is one of the m-sequences, then the other m-sequence ${\bf t} = t_0, t_1, t_2, \ldots $ is related to ${\bf s}$ via
$t_i = s_{qi}$ where $q = 2^{16}+1 = 65,537$. See the paper "Cross-correlation properties of pseudorandom and related sequences," Proc. IEEE, vol.68, pp.593-619, May 1980 (unfortunately behind IEEE's paywall, but copies can be found in the Internet) if you are unfamiliar with this idea. So, we crank up the LFSR that generates ${\bf s}$ and run it to produce the sequence ${\bf s}$, recording $$t_0 = s_0, ~~t_1 = s_{65537},~~ t_2 = s_{131074},~~ \cdots, t_{67} = s_{4390979}.$$  Once upon a time, cranking out 4 million+ bits would have taken some effort, but these days, it's no big deal.
What next? Well, the Berlekamp-Massey algorithm can be viewed as a shift-register synthesis algorithm: it finds the shortest LFSR that can generate any given sequence. Thus, applying the Berlekamp-Massey algorithm to the first $68$ bits of $\bf t$ which we have just found gives us the LFSR that generates the m-sequence $\bf t$, and away we go!!


The performance of line detectors are image and application dependent in generally. Typically, people tune the parameters of the detectors to get the best possible result. When such a result cannot be obtained, it is often because of the insufficient performance of the detector for the current problem. One approach would be to change the algorithm. For instance, OpenCV also has an LSD line segment detector. Have you tried that? Maybe it can also be used in combination with Hough lines.

Having said all those, I would also look at the hardware/system related issues, if they are under your control: Smartly selecting system parameters such as lighting, cameras and lenses generally solve half of the problem. If you would like to operate in the wild, then I would suggest more sophisticated software solutions.

Edit: This question looks very similar to the one recently posted:
MatLab - Can not detect all lines in a "simple" image using houghlines

Maybe this could help you out.


Interpolating the peak of a cross-correlation using a centroid

There are a lot of papers on bin interpolation and some of them have obvious connections to interpolating a parabola and others aren't' so obvious. Given what you presented, I'm sure that you can make the connection between cross correlation and spectral analysis. I suggest  you look at:

M. D. Macleod, "Fast nearly ML estimation of the parameters of real or complex single tones or resolved multiple tones," in IEEE Transactions on Signal Processing, vol. 46, no. 1, pp. 141-148, Jan 1998.
doi: 10.1109/78.651200
Abstract: This paper presents new computationally efficient algorithms for estimating the parameters (frequency, amplitude, and phase) of one or more real tones (sinusoids) or complex tones (cisoids) in noise from a block of N uniformly spaced samples. The first algorithm is an interpolator that uses the peak sample in the discrete Fourier spectrum (DFS) of the data and its two neighbors. We derive Cramer-Rao bounds (CRBs) for such interpolators and show that they are very close to the CRB's for the maximum likelihood (ML) estimator. The new algorithm almost reaches these bounds. A second algorithm uses the five DFS samples centered on the peak to produce estimates even closer to ML. Enhancements are presented that maintain nearly ML performance for small values of N. For multiple complex tones with frequency separations of at least 4π/N rad/sample, unbiased estimates are obtained by incorporating the new single-tone estimators into an iterative “cyclic descent” algorithm, which is a computationally cheap nonlinear optimization. Single or multiple real tones are handled in the same way. The new algorithms are immune to nonzero mean signals and (provided N is large) remain near-optimal in colored and non-Gaussian noise
keywords: {Fourier analysis;Gaussian noise;amplitude estimation;frequency estimation;harmonic analysis;interpolation;iterative methods;maximum likelihood estimation;optimisation;phase estimation;signal sampling;spectral analysis;white noise;AWGN;Cramer-Rao bounds;FFT;amplitude;cisoids;colored noise;complex single tones;discrete Fourier spectrum;fast nearly ML estimation;frequency;frequency separations;interpolator;iterative cyclic descent algorithm;noise;nonGaussian noise;nonlinear optimization;nonzero mean signals;parameter estimation;peak sample;phase;real single tones;resolved multiple tones;sinusoids;unbiased estimates;uniformly spaced samples;Amplitude estimation;Computational efficiency;Frequency estimation;Iterative algorithms;Maximum likelihood estimation;Noise level;Parameter estimation;Phase estimation;Phase noise;Signal processing algorithms},
URL: http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=651200&isnumber=14197

and if you can get into IEEE Explore, look at the papers that cite it. These are the State-OF-The-Art bin interpolation schemes, which you should be making your comparisons to.

The only comment that I have about you centroid technique, would be that the accuracy would seem to depend on there being only a single delay in your data.  My experience is that where there is one delay, there are often many.  Having too wide a window would seem to be more susceptible to additional delays.


Realize Simulink block "Ideal ADC Quantizer" in Matlab

Resolved.
The problem was in erroneous quantization.
Here is the correct quantization code:

% Quantization
dy = (2*A)/(2^res-1);
partition = (-A+A/2^res:dy:A-A/2^res);
codebook = (0:2^res-1);
[index,yq] = quantiz(y,partition,codebook);
yq = 2*A*yq/(2^res-1)-A;

where:

y - SIN waveform values

A - amplitude of SIN waveform

res - resolution of quantizer


First of all, you should totally fix the $x$-axis :)

Now, as you only have two sine waves, you should expect to have peaks at their exact frequencies (that is, at $\pm 60\textrm{ Hz}$ and $\pm 100\textrm{ Hz}$). The FFT function in MATLAB gives you the Discrete Fourier Transform of your signal within $0$ and $F_s$, where $F_s$ is your sampling frequency -and given that you represent your signal with samples, even if you are not aware of it, you have one. Anyway, to make visualization easier, we focus on the interval $[-\frac{F_s}{2},\frac{F_s}{2}]$; that way, your spectrum will be centered around the $0$ frequency -use the function fftshift for that. So, now, once you fix the $x$-axis according to your sampling frequency, you will find that the peaks are at the frequencies you have defined for your sine waves.

When it comes to the little ripples in the lower part of your spectrum, that's just disturbance caused by the noise. Since you are adding noise to the signal by the randn function, your noise follows a Gaussian distribution, translating into white noise for your signal. That is why your entire spectrum is showing sings of being slightly disturbed.

By the way, after fftshift, just get rid of the first half of your vector and you'll end up with positive frequencies.


sources for learning Image processing with MATLAB?

It does not update daily, but Steve Eddins' blog on image processing is an excellent source.

Steve is a co-author of Digital Image Processing using MATLAB, which is also a great source.

Another excellent source on image processing in general (not MATLAB-specific) are the lecture notes of Prof. Alan Peters from Vanderbuilt University.


Compute group delay of an audio file from STFT

The problem is that the stft function is splitting the signal up into different windows. That means that the signal from time $n$ to $n+N_{w}-1$ is multiplied by
$$ n, n+1, n+2, \ldots, n+N_w-1$$
instead of
$$
0, 1, 2, \ldots, N_w-1
$$
which is causing the scaling problem.

If I apply the group delay calculation from this derivation, I get:

where the top plot is the "data", the middle plot is the FFT of the data, and the bottom plot is the group delay. The group delay plot seems to be correct ($\tau_g = 49.5$, except for where the data is zero in the frequency domain.

To get this to work with the stft function, you're going to have to modulate the data with a sawtooth function rather than just np.arange(len(data)).

If I slice the sound file bachfugue.wav using the new, new code below I get the following image, which is more like what I'd expect.

New New Code Below

#!/opt/local/bin/python2.7
from scipy.signal import firwin
import numpy as np
import matplotlib.pyplot as plt
from scipy.io import wavfile

def group_delay(sig):
    b = np.fft.fft(sig)
    n_sig = np.multiply(sig, np.arange(len(sig)))
    br = np.fft.fft(n_sig)
    return np.divide(br, b + 0.01).real

rate, data = wavfile.read('bachfugue.wav')

N_w = 1024
L = np.floor(len(data)/N_w).astype(int)

full_delay = np.zeros((L, N_w))
for index in range(L):
    window = data[index:(index+N_w), 1].astype(float)
    full_delay[index][:] = group_delay(window)

plt.imshow(full_delay[0::10, 0::10], cmap=plt.get_cmap('hot'),aspect='auto')
plt.show()

New Code Below

#!/opt/local/bin/python2.7
from scipy.signal import firwin
import numpy as np
import matplotlib.pyplot as plt

N = 100
data = firwin(N, [0.05, 0.95], width=0.05, pass_zero=False)
n_data = np.multiply(data, np.arange(len(data)))
rate = 1

B = np.fft.fft(data)
Br = np.fft.fft(n_data)
group_delay2 = np.divide(Br, B + 0.01).real

plt.figure(1)
plt.subplot(311)
plt.plot(data)
plt.subplot(312)
plt.plot(np.abs(B))
plt.subplot(313)
plt.plot(group_delay2)

plt.show()

--

Original Answer

I've implemented it as I think it should be. I don't have your data, so I've just made a sinusoid.

The resulting plot is below. Code at end.

So I'm wondering what you should expect to see... more thought needed. I'll read that paper in more detail and see what I can make of it.

Now that I've had a chance to read through it, I'm wondering if that is correct.  Based on this derivation... I think the real way to do it might be different.  I'll try to code that up later.

Python Code

#!/opt/local/bin/python2.7
from scipy.signal import stft
import numpy as np
import matplotlib.pyplot as plt

N = 102400
rate = 1.0
data = []
for phi in ([2.0*3.14159265*x for x in range(0, N-1, 1)]):
    data.append(np.sin(phi))

n_data = np.multiply(data, np.arange(len(data)))
f, t, X = stft(data, rate, window="hamming", nperseg=800, return_onesided=False)
f_n, t_n, Y = stft(n_data, rate, window="hamming", nperseg=800, return_onesided=False)
# group_delay = (np.dot(X.real, Y.real) + np.dot(Y.imag, X.imag)) / np.power(np.abs(X), 2)

group_delay = np.divide(np.multiply(X.real, Y.real) + np.multiply(Y.imag, X.imag), np.power(np.abs(X), 2))

plt.imshow(group_delay)
plt.show()


16x interpolation using polyphase FIR - sum of the coefficients

Basically, you are up-sampling to a higher sample rate. The purpose of the interpolation filter is to eliminate the mirror-spectra, that you get from up-sampling and inserting zeros.

You need to chose the cutoff frequency for the interpolation filter accordingly. If you upsample by two you want a cutoff of fs/4 (or thereabouts). If you want to upsample by 8, your cutoff frequency needs to be lower, i.e. fs/16.


The following is an important paper about FFT algorithms:

Tran-Thong, “Algebraic Formulation of the Fast Fourier Transform,” IEEE Circuits and
Systems magazine, vol. 3, no. 2, June 1981, pp. 9-19.

In the above, the author distinguishes 16 basic types of FFT algorithm, characterized by: 1) decimation type (DIT - type T below, or DIF - type F below), 2) input order, 3) output order, and 4) geometry.  The author also shows graphs, but they will not be shown here.

Alg.    Input Order Output Order Geometry

T1      bit reverse sequential  in-place
T2      sequential  bit reverse in-place
T3      bit reverse bit reverse same output
T4      sequential  sequential  same output
T5      bit reverse bit reverse same input
T6      sequential  sequential  same input
T7      bit reverse sequential  isogeometric
T8      sequential  bit reverse isogeometric

F1      bit reverse sequential  in-place
F2      sequential  bit reverse in-place
F3      bit reverse bit reverse same output
F4      sequential  sequential  same output
F5      bit reverse bit reverse same input
F6      sequential  sequential  same input
F7      bit reverse sequential  isogeometric
F8      sequential  bit reverse isogeometric

Although many of the algorithms can be attributed to different authors, as stated in the above paper: "The Cooley-Tukey algorithm and its modifications are DIT algorithms.  The Gentleman-Sande algorithms and their modifications are DIF algorithms."

There are many more FFT algorithms other than the ones described above, including higher and mixed radix ones.  But they also exhibit similar characteristics.  For instance, take a look at the 6 point DIT radix-3/radix-2 graph shown on page 13 in the following reference:

http://www.ece.cmu.edu/~ee791/lectures/L18/FFTstructures.pdf

It is two radix-3 butterflies followed by three radix-2 butterflies.  It is bit-reversed input and sequential output.

Of course, given that the FFT can be viewed as a matrix transpose, one can also decompose a 6 point FFT into three radix-2 butterflies, followed by two radix-3 butterflies.  You can also have DIT or DIT (or mixed decimation).  So you have several options to choose from.


Audio Bit Depth "mapping" from DAW to Audio card: is this correct?

If your sound card has input range 2Vp-p (2 volts peak to peak = -1V to +1v) and has 24 bits of resolution (ENOB doesn't play any role in this, you still get 24 bits for each sample for each channel, but the higher ENOB the smaller noise you get), here is what voltage represents the output code of ADC:

00000000 00000000 00000000 (lowest possible result from ADC) => -1V
10000000 00000000 00000000 (midpoint) => 0V
11111111 11111111 11111111 (highest possible value) => +1V

Or with 3 bit ADC:

000 (0 dec) => -1V
001 (1 dec)
010 (2 dec)
011 (3 dec)
100 (4 dec) => (midpoint approximately) 0V
101 (5 dec)
110 (6 dec)
111 (7 dec) => +1V

So the voltage is proportionally distributed in the full ADC range, to calculate the voltage correspondig to specific code, you can use following formula:

U = (code/(2^bits-1)) * (Umax - Umin) + Umin

where Umin = -1V, Umax = 1V, bits = 24
For conversion between binary and decimal numbers refer to this document


Since the term linear does not appear in the question and the current answers, let me offer a complementary perspective.

A kernel in this acceptation (especially for images, which don't always follow linear rules, think about occlusion or saturation) is an array that is applied, somehow, on any input data. One often distinguishes linear and non-linear kernel (as one has linear and nonlinear filters, even is the terminology might seem improper).

Let us start from the linear point-of-view in the most specific sense: the filter array is applied as a convolution. Then @MarcusMuller's answer is perfect: a set, an array of linear filters, applied to input data as  convolutions to yield several separate output data.
Most  additional scalar linear operation (like the sum, the average, a weighted combination) on the output would be "useless": as they commute, summing the output is equivalent to summing the three filters in one single filter, and performing only one single convolution on the data.

Which leads us back to the objective in your comment; traditionally, a linear (analysis, I will come back on that later) filter-bank (FB) is used to split or separate data into components, often with separate spectra or a narrower ocontent (low-, mid- or high-frequencies for a three-band filter-bank). Or to merge different data streams into  others, with a wider spectrum. So a generic multi-input-multi-output (MIMO) FB takes one or several inputs, filters them into one or several outputs. One then distinguishes analysis or synthesis filter banks.

Generally, recombining the outputs from an analysis FB drifts away from the separation objective. But a single filter is also a filter-bank (not very interesting per se though).
But sometimes, this can be more efficient (computationally for instance).

Now, having narrower/wider outputs invites to rate variations, like downsampling and upsampling before or after the filters. To me, the most accepted sense of a filter bank is a bank of linear filters optionally combined with (linear, but not shift-invariant) upsampling or downsampling operations. And it is somewhat related to linear transforms, allowing expansion or shrinkage on the number of coefficients (they can be critical, oversampled or undersampled).

Then, people extend the notion to nonlinearity: filters can be nonlinear (like the median) and the kernels are interpreted as weights applied to a chunk of data. Or the data can be combined in non-linear ways, with $\min$, $\max$, AND or OR...

But in your case, as Marcus said, I'd bet on three standard filtered outputs. But in this case, there is not relationship between filters (except their kernel size), and what is powerful in the filter-bank theory is the linkage between the filters, and how one can optimize them. Now a couple of pointers:

Filter bank (wikipedia)
What is a Filter Bank?
Multirate filter bank and multidimensional directional filter banks


Convert complex valued sinusoid to real valued sinusoid

HINT:
It is based on the fact that
$$\sin(x) = \frac{e^{jx} - e^{-jx}}{2j}\quad\text{and}\quad j^2 = -1$$


This is most likely related to the comparator function and the waveform.

Long story short:

from scipy.signal import argrelextrema
from matplotlib import pyplot as plt
from scipy import linspace,pi,cos

#Time vector
t = linspace(0,4,4*1000); #4 seconds at Fs=1000Hz
#Phase vector
p = 2.0 * pi * t;
#Signal
S = cos(4.0*p); #A 4Hz sinusoid
# Have a look at it, focus on where the maxima and minima occur.
plt.plot(S);plt.show()
#Let's find some peaks
Pmin = argrelextrema(S, lambda x,y:x<y)
Pmax = argrelextrema(S, lambda x,y:x>y)

Pmin identifies all peaks: [500, 1500, 2499, 3499]

Pmax has skipped some    : [1000,2999]

Try it again, by going for Pmax2 = argrelextrema(S, lambda x,y:x>=y). Pmax2 now is [0,1000,1999,2000,2999,3999]

Now, peak detection is relatively straightforward. The mathematical way to think about it is to get the first derivative and look at where it drops to zero. Every zero of the first derivative is either a local maximum or local minimum. To decide if you are at a peak or a trough, you look at the value of the second derivative at that point. The second derivative of your signal, is the derivative of the first derivative. If that is positive, your candidate point is a local minimum, if it is negative, your candidate point is a local maximum.

That's fine, but costly and there are alternative ways to achieve it. For example, a peak could be defined as the pattern increase, decrease of signal values and a trough could be defined as the pattern decrease, increase of signal values. This is straightforward to understand but in reality, extrema might be "wide". In other words, the pattern of values now is increase, same, same, same, same, same, same, same, decrease. Where is the peak there? (or, conversely, the trough).

You might say, "That's not a problem, the peak is in the middle of the wide top" and that would be a reasonable assumption. In fact, you could achieve this by first "blurring" your data (taking the average) in such a way as to reduce the "wideness" and now you are back in a case of "increase, decrease".

This is achieved in argrelextrema by modifying the order parameter.

So, having said this, why does argrelextrema fails to identify those peaks in the extremely simple case of a cos?

If you take a close look at what the cos is doing, you will observe that:

When the peak is identified:

...,0.99993952,  0.99998889, 0.99999877,  0.99996914...

Did you spot the pattern? (increase,decrease)

When the peak is NOT identified:

..., 0.99995556,  0.99999506, 0.99999506,  0.99995556,...

Did you spot the pattern? (increase,same,decrease)

So, by modifying the comparator function, you enable argrelextrema to widen its search for what an extrema looks like.

Try introducing the <=, >= comparator function.

Hope this helps.


How does an image rotation affect the magnitude of a DFT?

Is there any way we could recover the angle the image was rotated by based on the magnitudes of the DFT?

Yes.
But some of these alternatives are very similar to the key idea behind the Hough transform.
It very much depends on how the "text image" is produced, but for a typical picture of a page of text (black letters on white background with margins), the "strong" coefficients will be vertical to the direction of the lines of text. The reason for this is not so much the text itself but the regular rectangles that are formed between the lines of letters. These appear as a regular square wave, the Discrete Fourier Transform (DFT) of which forms the "line of strong coefficients".
In general, after getting the spectrum, you would apply a threshold so that you retain only those strong peaks.
Without using the Hough Transform, the most straightforward thing to do here would be to fit a line to the points that make it through the threshold and get a vector that gives you the line's orientation.
The line fitting though would be sensitive to maxima that are not on the same line but still manage to go through the threshold. And this might distort the measurement.
This is very straightforward for the Hough Transform because it "converts" lines to points and due to its summation from different view points is much less sensitive to noise. Detecting line slopes with the Hough Transform is simply a matter of detecting local maxima in the Hough space.
But, this can be slow.
You can reduce the amount of work that the Hough Transform has to do by downsampling your original image in a controlled way. To understand this better, consider this:
Another alternative to the Hough Transform is to produce sums over a wide rectangle that rotates about the centre of the image and select the direction with the highest sum. This "rectangle counting" looks like this:

The gain here comes from summing over less orientations and summing over a band that has some width.
This is equivalent to downsampling the original image by a factor that is equal to the width of the band and reducing the range of angles that the Hough Transform has to produce sums over.
Finally, if we take this to the extreme and assuming that only 1 line is isolated from the thresholding operation, we could get away with just 2 "maxima scan lines".
These now look like this:

After thresholding the DFT spectrum, you are left with a handful of strong "peaks" that (ideally) form a line. Therefore, you could scan the perimeter of the thresholded spectrum for a local maximum. The position of that local maximum corresponds to the slope of the spectrum line.
If the maximum of the horizontal scan line is found in the middle, the line would be at 90 degrees (or, the page would be upright or upside down). If the maximum of the horizontal scan line is found in the beginning, the page would be rotated counter clockwise by 45 degrees. The same holds if the maximum occured at the end of the horizontal scan line which now denotes a 45 degree rotation clockwise.
The same applies for the vertical scan line now because there are some extreme rotations that would bring the line of coefficients entirely off the reach of that horizontal scan line. For those, we have the vertical scanline. And so on.
This is extremely fast because it does not rely on rotations at all, you just scan along the x-axis and y-axis, both of which are simple counters. But it does assume that there is only one strong "grid" in the image that will provide that line of DFT coefficients and it would still be sensitive to noise (or non-typical text images).
Ultimately, this comes down to balancing the accuracy of the assessment to the complexity of producing it.
Hope this helps.


Is there a hypothetical situation where once can use Wiener filter prior Matched filter to get better result?

A matched filter is matched to a noiseless known waveform.

A Weiner Filter is based on the correlation structure of a random waveform.

The signal models are different but often a signal has deterministic and random characteristics.

It would be a stretch to say that a Kalman Filter could reduce to a matched filter, but it does work in models without process noise.

The short answer is that one can mix detection and estimation together but the likelihood function is the better way to think about it than terms like “matched” or “Weiner” filter.


The number of FFTs is roughly N/(w-m). An fft costs you nfft*log2(nfft) so the total is proportional to N/(w-m)*nfft*log2(nfft).

A typical application would make nfft = w and m = w/2, so this would come out to be proportional to N*log2(nfft). So it's proportional to the length or you data set (obviously) and the log of the FFT size.


Matlab version rpmordermap and rpmfreqmap

At the very bottom of MATLAB manual pages there is always an info about the time when given function was first introduced. For both of rpm*map functions:


What does the term "railed" mean in signal processing?

A railed signal, or a railing signal, seems to indicate a flatline. On BIOPAC, Railing signal (flatline) says:

  When the amplified signal for any given channel exceeds the range -10
  to +10 volts, the signal will rail. You will see a straight line at
  -10 or + 10 volts (more likely the reading will be close 9.99 volts). The MP system is designed to work only in the range -10 to +10 volts.
  The signal could rail for several reasons (which are not exclusive)...

From Amplifiers: What do rail-to-rail and single supply mean?

  With respect to analog signals, a “rail” is a boundary that a signal
  has to work within.


can MFCC coefficients alone be used to do voice detection?

In essence, yes.

You would have to use something like Dynamic time warping to measure the similarities between two given signals, as long as they vary in time or speed.

From what I can remember, you can use DTW for speaking identification. BUT it is important to note that for single word recognition and single speaker identification this method would work but would not be suitable for commercial or large scale applications.

The DTW algorithm calculates the distances from the two signals, based on the Euclidean distances.. Whilst studying for my undergraduate degree, I did a simple word recognition system (In C++) and the machine learning algorithms (for instance, Hidden Markov Models) provided a much higher success rate than DTW.

Feel free to take a look at some of the code that I wrote for this project. Here

double Euclidean_distance(const vector<double> &actual, const vector<double> &training) {

    double distance = 0.0;

    for(unsigned i=0; (i < actual.size()); i++)
    {
        distance = pow((actual[i] - training[i]), 2);
    }

    return sqrt(distance);
}

double Distance(const std::vector< std::vector<double> > &actual, const std::vector< std::vector<double> > &training)
{
    int m = actual.size();
    int n = training.size();

    double cost[m][n];

    cost[0][0] = Euclidean_distance(actual[0], training[0]);

    for(int i = 1; i < m; i++)
            cost[i][0] = cost[i-1][0] + Euclidean_distance(actual[i], training[0]);
    for(int j = 1; j < n; j++)
            cost[0][j] = cost[0][j-1] + Euclidean::Euclidean_distance(actual[0], training[j]);

    for(int i = 1; i < m; i++)
            for(int j = 1; j < n; j++)
                cost[i][j] = std::min(cost[i-1][j], std::min(cost[i][j-1], cost[i-1][j-1]))
                    + Euclidean_distance(actual[i],training[j]);

    return cost[m-1][n-1];

}

(Note this is Academic code, not tested and just illustrate the algorithm)


Calculating the SNR from the power spectrum. Is my suggestion relevant?

Assuming the noise is evenly distributed among all frequencies, this approach makes sense in principle. However, some corrections:

You should not integrate over $pxx(n)f(n)dn$, since n is already the frequency. I think, you mean $\Delta_f(n)$, which is the frequency step between two bins. Assuming $\Delta_f(n)$ is equal and independent of n, then you can even cancel it in the equation.
You need to divide the noise and signal components by the respective bandwidth of the integral. Otherwise you will calculate more noise, when the bandwidth increases.
currently, you calculate (signal+noise)/noise. Though, you want to calculate signal/noise, i.e. (signal+noise-noise)/noise


Both indicators are used. It depends what you mean in speaking of your indicators.
One of common usages is to compare systems where we use indicators of the same type. For example,

to compare 5G NR and LTE RANs, we usually use max throughput which is defined as the multiplication of number of antennas, number of subcarriers and QAM order per subframe. There is no error probability here.
To compare channel decoders e.g. iterative belief propagation LDPC decoder and Ordered Statistic decoder, it does not make sense anymore if error probability is ignored and, therefore, the code dimension (number of bits before encoding) must be used.


Computer vision, Coupon code detection, sugested algorithms

You may try an OCR detector like Teseract with OpenCV:

https://github.com/tesseract-ocr

Train the Teseract algorithm with text samples obtained from the numbers you want to detect and it will work.

EDIT: If you want a similar implementation of what you want you may check this:

http://answers.opencv.org/question/10985/how-can-i-determine-the-location-of-lcdled-display-in-an-image-with-opencv-and-perl/


Change pitch of mp3 audio on a microcontroller

Use a "phase vocoder" which will use short term FFT properties to change the pitch without changing the duration.

There is plenty of information and links to "phase vocoder" on Wikipedia and Google in general.

It is a great educational exercise to go through its operation as you learn from it the critical importance of the phase information in the FFT for proper reconstruction.


Frequency domain filtering: must filter coefficients always be conjugate symetric?

I noticed that the coefficients, complex numbers, are not conjugate symetric.

Check what is the frequency associated with each coefficient, usually the coefficients are defined such that $a_i$ corresponds to the frequency $f_s i / N$, where $f_s$ is the sampling frequency, and $N$ is the number of coefficients. So what you have is $a_i = conj(a_{N-i})$, in some cases they may be circularly shifted, and so on.

the signal (audio_filtered) should again be real?

Yes, if you think in the time domain, the convolution involves only real values, so the result is a real signal. If you think in the frequency domain, both the filter and the signal have Hermitian symmetry, so the product of them will have hermitian symmetry as well.
Proof, let the signal with FFT $X$, and the filter with FFT $H$, the FFT of the filtered signal $Y$ is
$Y_{i} = X_i H_i = conj(X_{N-i})conj(H_{N-i}) = conj(X_{N-i}H_{N-i}) = Y_{N-i}

If the input filter is IIR, are things different with respect to conjugate symetry?

No, if you look at the fourier transform of an infinite signal
$$X(\omega) = \sum_{k=0}^{\infty} u_k exp(-2j\pi k \omega)$$
Assume $u_k$ is real, and will will notice that since $e^{-j x}$ is the complex conjugate of $e^{+jx}$, each term will be conjugated and so, $X(\omega)$ is the complex conjugate of $X(-\omega)$


How Many Samples Are Identical In Linear and Circular Convolution Result?

assuming $y[n]$ is padded with $L-P$ zeros (to make it the same length as $x[n]$ and your circular convolution length), then the number uncorrupted samples is $L-P+1$.

you might want to Google "Overlap Save" (OLS, a better name i have read is "Overlap-Scrap"), which is directly about this issue of uncorrupted samples in circular convolution.  the alternative is "Overlap Add" (OLA) which has no corrupted samples, but you have to overlap the tails of each block and add them.

this is well covered in books like O&S.


The Least Norm Solution of Under Determined Linear System

Least Squares solution is always well defined for Linear System of Equations.
In your case, which is under determined it means there are many solutions to the Linear Equations.
The Least Squares solution has nice property, it also minimizes the $ {L}_{2} $ norm of the solution (Least Norm Solution) hence it is well defined.
In practice, in order to solve this equation you need to use the pinv() function in MATLAB.
This function uses the Moore Penrose Pseudo Inverse Matrix of $ A $ to solve the equation:
x = pinv(A) * b

It basically decompose the Matrix using the Singular Value Decomposition (SVD) such that $ U S {V}^{T} = A $ and the rebuild the matrix $ {A}^{\dagger} = U {S}^{\dagger} {V}^{T} $ where $ {S}^{\dagger} $ is the result of inverting any non zero element of the diagonal matrix $ S $.
Update
The user @percusse said the following:

You just need A\b
You shouldn't use pinv for linear system solutions that's what I mean. Not only inefficient but also not accurate. There is no unique solution hence you can't talk about a least.
Backslash operator is optimized for all possible linear eq sets. Try it.
There is no least squares solutions for underdetermined systems because of the nontrivial nullspace. Anyways nevermind.

So, to make things clear.
The Least Squares solution always exists and it is also the Least Norm solution.
Namely, from the set $ \mathcal{B} = \left\{ x \mid \min_{x} \left\| A x - b \right\| \right\} $ there is one unique $ \hat{x} $ such that $ \hat{x} = \arg \min_{x} \left\| x \right\|_{2} \text{  subject to  } \left\{ x \mid x \in \mathcal{B} \right\} $.
For instance:
for ii = 1:1e12
    numRows = 3;
    numCols = 27; %<! Must be bigger than 'numRows' for Under Determined System

    mA   = randn([numRows, numCols]);
    vB   = randn([numRows, 1]);
    vXLs = pinv(mA) * vB;
    vX   = mA \ vB;

    if(norm(vXLs) >= norm(vX))
        keyboard();
    end

    if(mod(ii, 1e4) == 0)
        disp(['Iteration Number - ', num2str(ii)]);
    end

end

Remarks:

The difference arises for Under Determined System. For over / well determined system the solutions collide (Namely Tall Matrices which are the "Classic" case in LS problems).
The Psuedo Inverse is much slower to use then the \ operator. If you're looking for any solution, use \. If you're looking for the Least Squares solution (The Least Squares solution is always the least norm solution) for Under Determined System you need to use the pinv() function.
Proof and deeper analysis can be found in my Singular Value Decomposition (SVD) Presentation. Specifically this issue is on pages 50-60.
Cleve Moler's (Inventor of MATLAB) post on Pseudo Inverse vs. \ - The World’s Simplest Impossible Problem.


Looking for an algorithm that computes "Similarity" between two signals

Coherence? As in: http://en.wikipedia.org/wiki/Coherence_(signal_processing). By comparing the signal spectrums, you get real number numbers between 0 and 1, which each number is a measure of the similarity of the two signals at each frequency.  But it might be difficult for a beginner to understand and program.  You might look for simpler explanations and/or examples by doing a Google search for 'signal coherence', which results in about 17 million hits.


What is the difference between PSD and squared magnitude of frequency spectrum?

The power spectral density describes the density of power in a stationary random process $X(t)$ per unit of frequency. By the Wiener-Khinchin theorem, it can be calculated as follows for a wide-sense stationary random process:

$$
S_{xx}(f) = \int_{-\infty}^{\infty} r_{xx}(\tau) e^{-j2\pi f \tau} d\tau
$$

where $r_{xx}(\tau)$ is the autocorrelation function of the process $X(t)$:

$$
r_{xx}(\tau) = \mathbb{E}\left(X(t)X(t - \tau)\right)
$$

This is only valid for a wide-sense stationary process because its autocorrelation function is only a function of the time lag $\tau$ and not the absolute time $t$; stated differently, this means that its second-order statistics don't change as a function of time.

With that said, if you have a sufficiently-detailed and accurate statistical model for your signal, then you can calculate its power spectral density using the relationship above. As an example, this can be used to calculate the power spectral density of communications signals, given the statistics of the information symbols carried by the signal and any pulse shaping employed during transmission.

In most practical situations, this level of information is not available, however, and one must resort to estimating a given signal's power spectral density. One very straightforward approach is to take the squared magnitude of its Fourier transform (or, perhaps, the squared magnitude of several short-time Fourier transforms and average them) as the estimate of the PSD. However, assuming that the signal you're observing contains some stochastic component (which is often the case), this is again just an estimate of what the true underlying PSD is based upon a single realization (i.e. a single observation) of the random process. Whether the power spectrum that you calculate bears any meaningful resemblance to the actual PSD of the process is situation-dependent.

As this previous post notes, there are many methods for PSD estimation; which is most suitable depends upon the character of the random process, any a priori information that you might have, and what features of the signal you're most interested in.


What does the $n$ do in $1/(j \omega)^n$

What do you mean with Beam in this context.... Particle Beams in LINACs? Structural Beams in Buildings?.

$j\omega$ is a phasor term from frequency analysis.

$\frac{1}{j\omega}$ is the integration as @msm showed.

If you want to integrate two times, just use cumtrapz two times.

But as you said, this is only a "term" complementing a frequency model, which it would be really good to have here in order to help you more... You never simulate that term alone.


Is it necessary to downsample data for filtering of low frequencies?

Yes, you can get much better results if you downsample your data first. The problem occurs because the width of your transition band relative to the original data sample rate is very small. This can cause numerical problems when designing or implementing digital filters. The problem only gets worse as the filter order increases, especially with IIR filters.

With that said, in some cases, you might have more success (i.e. better numerical stability) if you convert the filter to second-order sections. However, for this case, I think the right answer is to downsample the data first.

The problem likely occurs because of the very small cutoff frequency on the low end of your bandpass filter. There is very little room for the transition band (less than 0.1 Hz, which is 1/320000-th of the sample rate), which is difficult to implement. If you downsample the data as much as possible, then you can at least increase this ratio by a good amount, which should give you better results.

I would also try using FIR filter implementations if possible, as they typically have superior numerical stability due to their lack of feedback.


Low data rate with high frequency offset

The answer is yes it is possible and any amount of frequency offset can be corrected to the extent your sampling rate is sufficiently high enough for unambiguous frequency resolution over ±10 KHz. A very simple frequency discriminator that can be used in a carrier tracking loop is further described at this post:Carrier frequency offset estimation for QPSK and asymmetrical spectrum and results in determining the frequency offset simply from the product of two successive (complex) samples given as $I_1Q_2-I_2Q_1$. For BPSK you can also remove the modulation from the frequency error result by first squaring the signal, in which case the result will be twice the actual frequency offset; which is of no consequence since the loop will drive the offset to 0 in any event).
This discriminator has unambiguous frequency resolution for IQ rotation up to 90° between samples. A 10 KHz frequency offset will rotate 90° if the sampling rate is 40 KHz, so as long as your sampling rate can be greater than this rate, this should be a very simple approach to carrier recovery / Doppler wipe-off (80 KHz if using the frequency doubling approach to strip the modulation).
The maximum possible offset sets your search window and the header duration if used for acquisition sets the carrier recovery Loop BW. The carrier Loop BW sets the acquisition SNR required (wider BW is faster but requires a higher SNR). So use a loop BW that is wide enough to settle for a step response in frequency change from 0 to 10 KHz within the time duration of 50 symbols. (And then the loop bandwidth can be further tightened once acquired if channel dynamics allow).


To better extend the signal at the border, you can consider a constant extension: $[x_0,x_1,x_2]$ can be extended to the left (and similarly to the right) as $[x_0,x_0,x_0,x_1,x_2]$. However, this tends to break the slope, as it assumes a zero-derivative (flat signal).

A more advanced option resides in anti-symmetric extensions: they are sometimes called whole-sample or half-sample extensions, the two main options are:

$$[2x_0-x_2,2x_0-x_1,x_0,x_1,x_2]$$
and
$$[2x_0-x_2,2x_0-x_1,x_0,x_0,x_1,x_2]$$
In other words, the average of the extreme points is  $x_0$. This better respects derivatives in low noise conditions. Other types of extrapolation exist, more ad hoc (linear or parabolic fits), but may help. In the Matlab wavelet toolbox, Border Effects are described, and some Discrete wavelet transform extension modes are:

'asym' or 'asymh': Antisymmetric-padding (half-point): boundary value antisymmetric replication
'asymw': Antisymmetric-padding (whole-point): boundary value antisymmetric replication
'spd' or 'sp1': Smooth-padding of order 1 (first derivative interpolation at the edges)
'sp0': Smooth-padding of order 0 (constant extension at the edges)

As said by @DaBler, such extrapolations may require an increase in the transformed coefficients but:

for mere processing, you can crop the extended parts afterward,
there exists symmetric/antisymmetric filter-banks (e.g. $M$-band wavelets, $M >3$ and even) for which the extension can be made non-expansive.

Aside, if you really care about signals close to the border, you should take care of the sub-sampling, and possibly avoid it, as you cannot afford loosing samples. And more importantly, make sure that your wavelet is (close to) symmetric as well, unless you can lose the positive effect of your extension.

Finally, wavelets build on the interval adapt the wavelet size close to the border.


How is the Laplacian used for Unsharp Mask?

The regular sharpening i.e. high pass filter, essentially applies the inverse filter of blurring.

Basically, when you do unsharp masking - instead of subtracting the blurred image, you add the edge image as a portion of it. we can have $g(x,y)$

$$g(x,y)=f(x,y)-f_{smooth}(x,y)$$

which is nothing but an edge image using LoG.

The sharp image is produced by adding portion of this is added back to original image and hence contrast near edge increases drastically. This is how you produce the sharp image $f_{sharp}$

$$f_{sharp}(x,y)=f(x,y)+k*g(x,y)$$

Courtesy http://homepages.inf.ed.ac.uk/rbf/HIPR2/unsharp.htm.

More details are given in that page.


How do I reconstruct a covered-over part of a regular pattern

You can try using the Hough transform to detect the disk and the visible ellipses.

Also, you may want to do a search for "texture generation".

Most importantly, keep in mind that just because a person can do something effortlessly, it doesn't mean that it is easy to program a computer to do it. Virtually any vision task is very easy for people, but are extremely hard for computers.


convolution of a real and complex signal

You are confusing several different notions. First of all, two signals $x(t)$ and $y(t)$ are orthogonal if

$$\int_{-\infty}^{\infty}x(t)y^*(t)dt=0\tag{1}$$

From (1) it is clear that if two real-valued functions $x(t)$ and $y(t)$ are not orthogonal, i.e. the integral in (1) is not equal to zero, then also $x(t)$ and $jy(t)$ are not orthogonal. So your argument that the two given signals must be orthogonal because one of them is real-valued and the other one is imaginary is wrong.

Furthermore, even if two signals are orthogonal to each other, their convolution is usually non-zero. Just take two rectangular functions that do not overlap. Obviously, they are orthogonal, but their convolution is a triangular function.

Consequently, one way to interpret the convolution of the two sinc functions is that a low pass signal with cut-off frequency $W$ is filtered by a low pass filter with cut-off frequency $2W$, and multiplied by the (irrelevant) factor $-j$. Since the low pass filter has no effect on the low pass signal, the only change to the signal is the constant factor $-j$.


Distortion in phase difference measurement using FFT

The OP mentions "demodulating" but I don't actually see that in the process (meaning the phase modulation is not actually removed).  Without demodulation the  signal can possibly be spread across multiple DFT bins thus reducing the SNR in the bin that was used to measure the phase. If the signal bandwidth is less than one bin, and the noise is still white, then there would be further increase in SNR by the rate of $\sqrt{N}$ by increasing the number of samples by $N$.
The extent of this spreading due to modulation isn't clear since only the modulation index was provided; we would also need to know the modulation rate and type of waveform in order to predict the occupied bandwidth. With both the modulation index and rate with a continuous modulation waveform the occupied bandwidth can be approximated using Carson's Rule which is applicable to FM and PM given as:
$$BW \approx 2(\Delta F + F_m) \tag{1} \label{1}$$
With:
$BW$: Approximate occupied BW
$\Delta F$: frequency deviation $F_m$: the modulation rate.
And using the modulation index $\beta$ which is the peak phase angle (for sinusoidal PM) and related to $\Delta F$ and $F_m$ according to:
$$\beta = \frac{\Delta F}{F_m} \tag{2} \label{2}$$
From \ref{1} and \ref{2} we can establish the occupied bandwidth from the modulation rate and modulation index directly as follows:
$$BW \approx 2F_m(\beta + 1) \tag{3} \label{3}$$
From this we can confirm if the occupied bandwidth exceeds the resolution bandwidth of the measurement, at which point the SNR from a single bin will be maximized as limited by the DFT noise floor. If the occupied bandwidth is less, and the DFT noise floor still white in close proximity to the signal (meaning the signal is sufficiently stationary over that duration of time), then there will be a further improvement simply by increasing the number of samples.
The resolution bandwidth for each bin in the DFT given in Hz is given by the following formula as the "Equivalent Noise BW":
$$\text{ENBW} = \frac{f_s}{N}w_b \tag{4} \label{4}$$
With:
$\text{ENBW}$: Equivalent Noise Bandwidth (Hz)
$f_s$: Sampling rate (Hz)
$N$: Number of samples in DFT
$w_b$: ENBW of the window function in units of DFT bins
Computing the ENBW of window functions is detailed here: Find the Equivalent Noise Bandwidth.  The equivalent noise BW of the Blackman Harris window is approximately 1.7 bins, so the equivalent noise bandwidth in Hz with the 10 MHz sampling rate and 32768 bins would be:
$$\text{ENBW} = \frac{10e6}{32768}1.7 = 518 \text{ Hz}$$
Using this result in \ref{3} results in the maximum modulation rate of
$$F_m \le \frac{\text{ENBW}}{2(\beta+1)} = 123.3 \text{ Hz}$$
Therefore, if the modulation rate is not significantly higher than 123 Hz, the signal will be entirely within the resolution bandwidth of the DFT and no further phase demodulation would be needed, within the SNR provided by the DFT as was done. If the modulation rate is significantly less, and the signal still stationary over longer duration (eventually phase noise of the signal source will limit this condition), then SNR can be further improved by $\sqrt{N}$ by increasing the number of samples by $N$, reducing the DFT noise floor. If the modulation rate is significantly more than 123 Hz, then there would be additional benefit in SNR to first demodulating the signal to remove the phase modulation. (Also to note that at higher modulation indexes, more signal energy is concentrated in the outer portions of the bandwidth such that there would be a decrease in SNR with further time duration or increasing $N$). When the modulation rate is less than the ENBW of the DFT, even though the signal is phase modulated, given the "courseness" of the DFT, the signal would appears just as a pure tone in this case since the DFT is providing us the average frequency, average phase over the duration of the capture.
So continuing the assumption of $F_m \le 123 \text{ Hz}$, with a phase variations of 0.15° and no motion (meaning this is phase noise); assuming the 0.15° was an rms measurement, that would be equivalent to an SNR as limited by noise due to phase of:
$$-20log_{10}((0.15) 2\pi/360) = 51.6 \text{ dB} $$
This is a 51.6 dB SNR in the resolution bandwidth of 518 Hz; if the signal for which the DFT was computed was real, then the SNR of the signal would be 3 dB higher given it would be at both bins 327 and 32441 and only one of those two bins was used for the measurement. The OP has not stated what bandwidth was used for the original 30 dB SNR metric of the signal, but this approach as outlined can be used to further quantify if the process is noisier than expected, or not. If the variation in phase was anything else such as peak to peak or peak measurements would translate to an even higher SNR, but that should be clarified as well for an accurate assessment. This should also help quantify the SNR requirements for a certain phase accuracy as needed for the final distance measurement.

When such phase noise is small, the small angle approximation $\sin(\theta) \approx \theta$ applies, thus the rms of any phase noise when normalized to the signal (typically given as "dBc") is the SNR when both units are squared to be power terms (hence $20Log_{10}(\theta_{rms}$))


well, there is always autocorrelation $$ R_x(\tau)=\sum x[n] x[n+\tau] $$ or AMDF $$ Q_x(\tau) = \sum |x[n] - x[n+\tau]| $$ or ASDF $$Q_x(\tau) = \sum (x[n] - x[n+\tau])^2 $$  with the latter there is the relationship between autocorrelation and ASDF $$ R_x(\tau)=R_x(0) - \frac12 Q_x(\tau) $$.  a measure of periodicity might be $\frac{R_x(P)}{R_x(0)} $ where $R_x(P)>R_x(\tau)$ for all $\tau\ne 0$.  sometimes we call "$P$" the "period".  if $\frac{R_x(P)}{R_x(0)} \approx 1 $ we say it's pretty periodic.


How to programmatically find a sharp change in slope

This looks like an interesting problem.

First, if your timing is evenly spaced, your idea of having two "pin" (aka calibration) points should give you an adequate linear relation between your two time scales.  If you have more than two points then you can use linear regression to find a best fit.

Second, I think you will want to use the raw values and values of a calculated first derivative to classify parts of your graph.  The actual value of a second derivative doesn't seem to me to be particularly helpful.  Because of the noise you have, using immediate first difference calculations may lead to bad results.

I have a suggestion for a different approach which may work quite well.  Check out the difference technique in my blog article Exponential Smoothing with a Wrinkle.  You need not necessarily use exponential smoothing, boxcar smoothing or some other would also work.  The idea is to smooth in a forward direction and also smooth in a backward direction.  At every point, the forward smoothing gives you the "expected value" from the past and the backward smoothing gives you the "expected value" from the future.  Substracting the former from the latter gives you a first derivative like signal.  The coding is rather trivial so I will not provide it.  I suggest you give it a shot on your data and see what the graphs look like.  The smoothing will solve your noise issue.

I'm not sure the best way to calculate your probability values, but I think with a few graphs and calculations under your belt you will be able to provide an estimate that will be valid and fairly consistent across all your cases.


Is blending of colour lights a linear system?

Your mistake is that you are describing subtractive color mixing, not applicable in this scenario.

If you shed two lights onto a diffuse surface that does not absorb any light, and collect all the light diffusely reflected, the sum will be an additive mixture of the two. Incoherent (non-laser) light can be thought of as noise that has some electromagnetic spectral power distribution. From the signal processing point of view, long-time average power spectra of independent noise signals mix additively when they are summed