Are both halves of a symmetrical FFT result required to generate the original waveform? Performing an FFT on a purely real signal produces a symmetrical result. I want to recreate the $50\textrm{ Hz}$ mains hum to actively subtract it from my signal, but:

are both halves of the FFT result required to do this?
Is there any drawback to only using one half of a symmetrical FFT result?

You can perfectly use, $X_h[k]$, one half of the conjugate symmetric DFT $X_f[k]$ which belongs to a real signal $x[n]$ of length $N$.

Of course in doing so, you must symmetrically pad the half length DFT $X_h[k]$ to the full length DFT $X_f[k]$, which is obviously necessary for an unaliased reconstruction of $x[n]$ from its proper length DFT $X_f[k]$ of length N.

Note that during the creation of the mirrored padded full DFT $X_f[k]$ from the symetric half $X_h[k]$, the first sample $X_h[0]$ is not repeated and the remainig samples $X_h[1]$ to $X_h[M-1]$ are mirrored with their phase being negated (i.e. they are conjugated), where $M$ is the length of the half length DFT $X_h[k]$.

The length $M$ of the half sized DFT $X_h[k]$ is given by the following: If the signal's (or its natural DFT $X_f[k]$'s) length $N$ is even then $M= 1 + \frac N 2 = \frac{N+2}{2}$, otherwise, if $N$ is odd $M = 1 + \frac{N-1}{2} = \frac{N+1}{2} $.

Also when $N$ is even, the last sample of the half sized DFT, $X_h[M]$ is not repeated, and when $N$ is odd then the last sample $X_h[M]$ is conjugate mirrored in formation of the full length DFT $X_f[k]$, as examplified below Matlab/Octave code:

x = [1:8];      % N = 8 , a signal of length N = 8 (even)
X = fft(x,8);   % X[k] is the natural DFT of x[n]
Xh = X(1:5);    % M = (8+2)/2 = 5, also matlab indexing starts from k=1
Xf = [ Xh conj(Xh(end-1:-1:2))];    %N:even => Xh[0] and Xh[5] not repeated
xr = ifft(Xf,8); % reconstruct the original signal from Xf[k]
xr = 1     2     3     4     5     6     7     8

as expected.


How can a linear operator on DFT vector produce the same vector using only half of the DFT vector? Suppose there is a DFT vector $\mathbf{X}$ with length N, which presents complex conjugate symmetry around its middle point, i.e., $X(1) = X(N-1)^*$, $X(2) = X(N - 2)^*$ and so forth.  $X(0)$ and $X(N/2)$ are the DC and Nyquist frequency respectively, therefore are real numbers. The remaining elements are complex.

Now, suppose there is a matrix $\mathbf{T}$, with size $N \times N$, which multiplies vector X.

\begin{align}
\mathbf{Y} = \mathbf{T}\mathbf{X}
\end{align}

Assuming the operator $\mathbf{T}$ preserves the complex conjugate symmetry in $\mathbf{Y}$ (see Conditions for precoding matrix to preserve complex conjugate symmetry on DFT vector for more details), how can I obtain the same "positive frequencies" (from $X(0)$ to $X(N/2)$), using an $(N/2 + 1) \times (N/2 +1)$ matrix $\mathbf{T}$, instead of an $N \times N$?

The motivation for this question is to reduce complexity in the operation.

Consider that $X$ is a DMT symbol and its elements are mapped from QAM constellations with different sizes. Consider also that only the positive frequencies carry useful information, the "negative" frequencies are only used for the IFFT to be real.

EDIT:
Is there a low-complexity solution other than using the $N/2+1$ upper rows of $\mathbf{T}$ multiplying the length-$N$ vector $\mathbf{X}$? This reduces the complexity to $N(N/2+1)$ CMACS plus the operations required to assemble the hermitian symmetric symbol

EDIT(2):
I suggest the following code to illustrate the problem:

N = 8;
A = rand(N,N); %must be real-valued
w = exp(-1j*2*pi/N); % twiddle factor
W = w.^(repmat(0:N-1,N,1).*repmat(0:N-1,N,1).'); % DFT matrix
T = W*A*W' % Linear operator
x =  fft(rand(8,1));

It is obviously clear that:
   T(1:N/2+1,:)x
Returns the positive frequencies from:
   Tx
However,

Using x(1:N/2+1), what linear operator (with size $(N/2 + 1) \times (N/2 + 1)$) should I use to obtain the positive frequencies?

I've suggested in the comments above (and I think Matt L is trying to say the same thing also) that you should be able to just discard the negative frequency elements entirely for your analysis. They don't seem to be relevant to what you're doing (which makes sense, since they don't carry any information, as the time-domain signal is real). Look at your equation again:

$$
\mathbf{y} = \mathbf{Tx}
$$

(I changed the casing so that lower-case symbols are vectors and upper-case are matrices, as is the typical convention)

Assume for the time being that $\mathbf{x}$ and $\mathbf{y}$ are length-$N$ and thus $\mathbf{T}$ is $N$-by-$N$. This would correspond to the general case, where the time-domain signals corresponding to $\mathbf{x}$ and $\mathbf{y}$ need not be real. If $\mathbf{y}$ is a column vector, then you can think of the matrix operation instead as:

$$
\mathbf{y} = \sum_{k=1}^N x_k\mathbf{t_k}
$$

where $x_k$ is the $k$-th element of the vector $\mathbf{x}$ and $\mathbf{t_k}$ is the $k$-th column of the matrix $\mathbf{T}$. My question at this point would be, instead of using the full $N$-length vector and $N$-by-$N$ matrix, why not just use the $N/2$ rows that contain information? Then, the calculation would turn into:

$$
\mathbf{y^+} = \sum_{k=1}^{N/2} x_k\mathbf{t_k}
$$

(neglecting the details of the boundary cases for the zeroth and Nyquist bins, which depend upon whether $N$ is even or odd)

The resulting vector $\mathbf{y^+}$ contains the positive frequencies for the desired result. You can then just symmetrically extend it to get a full length-$N$ vector that has the conjugate symmetry that you seek. This should make intuitive sense in the following ways:

If there is no information in the negative frequencies in the input vector $\mathbf{x}$ (which there isn't, since the entire vector is defined by its positive frequency components), then you shouldn't need to use them in your calculations.
If there is no information in the negative frequencies in the output vector $\mathbf{y}$ (which there isn't, since the entire vector is defined by its positive frequency components), then you shouldn't need to waste time calculating what their values would be.

This should free you from worrying about any bizarre constraints on the properties of the matrix $\mathbf{T}$. I think you're trying to think of the action of the precoder as too coupled with the OFDM modulation. Instead, think of it just as a preprocessing operation that happens on the encoded symbols before you assign them to frequency bins in the OFDM modulator's IDFT.


FMCW Radar Signal Processing Flowchart / Ambiguity Function I Read many papers on signal processing for pulsed and FMCW radars, taking the FMCW type, the signal processing start after mixing the received and transmitted signal, the result of this operation is a signal called FMCW beat signal. We apply 2D-FFT to extract targets from the beat signal then we apply a detection algorithm to recognize true targets and reject others.
I read also from some papers that they use a tool called ambiguity function, I don't know where this function is situated in fmcw signal processing flowchart or this function is just with pulsed radar!
My question is :
What's the role of ambiguity function in radar signal processing? and when we should use it?

The Ambiguity Function is just a name for 2D Correlator.

It is known that given a shifted 1D Signal the optimal estimator for the shift (Range in RADAR / LIDAR, etc...) it the correlation with the signal.

What if our signal had 2 parameters and both are shifted?
In the case of RADAR we have a shift of the signal due to the range (Time Shift) and due to velocity (Frequency Shift - Doppler Effect).

So how would you solve this?
Well, correlation with the transmitted signal shifted on both parameters.
The result of this operation is the Ambiguity Function which is a grid (Discrete) saying the correlation value between the received signal and a shifted (In both parameters) transmitted signal.


Common Use Cases for 2D Non Separable Convolution Filters In the image processing world, I've noticed that a lot of the popular convolution filters are separable. Here's a quick list of common separable filters:

Sobel
Gaussian blur
Box filter (all ones, for example)

Not separable:

Laplacian
What other common filters are not separable?

I'm a little bit worried about this post being closed for being too broad or not constructive enough. But, some "what is this good for?" questions have been well-received here.

Background: I've developed some high-performance 2D convolution implementations. I'd like to have a better handle on the motivating applications.

There are many.  In general, steerable filters are not separable.  Few are, but many are not.

Here is a short list:

Gabor filters in certain directions.  Of course, one can use the bases that are separable due to the fact that Gabor filters are steerable.  Gabor filters have the advantage that they are optimally localized in time and frequency.  As a degenerate case, consider a non-symmetric 2D Gaussian where the principal components are not horizontal or vertical.
Difference of gaussian (DoG) used to model retinal ganglion and thalamic receptive fields are also not separable.
Marr-Hildreth for edge detection is not separable.
A bunch of filters used for segmentation are not separable: http://dl.acm.org/citation.cfm?id=720464&CFID=260338827&CFTOKEN=99199893

The list can go on and on, and these are just common examples!  You can certainly design FIR filters to fit your need, and there are methods to do this.  That is the whole purpose of a large sub-field of DSP.

In general when you design a 2-D filter, if it is anisotropic and not oriented horizontally or vertically, chances are that it is not separable.  For example, if you were to design a low pass (pre-)filter to remove high frequency that may cause aliasing in a hexagonal sampling system.  This would have to be an analog filter, but for digital, consider a down sampling scenario.  For optimal design, you'd want to keep as much information as possible, so your filter would most likely match the hexagonal sampling grid.  This will definitely give you a filter that is not separable.


Energy in baseband and bandpass signal Is there any difference in energy content of baseband and bandpass signals? When we see their representation in frequency domain, it looks like a bandpass signal has double the energy content compared to basband. But my friend says there's no difference. Someone please clarify.

Let $x(t)$ denote a real-valued finite-energy signal of energy $E_x$, that is,

$$\int_{-\infty}^\infty |x(t)|^2\,\mathrm dt = E_x < \infty \tag{1}$$

and let $X(f)$ denote the Fourier transform of $x(t)$.
Parseval's theorem assures
us that

$$\int_{-\infty}^\infty |x(t)|^2\,\mathrm dt = \int_{-\infty}^\infty |X(f)|^2\,\mathrm df \tag{2}$$

and thus the energy $E_x$ can also be computed via the right side of $(2)$.

Now consider the DSB-AM (Double-Side-Band Amplitude-Modulated)
signal $x(t)\cos(2\pi f_c t)$. Its energy is

$$\begin{align}
\int_{-\infty}^\infty |x(t)\cos(2\pi f_c t)|^2\,\mathrm dt
&= \int_{-\infty}^\infty |x(t)|^2|\cos(2\pi f_c t)|^2\,\mathrm dt\\
&= \frac 12 \int_{-\infty}^\infty |x(t)|^2[1 + \cos(2\pi 2f_c t)]\,\mathrm dt\\
&= \frac 12 E_x + \int_{-\infty}^\infty |x(t)|^2\cos(2\pi 2f_c t)\,\mathrm dt\tag{3}
\end{align}$$
Now, if the value of the integral on the right side of $(3)$ is exactly zero,
then the modulated signal $x(t)\cos(2\pi f_c t)$
has exactly half the energy of the unmodulated signal $x(t)$.
In the time domain, it can be plausibly argued that if the integral is
not identically zero, it is at least very small when $x(t)$ is a
low-pass signal and $f_c$ is very large compared to the highest frequency
in $x(t)$.  Denote by $T = (2f_c)^{-1}$ the period of the double-frequency
sinusoid in $(3)$ and consider that
$$\int_{-\infty}^\infty |x(t)|^2\cos(2\pi 2f_c t)\,\mathrm dt
= \sum_{n=-\infty}^\infty \int_{nT}^{(n+1)T}|x(t)|^2\cos(2\pi 2f_c t)\,\mathrm dt.
\tag{4}$$
Now, $x(t)$ is a slowly-varying signal that does not change very much
during $[nT, (n+1)T]$ and so
$$\int_{nT}^{(n+1)T}|x(t)|^2\cos(2\pi 2f_c t)\,\mathrm dt
\approx |x(nT)|^2 \int_{nT}^{(n+1)T}\cos(2\pi 2f_c t)\,\mathrm dt \approx 0.$$
Note also that
$\displaystyle \int \cos(2\pi 2f_c t)\,\mathrm dt
= -\frac{1}{4\pi f_c}\sin(2\pi 2f_c t)$ so that any contribution from
$\int_{nT}^{(n+1)T}|x(t)|^2\cos(2\pi 2f_c t)\,\mathrm dt$, which we have
already argued is going to be small is further tamped down towards $0$ by
that $\frac{1}{4\pi f_c} << 1$ factor.

Calculations such as these are more easily handled in the frequency
domain. We have the Fourier transform pair
$$x(t)\cos(2\pi f_c t) \leftrightarrow
\left.\left.\frac 12\right[X(f+f_c)+X(f-f_c)\right]$$
so that Parseval's theorem gives us that
$$\begin{align}
\int_{-\infty}^\infty |x(t)\cos(2\pi f_c t)|^2\,\mathrm dt
&= \frac 14 \int_{-\infty}^\infty |X(f+f_c)+X(f-f_c)|^2\,\mathrm df\\
&= \frac 14 \left[\int_{-\infty}^\infty |X(f+f_c)|^2\,\mathrm df
+ \int_{-\infty}^\infty |X(f-f_c)|^2\,\mathrm df
\right.\\
&= \qquad \left. + \int_{-\infty}^\infty X(f+f_c)X^*(f-f_c)
+X^*(f+f_c)X(f-f_c)\,\mathrm df\right]\\
&= \frac 12 E_x + \frac 14 \left[\int_{-\infty}^\infty X(f+f_c)X^*(f-f_c)
+X^*(f+f_c)X(f-f_c)\,\mathrm df\right].\tag{5}
\end{align}$$
We now argue (as in the comments) that if the support
of $X(f)$ is $[-W,W]$ and $f_c > W$, then the supports of
$X(f+f_c)$ and $X(f-f_c)$ are $[-f_c-W, -f_c+W]$ and
$[f_c-W, f_c+W]$ respectively, and since $-f_c+W < 0 < f_c-W$,
the supports do not overlap so that the integrand in $(5)$
has value $0$ for all $f$. Thus, $f_c > W$ is an explicit sufficient
condition for the energy in $x(t)\cos(2\pi f_c t)$ to be half of
the energy in $x(t)$.  The modulation process "loses" half the
energy because $\cos(2\pi f_c t)$ is a carrier signal of power $\frac 12$;
if we had used a unit-power signal such as $\sqrt{2}\cos(2\pi f_c t)$,
there would be no "loss".


Python audio analysis: which spectrogram should I use and why? I am doing my final project at university: pitch estimation from song recording using convolutional neural network (CNN). I want to retrieve pitches existed in a song recording. For CNN input, I am using a spectrogram.
I am using MIR-QBSH dataset with pitch vectors as data label. Before processing the audio to CNN (each audio has 8 sec duration in .wav files of 8 KHz, 8 bit, mono), I need to pre-process the audio into a spectrogram representation.
I have found 3 ways to generate a spectrogram, the code are listed below.
Audio example I am using in this code is available here.
Imports:
import librosa
import numpy as np
import matplotlib.pyplot as plt
import librosa.display
from numpy.fft import *
import math
import wave
import struct
from scipy.io import wavfile

Spectrogram A
x, sr = librosa.load('audio/00020_2003_person1.wav', sr=None)

window_size = 1024
hop_length = 512
n_mels = 128
time_steps = 384

window = np.hanning(window_size)
stft= librosa.core.spectrum.stft(x, n_fft = window_size, hop_length = hop_length, window=window)
out = 2 * np.abs(stft) / np.sum(window)

plt.figure(figsize=(12, 4))
ax = plt.axes()
plt.set_cmap('hot')
librosa.display.specshow(librosa.amplitude_to_db(out, ref=np.max), y_axis='log', x_axis='time',sr=sr)
plt.savefig('spectrogramA.png', bbox_inches='tight', transparent=True, pad_inches=0.0 )

Spectrogram B
x, sr = librosa.load('audio/00020_2003_person1.wav', sr=None)
X = librosa.stft(x)
Xdb = librosa.amplitude_to_db(abs(X))
# plt.figure(figsize=(14, 5))
librosa.display.specshow(Xdb, sr=sr, x_axis='time', y_axis='hz')

Spectrogram C
# Read the wav file (mono)
samplingFrequency, signalData = wavfile.read('audio/00020_2003_person1.wav')
print(samplingFrequency)
print(signalData)

# Plot the signal read from wav file
plt.subplot(111)
plt.specgram(signalData,Fs=samplingFrequency)
plt.xlabel('Time')
plt.ylabel('Frequency')

Spectrogram results are displayed below:

My question is, from the 3 spectrograms I have listed above, which spectrogram is best to use for input to CNN and why should I use that spectrogram type? I am currently having difficulty to find their differences, as well as their pros and cons.

I recommend a Synchrosqueezed Continuous Wavelet Transform representation, available in ssqueezepy. Synchrosqueezing arose in context of audio processing (namely speaker identification), and there's much literature on applying CWT for audio tasks.
Advantage over STFT is the inherently logarithmic nature of the feature extractor, matching audio structuring. That is, we perceive sound differences in relative terms, so 6kHz is to 3kHz what 24kHz is to 12kHz. STFT is unable to map such features without being very wasteful, as it lacks an adaptive window (your Spectrogram A is not logarithmic; it's a logarithmic display of linear data). Analytic CWT (CWT with analytic wavelet) further enjoys superior instantaneous frequency and amplitude representations (see below references).
Synchrosqueezing further enhances these representations, and can be thought of as an attention mechanism. Experiments have shown it to enhance EEG classifier performance, for example. Currently only CWT is supported, but I'll be adding STFT within a week or so. Also see:

Excellent CWT tutorial (also comparing with STFT)
Synchrosqueezing friendly overview
CNN + CWT vs STFT comparison
Synchrosqueezing paper

Examples on EEG data:


'xcorr' function in matlab not working for zero signal I came across an interesting problem when using 'xcorr' in MATLAB. The expected signal and generated signal both are at zero DC. So, I expect correlation coefficient to be 1 or perfectly correlated. To my surprise, It is showing the highest error and time lag !!! Pls help !

Thanks,

the xcorr function is defined as
$$R_{xy}(m)=
\left\{
 \begin{array}{ll}
  \sum\limits_{n=0}^{N-m-1}x_{n+m}y^*_n  & \mbox{if } m \geq 0 \\
  R^*_{yx}(-m)& \mbox{if } m < 0
 \end{array}
\right.$$
The output is given as $$c(m)=R_{xy}(m-N)$$ where N is the length of the signal by the above definition the output for zero signal is zero


Fourier Series Representation of Continuous-Time Periodic Signals As a novice in signal processing, I have been going through Signals & Systems by Oppenheim to try and understand how continuous time periodic signals are represented by Fourier series coefficients.

The book defines a pair of equations which represent the Fourier series of a periodic continuous time signal.

I cannot understand when one should use the synthesis equation (3.38) and when one should use the analysis equation (3.39) to  determine the set of Fourier series coefficients.

Can someone give examples/explain as to when I should use one equation over another? Why is this so?

When you are given a continuous-time periodic signal $x(t)$ and you want to find out the corresponding CTFS (continuous-time Fourier series) coefficients $a_k$ associated with $x(t)$, then you use the analysis equation; i.e, analyse $x(t)$ to find out $a_k$ use

$$ \boxed{ a_k  = \frac{1}{T} \int_{0}^{T} x(t) e^{-j \frac{2 \pi}{T} k t} dt  }\tag{1}$$

On the other hand, when you are given a set of CTFS coefficients $a_k$ and you want to obtain the corresponding continuous-time periodic signal $x(t)$, then you use the synthesis equation; i.e., sum up exponentials weighted by $a_k$ and create $x(t)$ as in

$$ \boxed{ x(t) = \sum_{k=-\infty}^{\infty} a_k e^{j \frac{2\pi}{T} k t } } \tag{2}$$


Circular Convolution in Scilab I am quite familiar with signal processing using Matlab. Now I am working in its open source alternative 'Scilab.' In Matlab, cconv() function is used to perform circular convolution between two sequences. I want to do circular convolution in Scilab. But I could not find an equivalent function for circular convolution in Scilab. Are there any other methods to perform circular convolution in Scilab (or, is there really a function for circular convolution in Scilab)?

I can't tell about SciLab but if it is not implemented, you can always get the circular convolution by

IDFT(DFT(x1) * DFT(x2))

where x1 and x2 are your signal vectors and the multiplication is element-wise. If the signals have a different length, zeropad the shorter one to the length of the longer one. This assumes that you want to do circular convolution with a modulus equal to the length of the longer vector.


Frequency conversion in IIR lowpass-to-lowpass transformation We use
$\ z^{-1} \xrightarrow{} \frac{z^{-1}-\alpha}{1-{\alpha}z^{-1}}$
and
$ \alpha = \frac{\sin(\omega_c^{'}-\omega_c)/2}{\sin(\omega_c^{'}+\omega_c)/2}$
for a lowpass-to-lowpass frequency transformation, where $\omega'_c$ is the cut-off frequency of the prototype lowpass filter, and $\omega_c$ is the cut-off frequency of the transformed lowpass filter.
In practice we know the desired edge frequencies $\ \omega_s$ and $\ \omega_p$ and we are required to find the prototype lowpass cutoff frequencies $\ \omega_s^{'}$ and $\ \omega_p^{'}$.
For calculating $\ \alpha$ we have to know $\ \omega_c^{'}$ . If we calculate $\ \alpha$ we can have $\ \omega_s^{'}$ by $\ \omega_s^{'} = \angle{(\frac{e^{-j\omega_s}+\ \alpha}{1-\alpha e^{-j\omega_s}})}$
But I don't know how select a suitable $\ \omega_p^{'}$.

The prototype filter is not uniquely defined by the two band edges of the filter to be designed. So you have to choose one of the two edge frequencies of the prototype filter. E.g., if you choose $\omega'_p$, you can compute the value of $\alpha$ according to
$$\alpha=\frac{\sin[(\omega'_p-\omega_p)/2]}{\sin[(\omega'_p+\omega_p)/2]}\tag{1}$$
where $\omega_p$ is the desired pass band edge. Now that you have $\alpha$, the stop band edge of the prototype filter can be computed from the desired stop band edge $\omega_s$:
$$\omega'_s=-\angle\left(\frac{e^{-j\omega_s}-\alpha}{1-\alpha e^{-j\omega_s}}\right)\tag{2}$$

EDIT: the lowpass-to-lowpass transformation
$$z^{-1}\longrightarrow\frac{z^{-1}-\alpha}{1-\alpha z^{-1}}\tag{3}$$
transforms a prototype edge frequency $\omega'_s$ in the following way:
$$e^{-j\omega'_s}=\frac{e^{-j\omega_s}-\alpha}{1-\alpha e^{-j\omega_s}}\tag{4}$$
From $(4)$ it follows that
$$-\omega'_s=\angle\left(\frac{e^{-j\omega_s}-\alpha}{1-\alpha e^{-j\omega_s}}\right)\tag{5}$$
which is equivalent to $(2)$.
Note that Eq. $(5)$ is different from Eq. $(8.71)$ in the book because $(8.71)$ is derived from a lowpass-to-highpass transformation.


In What Way Is the Difference of Gaussian (DoG) More Tunable than the Laplacian of Gaussian (LoG)? I was reading through this post "What is the difference between Difference of Gaussian, Laplace of Gaussian, and Mexican Hat wavelet?" and noted in Chris Luengo's answer that he mentioned that the Difference-of-Gaussians (DoG) had an additional property of being tunable in specific way different to that of the Laplacian of Gaussian (LoG) - I've been a lurker for quite some time, so I don't have the privileges to comment, unfortunately.

My intuition is that the DoG has two sigma parameters (standard deviations) that can be varied independetly, wheres the LoG only has one (the initial Gaussian), and that this gives the DoG some bandpass flexibility that the LoG does not have. Have anyone looked into this further? Is there any research dealing directly with this issue? What I want to know are some more specific details in which way the DoG is tunable compared to the tuning properties of the LoG.

Thank you for your attention.

I have delved deeper into computer vision since first asking that question, and I have a better handle on real-world use-cases of DoG vs LoG. All of this signal processing basically boils down to looking for specific features, when you do not know a priori the scale that feature exists at. First, let us look at SIFT, Scale-invariant feature transform, specifically scale-space extrema detection.

SIFT

The scale space of a n-dim signal $f(x_0,x_1,\dots,x_n)$ is the family of derived signals $L(x_0,x_1,\dots,x_n; t)$, where $t$ is a tunable parameter corresponding to the scale of the filter, $t = \sigma^2$, or the variance of the Gaussian filter. To be more concrete, let us take a digital image $f(x, y)$ and create a 3D array $f(x, y, t)$ of $n_x \times n_y$ pixels and $n_t$ "slices", where each slice is a blurred copy of the original image for $t=\{0\dots n_t\}$, the first being 0 blurring, a.k.a. original image. Now for SIFT (simplification), apply an edge/corner detection to each slice, and correlate each point $corner(x,y)$ through $t$. When you have a strong number of these points detected, you can be confident that this is a keypoint of the image.

If you take the difference of any two different slices of $f(x,y,t)$, you get a DoG for $t_1$ and $t_2$.

Check out this page to see some examples of scale-space representation

SURF

SURF, or Speeded-up Robust Features, is an improvement on SIFT. This is a fantastic paper, and it, along with some of the papers it cites in Related Work, go much of the way to answer your question of "Have anyone looked into this further?"

Here we get into the heart of the matter, why one might use LoG/DoG/Wavelet for a particular problem. Take a look at fig 2 on page 3:

These kernels approximate second-order Gaussian derivatives, which are objects normally of infinite width and resolution. However, because we are dealing with digital images, and also want to be fast, they are cropped, quantized, and discretized. However, once again, there are many choices for how we choose to crop and quantize, with pros and cons for each.

Wavelets

Wavelets are yet another approach to scaling, edge detection, etc. Wavelets have some nice mathematical properties which are beyond this scope. The picture above shows successive stages of discrete wavelet transform (DWT). The simplest wavelet is the Haar wavelet, $[1,-1]$, which is essentially the 1D discrete derivative, very similar to numpy.diff(array, axis=n). You apply your wavelet in one axis, which gives you a detail image $d_0$ (the edges) and an average image $a_0$. $a$ and $d$ are basically low-pass and high-pass outputs, respectively. In contrast to continuous filters, each time you apply a DWT, you get 2 outputs, each with half the number of pixels along the given axis. With 2d images, you will typically take a $m \times n$ image, apply a DWT in the vertical and in the horizontal, to get 4 images, $(d,d)$, $(d,a)$, $(a,d)$, and $(a,a)$. Your $(a,a)$ is the low-pass output in both directions, at half the resolution. You can apply DWT again, and again, recursively, each time getting an "octave" in scale space.

One obvious advantage of DWT is that each low-pass output has $\frac{1}{4}$ the number of pixels to process. The downside is you are limited to only octaves of scale space.

Continuous wavelet transform behaves in much the same way as Gaussian scale-space, except you specify not only the scale, but the filter kernel as well. In discrete images, different wavelets have different properties. The Daubechies wavelet is quite popular, because, insofar as I can tell, it works on magic.

deep breath Ok bear with me, almost done here.

Deep Convolutional Neural Networks

Feature-detecting filters in scale space show up in deep learning, but in this case, the convolutional kernels are learned, rather than specified by the scientist. Now, instead of tuning just $\sigma_1$ and $\sigma_2$, you have the entire network worth of architectures, parameters, and hyperparameters to fiddle with. Also your results are dependent on your image dataset. Your bandpass properties are completely beholden to your NN. But I think it's really cool that the NN is learning this of its own accord. No one said "let's have a neural network for images which uses successive edge detection for feature detection", they just created the convolutional neural network architecture, trained it on image classification tasks, and those filters just emerged as optimal solutions to the task.

Which is best? Of course, there is no single best for a given application. But here are some pros and cons:

LoG: Simple, but has to be calculated for each scale. The discrete * LoG is a good approximation of the continuous. Basically a single parameter to tune.
DoG: similar to LoG, but when generating a scale space, there is a lot of computational redundancy that can be leveraged to save time. Two parameters to tune, or more precisely, you chose window widths in scale space.
DWT: Crazy fast, sometimes faster even than FFT. Limited to octave scaling. Choice of many discrete wavelets.
CWT: similar to LoG, but with more choices of kernels, thus more knobs to twiddle.
CNN: Great for real world datasets, but no direct over control tuning. Requires intensive training.

There is also a connection to Bessel functions to all this, if you want yet another rabbit hole.

(I've edited this answer and written it in several stages since there is a lot of content)


Image Shadow Removal Using OpenCV and Python How do I remove a shadow after MOG2 background subtraction using OpenCV Python? I used all morphological operations, gaussian and median blur, thresholding. it removes noises but deep shadow is resulting in foreground object.

fgbg = cv2.createBackgroundSubtractorMOG2(128,cv2.THRESH_BINARY,1)
masked_image = fgbg.apply(image)

in masked_image shadow will be grey color(pixel value= 127)
just replace 127 to 0, to convert grey pixel to black.

masked_image[masked_image==127]=0


Fourier Transform and Delta Function I am very new to Fourier analysis, but I understand that through the use of the Fourier transform a signal in the time domain is displayed in the frequency domain, where frequency values are normally displayed along the x-axis, and amplitude is displayed along the y-axis.  However, at one point in the textbook I am using, the following is stated:

Let us assume that we have the function $f(t) = \cos(\omega_0 t)$.  The spectrum then consists of two delta-functions

$$F(\omega) = \pi \delta(\omega - \omega_0) + \pi \delta(\omega + \omega_0)$$

__

This confuses me.  When we have $f(t) = \cos(\omega_0 t)$, then I would assume that the Fourier transform should yield an amplitude of $1$ at $\omega = \omega_0$ and $0$ elsewhere. But the delta function is defined as:

$$\delta(\omega - \omega_0) = \left\{ \begin{array}{1 1} \infty & \quad \omega = \omega_0 \\ 0 & \quad \omega \neq \omega_0 \end{array} \right.$$

So wouldn't this give an infinite value at $\omega = \omega_0$?

If anyone can explain the intuition behind the statement in my textbook, then I would be very grateful!

The textbook is right. A sine wave in the time domain has infinite energy since it continues over an infinite amount of time. When you transform into the Frequency domain all this energy is concentrated on a single (or two) frequency. Hence the value there is indeed infinite.

Of course these are all theoretical considerations. In the real world, ideal sine waves do not exist since they all have a beginning and an end.


parseval for a continuos but limited signal I have a question about the parseval relation written here
https://en.wikipedia.org/wiki/Parseval%27s_theorem (In the chapter Notation used in physics).

If I have a signal continuous but limited (so it does not go from $\infty$ to $\infty$ but from 0 to T, can The Parseval theorem be applied?

Sure. You can can just integrate in the time domain from $0$ to $T$ since the area outside  $[0, T]$ is 0.

Please note that you still must integrate from $-\infty$ to $+\infty$ in the frequency domain since finite support in the time domain implies infinite support in the frequency domain (and vice versa).

Things are different if the function is periodic.


Get the lag of a signal with respect to the other with cross correlation I have 2 vector of length 20 each. The first 15 entries of $a$ are exactly the same as the last 15 entries of $b$. There is random noise in the 5 last entries of $a$ as well as the 5 first entries of $b$. I would like to use cross correlation to detect the lag of 5 between those 2 signals and wrote the following code:

%Matlab Code
clear all; close all;
iteration = 128;
suma = 0;
for it = 1:iteration
  offset = 5;
  range = 100;
  vecLength = 20;
  a = randi(range,1,vecLength);
  b = zeros(1,vecLength);
  b(offset:size(b,2)) = a(1:size(b,2)-offset+1);
  noise = randi(range,1,offset-1);
  b(1:offset-1) = noise;
  c = xcorr(a,b);
  %[c,lags] = xcorr(a,b);
  [argvalue, argmax] = max(c);
  argmax = argmax - length(a);
  suma = suma + argmax;
end
lag = suma/iteration

However, I get an average of around -2.0 for the value of lag. Here are the result of 7 runs:

-1.7344
-2.1875
-2.1719
-2.1641
-1.9531
-2.3750
-2.1172

How can I get the back the value of offset = 5 instead?

The offset is actually 4 samples given the code generated since MATLAB starts indexing at sample 1.

The accumulation from the DC offset (the non-zero average value) creates a larger triangular correlation pattern and together with the random component can sometimes overcome what would otherwise be the strongest correlation at an offset of 4. The solution is to remove any average level prior to doing the correlation.

Given the uniform distribution between 1 and 100, the average will be close to 50 (and random itself so will vary from run to run, but the mean of the random process is 50.5). When both sequences align, the correlation due to the offset alone would on average be $20(\bar x ^2)$, where $\bar x$ is the average value. This result would be close to 50,000 in this case. Since fewer and fewer of the samples are included in the correlation process as the lag increases, this value will linearly decrease to 0 with each shift of the samples.  Due to the random element each sample will vary from this triangular shape. Similarly the correlation of the sequence itself when aligned would on average be $15(\bar x^2)$. By removing the mean we eliminate the DC offset effect and allow the repeated sequence to dominate the correlation result when those sequences are aligned.

This demonstrated in the plots below plotting the 20 correlation runs individually of the OP's code to show the variability, using

[c,lags] = xcorr(a,b);
plot(lags,c)

and

[c,lags] = xcorr(a-mean(a),b-mean(b));
plot(lags,c)

The other option derived from the cross-correlation theorem is to compute a circular cross-correlation using the FFT given as:

$$Corr = \text{ifft}(\text{fft}(a)\text{fft}^*(b))$$

Since it is a circular correlation, the correlation from the DC offset is the same throughout the correlation and not strongest at zero lag as in the linear correlation so need not be removed in this case.

Here is the general MATLAB code for providing the circular correlation results with FFTs given two sequences of equal length with $N$ samples each. Since the inverse FFT will provide a result starting at time =0 and is circular in time (such that the last sample is also the result for the first negative sample), the fftshift command is used to center time = 0 in the middle of the plot similar to the plots above.

# lags is different if odd or even length
if (mod(N,2)); lags = -(N-1)/2: (N-1)/2; else; lags = -N/2:N/2-1; end
c = fftshift(ifft(fft(a).*conj(fft(b))));

As RBJ points out in the comments, this results in time aliasing, as the correlation is effectively done by wrapping around to the beginning rather than correlating to zeros at the ends of the sequence as in the linear correlation that is performed using xcorr.  The sequences can be zero-padded (append 20 zeroes to each sequence) prior to taking the FFT to replicate the linear correlation result, but in this case the mean value would need to be removed for the same reasons given above. Given the correlated samples are random noise I actually prefer the circular correlation personally without zero-padding and removing the mean as it provides a constant noise floor in the result (rather than tapering to zero as you can see in the first two plots above).


Why is the last value of an RFFT always real? I am using numpy to do FFTs of real-valued data. And I don't understand why the Nyquist frequency is always real (or has zero phase).

So, say A = rfft(data) then A[-1] is always a real value, and not complex. Is this the correct value for that frequency? Or is this a computational artifact that can be fixed?

The explanation given in the documentation is:

  A[-1] contains the term representing both positive and negative Nyquist frequency (+fs/2 and -fs/2), and must also be purely real.

Is there a way to extract only the positive frequency and get rid of the degeneracy in the imaginary component caused by the Hermitian nature of the real FFT?

The real component of an FFT or DFT represents cosine or even sinusoids (symmetric around the middle), and the imaginary component represents sine, or odd sinusoids (anti-symmetric around the center).  An anti-symmetric waveform has to be zero at the center of an even length window.  If a sinusoid of frequency Fs/2 is zero at any sample, it is zero at all samples taken at a sample rate of Fs.  The imaginary component of the basis vector at N/2 for a DFT transform is also all zeros for the same reason. Thus the imaginary bin at Fs/2 is the sum of zero times zero, and has to be zero.  Leaving the result in bin N/2 strictly real for even length DFTs.  Perhaps this is the bin your RFFT implementation reports for A(-1).

The is also the reason why Nyquist sampling has to be at a rate ABOVE the_highest_frequency / 2, not equal to.  Sampling will leave out part of any anti-symmetric waveform or waveform component (in any even length set of samples) if that component is of exactly frequency F when sampled at F/2.

There are different rules for odd length FFTs.


Comparison/Difference between two distinct plots I am making a project where I  can identify leaves by taking a picture of them and a database. I will use their widths to determine their shape. So, I worked across the stem of the leaf each width and plotted them onto the following graphs.

The data of my graphs are located at: https://docs.google.com/spreadsheets/d/1HnYYA9keX7jjkhp4tocd2gy-i8VtHzYaS28fax-Nofc/edit?usp=sharing

To compare two graphs, I take the minimum of two functions:

The sum of the absolute difference between two graphs for every X-axis  point
Reverse one graph, and then run another sum of absolute difference between the reversed and the other graph. This is in case if one of the two leaves is measured in an opposite direction from the other one.

i.e. min( sum(abs(a-b)), sum(abs(a’-b)), where n’ is the reverse plot of n (In case a plot is reversed).

So, the problem is, the difference between two very similar leaves is actually greater than two very distinct leaves.

(On top) I have three graphs, two of which are similar-pos,pos2 (which should return a smaller difference), both of which are the acer ginnala leaf and another is different-neg, the betula alleghaniensis (which should return a greater difference than the "similar"):

However, this is not the case. On the 2nd layer, the two graphs are very distinct. (one is acer, other is betula). However, running the above algorithm, the sums are 70648 and 27362, returning a result of 27362.

On the last layer, when two very similar graphs are compared, (both acer) they return a much higher difference, with 40084 (normal) and 85664 (reversed), and returns the minimum of 40084, which is a higher difference than the first comparison.

I’ve tried standard deviation of the differences, squaring differences, Intersection Correlation, but had no luck. So, since the current measurement metric doesn’t work on this, what is a concise, straightforward, and better way to measure the difference between two graphs?

One way to compare similarly-shaped data that is translated in time is to look at the correlation between them.

If I do this in scilab:

ddata = diff(data',1,'c');

XC = xcorr(ddata(1,:),ddata(2,:));
CC1 = xcorr(ddata(1,:),ddata(1,:));
CC2 = xcorr(ddata(2,:),ddata(2,:));

den = sqrt(max(CC1)*max(CC2))

clf
plot(XC/den);
plot(CC1/max(CC1),'r')
plot(CC2/max(CC2),'m')

disp(max(XC)/den)

then I get this:

where the blue plot is the (normalized) cross-correlation of the time-differenced data and the red and magenta plots are the (normalized) auto-correlation of each individual time series.

The result is that the "correlation" between them is 0.8416782.


Why use complex envelope to model raw measurements? The usual answers when asking "why use complex numbers to model measurements?" is usually: (a) Because you are measuring a phenomenon like electricity and magnetism. (b) Because it is super useful for analysis- Fourier, convolutions, differentiation: these are easier on the complex plane.
When discussing seismology, acoustics, and other real-valued measurements, authors still use complex numbers, even if analysis in the time-domain, and no convolutions involved.
The usual argument is that the complex numbers represent the complex envelope of the real valued measurement. But why do that? Is this a sociological/communal habit, or is there a practical advantage in working with the convex envelope of a signal, and not its real version?

'Imaginary' is a misnomer; they 'exist' as much as reals do. As for the envelope: it provides a meaningful representation of phase and amplitude over time, not otherwise doable with reals alone.
To illustrate, consider a real $x(t) = a(t)\cos(\phi(t))$, where $a(t)$ is amplitude and $\phi(t)$ is phase at time $t$. Suppose we observe $x(t)$, but don't know its function (i.e. righthandside of above). For any $t=t_0$, how do we tell its amplitude or phase? At $t_0$, all we have is $x(t_0)$, which can be at any point in the cycle:

Enter the analytic representation:

with amplitude & phase computed as:
$$
a_a(t) = \sqrt{\Re^2 + \Im^2},\ \phi_a(t) = \tan^{-1}(\Im / \Re), \\
\Re = \text{real}(x_a(t)),\ \Im = \text{imag}(x_a(t))
$$
Example with trickier $\phi(t)$:

But I seek $x$'s amplitude & phase, not $x_a$'s: that's the thing, they're the same: $x_a(t)$ is constructed precisely so to keep $x(t)$'s $a(t)$ and $\phi(t)$. Not only can we now compute these quantities, but also instantaneous frequency, and with some luck, defy Heisenberg's limitations. Particularly once entering the discrete realm however, the method isn't flawless, but works well enough in general.

Why is this appropriate for measured data? Note that we're not transforming the data in any fundamental way, but rather expressing it in an equivalent but more useful form. So we don't say "H2O is liquid", which assumes certain temperatures and pressures, but more like "H2O has two hydrogen and one oxygen atoms", i.e. exact same thing but reformulated.
Whether the results are meaningful, however, does use assumptions: we assume the data is mainly composed of (1) periodic (repeating over time), and (2) sinusoidal processes. For acoustics, both hold very well; air compresses and expands periodically and sinusoidally over a local time segment. Then it's meaningful to ask: what's $a(t)$ and $\phi(t)$ of this physical sinusoidal process?
Assumption (3) is equally crucial: the data has a single intrinsic mode or 'major frequency component'. For example, if we have two A.M. sines, $f=2\text{Hz}, 4\text{Hz}$ out of phase, then we need two amplitudes and phases, $a_1(t), a_2(t)$ and $\phi_1(t), \phi_2(t)$, so a simple analytic transform won't suffice. For this we use methods like STFT, CWT, and synchrosqueezing to decompose the signal into its components.
Lastly, assumption (2) is flexible: if signal is, say, triangular waves instead of sinusoidal, then we simply re-interpret resulting $a(t)$ and $\phi(t)$.


BER and distance relation I have a system using 16-QAM, with 2.4 GHz, working up to 2 km at different bit rates, e.g. 1Mbit/s.
My task is to provide graph, which will show relation between BER and distance for different bit rates. However, based on that article:
https://www.mathworks.com/help/comm/ug/bit-error-rate-ber.html#bq4230l and adding bit rates as errors in dB, I am able to provide the graph for BER vs. Eb/N0, I wonder- what is relation BER to distance?

In short, increasing the distance results in a decrease in $E_b$.
Consider two antennas in an otherwise empty universe (this is known as "free space"). The transmitter antenna radiates $P_t$ watts, and has gain $G_t$ in the direction of the receiver. The transmitted signal has wavelength $\lambda_0$. Then, the receiver antenna, with gain $G_r$ and at distance $d$ from the transmitter, receives a signal with power given by $$P_r = P_tG_rG_t\left(\frac{\lambda_0}{4\pi d}\right)^2.$$
You are assuming transmission at 1 Mb/s. Then, $E_b$ is given by $$E_b = \frac{P_r}{10^6}\frac{\text{joule}}{\text{bit}}.$$ You can see that, as $d$ increases, both $P_r$ and $E_b$ decrease. This causes an increase in the bit error rate.
In environments other than free space, the calculation of $P_r$ is different, since the presence of reflectors, obstacles, etcetera will influence how much power is received. A popular model for urban/suburban/rural environments is the Okomura-Hata model, which is valid for the range of frequencies, distances and antenna heights commonly used in cellular communications. There are many other models that apply in different scenarios, and with varying degrees of complexity and accuracy.


What Are the Differences between Super Resolution, Denoising and Deblurring? In the fields of computer vision and image processing, what are the differences between Super Resolution, Denoising and Deblurring?

All 3 of them fall into the category of Inverse Problems in the Image Processing world.

Lets assume a Linear Model and then we will be able to show all 3 of them as parameters of the same framework. Then the differences will be clear and one could generalize it into Non Linear settings as well.

Have a look on the following model:

$$ y \approx A x + w $$

Where $ x $ is the image we want to infer from the given image $ y $, $ A $ is the degradation operator applied on $ x $ and $ w $ is white additive noise.

Denoising
For denoising only noise is added hence $ A = I $ and we left with estimating $ x $ from a noisy measurements.
Deblurring
In that case $ A $ is a matrix form of some Low Pass Filter (Circulant Square Matrix) which applies a blur on the image. We try to infer $ x $ from a blurred and noisy version of it given by $ y $.
Super Resolution
In this case take $ A $ from (2) and remove few rows of it. Namely we get only some of the data in $ x $ and the data is both blurred and noisy and we try to infer the whole data of $ x $ given $ y $.

Usually the problems above aren't well defined which means they require some kind of regularization.

Have a look on the following problem:

$$ \arg \min_{x} \frac{1}{2} {\left\| A x - y \right\|}_{2}^{2} + \lambda r \left( x \right) $$

This is the Linear Inverse problem where you can adapt to any of the cases above. The $ r \left( \cdot \right) $ is the regularization function which means we have some model on the data. For instance for Images many use the Total Variation model which means the Image is Piece Wise Smooth (Sparse Derivatives).


How to remove text region in a document image? or extract image from a document image Given a document image (i.e. newspaper), how to extract photos in it or remove text region?

I think traditional OCR methods may not be suitable here, as I don't need to recognize the text, and OCR is not accurate and slow. I believe text region (i.e. text blocks) and image region should be distinguishable by some threshold based methods in image processing. Any suggestions or example codes in OpenCV will be appreciated. Thanks！

BTW, what if the background color is not white, or the background color of certain blocks are not white?

Example image:

What you are asking is called 'page layout analysis' in OCR. Take a look at https://en.wikipedia.org/wiki/Document_layout_analysis.


Phase difference between signals sampled at different frequencies I want to know that if it is possible to measure the relative phase difference between a signal that has been sampled at two different locations with different sampling frequencies? Also can that method be extended to undersampled cases as well?

Edit: Adding Matlab script to test possible solution (Eq.3) provided by Dan Boschen

clear all
close all
clc

Len = 768/121e6;
Fs1  = 157e6;
t1 = 0:1/(13*Fs1) :Len-1/Fs1; %Time vector for Channel 1
Fs2 = 121e6;
t2 = 0:1/(13*Fs2) :Len-1/Fs2; %Time vector for Channel 1

f=25e6; % Incoming signal frequency

phase_diff_in=0; % Modelling the actual phase difference taking In-Phase for now

% Creating signals
sign1 = cos(2*pi*f*t1);
sign2 = cos(2*pi*f*t2 + deg2rad(phase_diff_in) );
sign1 = sign1(1:13:end);
sign2 = sign2(1:13:end);

% Adding a reference cosine
sig_ref=cos(2*pi*Fs1*t2);% Fs1 sampled by Fs2
sig_ref =sig_ref(1:13:end);

% Test of phase difference formula in time domain
phi1=acos(sign1(1:256));% In first window of 256 points
phi2=acos(sign2(1:256));
phi3=acos(sig_ref(1:256));

T1=1/Fs1;
n=0:255;
phase_diff=2*pi*n*f*( ((T1*phi3(n+1))/(2*pi*n)) -T1)...
    - (phi2(n+1) - phi1(n+1));
phase_diff=wrapToPi(phase_diff);
figure;plot(rad2deg(phase_diff),'-*r')

As far as I understood the phase difference in this case should have been 0 but that is not the case. The phase difference (in deg) is as shown below:

Update: Simulating the code added by Dan

Fs1  = 157e6;
Fs2 = 121e6;
f=500e6;%25e6
samples = 400;
Len = samples;
Phi = 45;
phase_out=phase_scale(Fs1,Fs2,f,Phi,Len);
figure;
plot(phase_out)
mean(phase_out)

for the case when f=25e6 and phi=45 the following was obtained:

And for the case when f=500e6 and phi=45 the following was obtained:

The error increases significantly as the frequency is increased further.

Update #2: Simulation results after the code modifications by Dan

for the case when f=25MHz and phi=45 the following is obtained:

Which shows that the phase difference was measured very accurately.

Also for the subnyquist case as well @f=600MHz and phi=75, the following is obtained:

which shows that this works in the subnyquist cases as well. Hence the given solution works under the assumptions stated by Dan in 'Practical Limitations' section of the answer.

SOLUTION

Bottom Line

$$(\theta_2-\theta_1) = 2\pi f(T_2-T_1)n -(\phi_2[n]-\phi_1[n]) \tag{1}$$

$f$: frequency in Hz of two tones of the same frequency and fixed phase offset

$(\theta_2-\theta_1)$: phase difference in radians of tones being sampled

$T_1$: period of sampling clock 1 with sampling rate $f_{s1}$ in seconds

$T_2$: period of sampling clock 2 with sampling rate $f_{s1}$ in seconds

$\phi_1[n]$: phase result from sampling tone with $f_{s1}$ in radians/sample

$\phi_2[n]$: phase result from sampling tone with $f_{s2}$ in radians/sample

This shows how any standard approach of finding the phase between two tones of the same frequency that are sampled with the same sampling rate (common phase detectors approaches including multiplication, correlation etc) can be extended to handle the case when the two sampling rates are different.

Simpler explanation first:

Consider the exponential frequency form of equation (1):

$$e^{j(\theta_2-\theta_1)} = e^{j2\pi f(T_2-T_1)n}e^{-j(\phi_2[n]-\phi_1[n])} \tag{2}$$

The term $e^{j2\pi f(T_2-T_1)n}$ is the predicted difference in frequency between the two tones that would result from sampling a single tone with two different sampling rates (when observing both on the same normalized frequency scale).

The observed difference in frequency between the two tones would be $e^{j(\phi_2[n]-\phi_1[n])} $.

Both terms have the same frequency with a fixed phase offset. This phase offset is to the actual difference in phase between the two continuous-time tones. By conjugate multiplication we subtract the two, removing the phase slope and the fixed phase difference results.

Derivation

The approach is to carefully work with all units with a time axis of samples. The frequency domain is thus in units of normalized frequency: cycles/sample or radians/sample corresponding to cycles/sec or radians/sec when the time axis is seconds. Therefore our sampling rate, regardless of what it is in time given in seconds, will be always equal to $1$ cycle/sample (or $2\pi$ radians/sample if working in normalized radian frequency).

The two signals with the same analog frequency once sampled each with a different rate in the time domain, will be two signals each with a different normalized frequency.

This simplifies the problem to gives us the following result:

Given our original signals as normalized sinusoidal tones at the same frequency with different phase offsets:

$$x_1(t) = \cos(2\pi f t + \theta_1) \tag{3}$$
$$x_1(t) = \cos(2\pi f t + \theta_2) \tag{4}$$

Once sampled but still with time in seconds:
$$x_1(nT_1) = \cos(2\pi f n T_1 + \theta_1) \tag{5} $$
$$x_2(nT_2) = \cos(2\pi f n T_2 + \theta_2) \tag{6}$$

Equation (5) and Equation (6) expressed time in units of samples is:

$$x_1[n] = \cos(2\pi f T_1 n+ \theta_1) \tag{7}$$
$$x_2[n] = \cos(2\pi f T_2 n+ \theta_2) \tag{8}$$

Convert to complex exponential form so that we can easily extract the phase terms using complex conjugate multiplication, (for a single tone we just need to split the input signal into quadrature components; $\cos(\phi) \rightarrow [\cos(\phi),\sin(\phi)]\rightarrow \cos(\phi)+j\sin(\phi) = e^{j\phi}$, this is described using the Hilbert Transform as $h\{\}$)

$$h\{x_1[n]\} =e^{-j(\phi_1[n])} =  e^{2\pi f T_1 n+ \theta_1} = e^{2\pi f T_1 n}e^{\theta_1} \tag{9}$$
$$h\{x_2[n]\} = e^{-j(\phi_2[n])} =e^{2\pi f T_2 n+ \theta_2} =e^{2\pi f T_2 n}e^{\theta_2} \tag{10}$$

The complex conjugate multiplication gives us the difference phase term we seek and its relation to our measured results:

$$e^{-j(\phi_2[n]-\phi_1[n])} = e^{2\pi f T_2 n}e^{\theta_2}e^{-2\pi f T_1 n}e^{-\theta_1} \tag{11}$$

Resulting in

$$e^{j(\theta_2-\theta_1)} = e^{j2\pi f(T_2-T_1)n}e^{-j(\phi_2[n]-\phi_1[n])} \tag{12}$$

Note that $e^{-j(\phi_2[n]-\phi_1[n])}$ represents the measurement which for single tones will result in a frequency and this frequency is predicted to be $\omega = 2\pi f(T_2-T_1)n$, given by the  $e^{j2\pi f(T_2-T_1)n}$ term. If we remove the frequency offset (by the multiplication above), the result is the phase difference of the original signal.

Taking the natural log of both sides reveals the result in units of phase (radians):

$$(\theta_2-\theta_1) = 2\pi f(T_2-T_1)n-(\phi_2[n]-\phi_1[n]) \tag{13}$$

So in summary, $\phi_1[n]$, $\phi_2[n]$ come from our measurements given as $cos(\phi_1[n])$, $cos(\phi_2[n])$ and we establish the difference that we need to get our answer through the complex conjugate multiplication of the Hilbert Transform of those measurements.

Demonstration

I demonstrate this with the script below similar to the OP's configuration with the results plotted below, which now includes the decimation and was tested for both f = 25 MHz and f = 400 MHz (undersampled) with similar results This shows each step to demonstrate the process above, and the operations can be further combined. The Hilbert Transform in implementation would be any approach of choice to delay the sampled tones 90° (A fractional delay all-pass filter is a reasonable choice).

Len = 10000;
phase_diff_in = 45;
f=400e6; % Incoming signal frequency
D = 13
Fs1 = 157e6*D;
Fs2 = 121e6*D;
t1 = [0:Len-1]/Fs1;  % Time vector channel 1
t2 = [0:Len-1]/Fs2;  % Time vector channel 2
phi1 = 2*pi*f*t1;
phi2 = 2*pi*f*t2 + deg2rad(phase_diff_in);
sign1 = cos(phi1);
sign2 = cos(phi2);

% emulation of perfect Hilbert Transform for each tone:
c1_in = 2*(sign1 - 0.5*exp(j*phi1));
c2_in = 2*(sign2 - 0.5*exp(j*phi2));

% create expected phase slope to remove
n = [0:Len-1];
comp_in = exp(-j*2*pi*f*(1/Fs2-1/Fs1)*n);

% decimation
c1 = c1_in(1:D:end);
c2 = c2_in(1:D:end);
comp = comp_in(1:D:end);
pdout = c1.*conj(c2);
result = pdout.*comp;

%determine phase_diff
phase_out = rad2deg(unwrap(angle(result)));
mean_phase = mean(phase_out);

Below is the result for two test cases, 0° as the OP was trying in his example and then a 45° phase shift.

Below shows the result for the copies of the input signal at frequency $f$ sampled by $f_{s1}$ as sig1 and $f_{s2}$ as sig2 for the case of zero degree phase between them. The real of the complex conjugate product pdout is the bold red sinusoid, and we note that it has zero phase offset.

To confirm the calculations, the plot below compares it directly to the real of the compensation term $cos(2\pi f(T_2-T_1)) to see that they are the same frequency consistent with the equation.

And repeating with $\theta_2-\theta_1 = 45°$

The result of the raw phase data for every sample shows that each sample individually has extremely low noise (limited by numerical precision, so the result can be determined with very few samples!). Such performance will depend on the actual quality of the Hilbert transform to accurately delay the input tone by 90° to create a qaudrature copy. Under conditions of noise the result can be averaged to the extent of waveform stability for a very robust solution.

Extended testing of performance with the undersampling case shows excellent results (f = 400e6):

Testing every difference angle in 1 degree steps:

RMS error of 10,000 samples (Note vertical axis is in increments of 0.5e-11)

Result of a greatly extended frequency sweep of the input frequency from 1e6 to 4000e6 in steps of 1e6 with a 45 degree phase shift with 10,000 points measured at each frequency showed a consistent result for phase determination at all frequencies (oversampling and undersampling). This is with the OP's configuration with the two frequencies including the decimate by 13. (Thus the sampling rate of each of the input tones after decimation for this test was at fs = 157e6 and 121e6, thus the far right of this graph with the frequency of the tone being sampled being 4e9 is significantly under-sampled. The RMS error is proportional to the frequency of the tone as shown, but even under this extreme condition, the error is still less than 5e-10 degrees. (8.7e-12 radians or -221 dB).

Practical Limitations

The accuracy of the above result is limited by knowledge of the exact frequencies and phase relationship given by $f_{s1}$ and $f_{s2}$, and knowledge of the frequency $f$ of the tone being sampled.

(As written the solution also assumes that the two sampling clocks both start at time $t=0$, but the sampling offset can be added starting with equation (8) with a similar result; bottom line is the starting phase relationship between the two sampling clocks must be known or measured as it will introduce an additional offset).

The reality is that no two free-running clocks will stay in perfect synchronization; there will be an inevitable drift in the actual frequency and phase difference between the sampling clocks that are not locked to a common reference (see Segal's Law https://en.wikipedia.org/wiki/Segal%27s_law). One of the clocks must be declared our reference of time (and our measurement will be to the accuracy of that one clock). If the clocks are not physically co-located, two-way time transfer techniques (see https://tf.nist.gov/time/twoway.htm) can be used to measure one clock versus the other. If they are physically co-located, then the simple thing to do would be to sample one clock with the other.

Below I show how this approach can completely eliminate one of the sampling clocks from the equation for our solution:  (I haven't yet tested this so may contain math errors)

Consider sampling $f_{s1}=\frac{1}{T_1}$ with $f_{s2}=\frac{1}{T_2}$. This will ultimately remove $f_{s2}$ from the equation entirely by using $f_{s1}$ as the common reference (we essentially measured $f_{s2}$ with $f_{s1}$ by sampling $f_{s1}$ with $f_{s2}$ allowing us to put the samples of $f_{s2}$ in units of $f_{s1}$ counts.):

$f_{s1}$ as a cosine:

$$x_{s1}(t) = cos(2\pi f_{s1}t) \tag{14}$$

When sampled with $f_{s2}$ given the constraint they both start at t=0 becomes:

$$x_{s_1}(nT_2) = cos(2\pi f_{s1}nT_2) =  cos(2\pi nT_2/T_1) \tag{15}$$

Which in units of samples is:

$$x_{s_1}[n] =cos(2\pi T_2/T_1 n) \tag{16}$$

Resulting in a third phase measurement in units of samples that we can get by sampling $f_{s1}$ with $f_{s2}$ (importantly to be done at the same time $x_1(t)$ and $x_2(t)$ are sampled!):

$$\phi_3[n] = 2\pi T_2/T_1 n \tag{17}$$

Thus if we don't know $T_2$ but have $\phi_3$ we can substitute the above  equation  to get:

$$T_2 = \frac{T_1 \phi_3[n]}{2\pi n} \tag{18}$$

substituting into (4):

$$
\phi_2[n]- \phi_1[n] = 2\pi nf\bigg(\frac{T_1 \phi_3[n]}{2\pi n}-T_1\bigg) + (\theta_2-\theta_1)  \tag{19}
$$

Resulting in the following solution for the desires phase difference of the original input signals:

$$
\theta_2-\theta_1=  2\pi f\bigg(\frac{T_1 \phi_3[n]}{2\pi n}-T_1\bigg)n - (\phi_2[n]- \phi_1[n])  \tag{20}

$$

Where

$f$: frequency of tone being sampled

$T_1$: period of sampling clock 1 with sampling rate $f_{s1}$

$\phi_1[n]$: result from sampling tone with $f_{s1}$, values will be $cos(\phi_1[n]) $

$\phi_2[n]$: result from sampling tone with $f_{s2}$, values will be $cos(\phi_2[n]) $

$\phi_3[n]$: result of sampling $f_{s1}$ with $f_{s2}$, values will be $cos(\phi_3[n]) $

Thus by only knowing $T_1$ which is $1/f_{s1}$, we can measure $f$ from the samples of $x_1(t)$ directly, measure $\phi_1[n]$ by sampling $x_1(t)$ with $f_{s1}$, measure $\phi_2[n]$ by sampling $x_2(t)$ with $f_{s_2}$ and measure $\phi_3[n]$ by sampling  $f_{s1}$ with $f_{s2}$ and from those measurements resolve $\theta_2-\theta_1$.

Similarly if your application is for a phase offset that would not be changing, then you can measure $f_{s2}$ error using the slope of the result without having to sample $f_{s1}$ with $f_{s2}$.

The true results will depend on the actual clock accuracy of $f_{s1}$ but we have completely removed $f_{s2}$ from the equation. If you can consider $f_{s1}$ to be your true reference of time, meaning it is accurate enough for the precision and accuracy of your measurement, then the result will be the phase difference of your two waveforms being sampled. This means that ultimately you need something to be your common reference of time.


Variance of an Implicit Function of Kalman State Vector Given a state vector given by $ x = {[r, v, a]}^{T} $ (Range, Velocity, Acceleration) the Time to Hit is the the time which holds the following:
$$ r + v {T}_{tth} + \frac{a {T}_{tth}^{2}}{2} = 0 $$
Now, given the State Vector covariance $ P $ what would the the Time to Hit Covariance?

If the Time to Hit was given by a regular function it would be easy to do given its Jacobian, yet how can it be for this kind of calculation (Implicit Function)?

I think I have the solution.
I'd be happy to hear others' thought.

Defining $ F \left(r, v, a, {T}_{tth} \right) = r + v {T}_{tth} + \frac{a {{T}_{tth}}^{2}}{2} $ which is the implicit function which connects all variables.
Since we're dealing with non linear function the variance is given by:

$$ var \left( {T}_{tth} \right) = J P {J}^{T} $$

Where $ P $ is the state vector covariance at a given time and $ J = \left[ \frac{\partial {T}_{tth} }{\partial r} \frac{\partial {T}_{tth} }{\partial v} \frac{\partial {T}_{tth} }{\partial a} \right] $

Using Total Derivative Law on the Implicit Function:
$$ \frac{\partial {T}_{tth} }{\partial r} = -\frac{{F}_{r}}{{F}_{{T}_{tth}}}, \frac{\partial {T}_{tth} }{\partial v} = -\frac{{F}_{v}}{{F}_{{T}_{tth}}}, \frac{\partial {T}_{tth} }{\partial a} = -\frac{{F}_{a}}{{F}_{{T}_{tth}}} $$

Where $ {F}_{x} $ is the partial derivative of $ F $ with respect to $ x $.

Each individual component would be:

$$
\frac{\partial {T}_{tth} }{\partial r}=-\frac{1}{v+aT_{tth}}, \frac{\partial {T}_{tth} }{\partial v}=-\frac{T_{tth}}{v+aT_{tth}}, \frac{\partial {T}_{tth} }{\partial a}=-\frac{T^2_{tth}}{2\left(v+aT_{tth}\right)}
$$

Transform the state covariance matrix with $J P {J}^{T}$ and you'd have the right solution.

Given a diagonal state covariance matrix, we should have:

$$
P_{tth}=\frac{\sigma^2_r}{\left(v+aT_{tth}\right)^2}+\frac{T^2_{tth}\sigma^2_v}{\left(v+aT_{tth}\right)^2}+\frac{T^4_{tth}\sigma^2_a}{4\left(v+aT_{tth}\right)^2}
$$


Can't compute real cepstrum of real signal I'm trying to compute the real cepstral coefficients of recorded telephone audio in Matlab using the rceps function.

On some audio frames (480 samples per frame -- 60ms of audio at 8kHz), I get a Matlab error:  "rceps:ZeroInFFT". The cepstrum does not exist because some of the DFT coefficients are 0. The frames in question are not zero, nor are any of the sample values complex. Going by the RMS of the frame, many of the frames that have errors have significant energy. About 10% of my frames give this error, so it seems like a bigger problem than a fluke.

I'm confused because this type of analysis is very common, especially for speech analysis, but I can't find record of anyone else having this problem.

Any advice would be greatly appreciated.

audio % is a 236000x1 vector of doubles containing PCM audio data at 8kHz
window_size = 3 %Working with multiples of 20-ms frames
errs = zeros(N_frames,1);
for i=0:N_frames-window_size
    %% Cepstral analysis
    s = audio((i*160+1):((i+window_size)*160)); %Grab 60ms of audio -- 480x1 matrix
    try
        c = rceps(s);
    catch err
        errs(i+1) = 1;
    end
end

It's common, when computing a cepstrum, to replace any zero's or tiny magnitudes in the 1st FFT result with some (noise) floor value to keep the scale and range of the log function "reasonable looking".

Huge negative spikes (or -inf) from the log() of tiny spectrum magnitudes don't usually provide that much added useful information to the rest of the results in cepstral-type analysis.


Histogram operation for improving visual object tracking I'm developping a people tracking with nearest neighbour data association technique and kalman filter for smoothing/predicting.

The tracker works quite good but sometimes it mismatch objects if they cross them-self.. (I think it is a common problem for a tracker).

I really believe that this tracker could be strongly improved with also histogram matching (I'm using OpenCV and I've seen an easy implementation of Bhattacharyya distance that seems quite popular in this field).

The fact is that I'd like to keep an average histogram of the tracked objects (because objects slightly changes their appearance if they rotate for example) and I'm a little bit confused on that: how to continuously adapt the average-histogram-object when new observation arrives?

observations could be different in size.. how can I handle this situation? using the same bins for each observation (even if they are 10 pixels or 100 pixels?) and sum up the histograms? using a num of bins proportional to the number of pixels?  and what about doing the average of two histograms?

I think I need all histograms with the same num of bins, and maybe use a weighted average where the weight is proportional to the num of pixels, but I'd like to ask for some more accurate advice from you!

I don't think you have any choice other than to use the same number of bins for each observation.  Otherwise not only will you not be able to average the histograms, you will also not be able to compare them.

And you definitely need to change the histogram slowly, i. e.
$$h = (1 - \alpha)h + \alpha h_{obs}$$ where $h$ is your "moving average" histogram, $h_{obs}$ is the histogram of the current observation, and $\alpha$ is a small factor, whose value you will have to tune experimentally.

As far as different sizes of observations, the simplest thing to do is to impose a hard constraint saying that the size of the object cannot change too much from frame to frame.


How does cascading biquad sections for higher order filters work? I'm trying to implement an 8th order IIR filter and every application note and textbook I've read says that it is best to implement any filter of order more than 2 as second order sections. I used tf2sos in MATLAB to get the coefficients for second order sections which gave me a 6x4 coeffs for 4 second order sections, as expected.
Prior to implementation as SOS, the 8th order filter required 7 previous sample values to be stored (and output values as well). Now when implementing as second order sections how does the flow work from input to output, do I need to store only 2 previous sample values? Or does the output of the first filter feed as x_in in to the second filter and so on?

It's the last thing you said ("Or does the output of the first filter feed as x_in in to the second filter and so on?"). The idea is simple: you treat the biquads as separate second-order filters that are in cascade. The output from the first filter is the input to the second, and so on, so the delay lines are spread out among the filters. If you need to optimize the structure in a memory-constrained environment, you can note that adjacent biquads have redundant delay memory (i.e. the last few output samples of stage 1 are the same as the last few input samples of stage 2, so you don't have to store them separately like you would if you just implement the filters in isolation).


How to match zero-pole diagrams to their frequency responses (Discrete Time) I get confused when there are a lot of zeros/poles in the zero-pole diagram and I find difficulty understanding their frequency response.

I know the following:
1. Complex conjugates cause double peaks at their respective angles. The larger their magnitude, the sharper their peaks.

A single pole causes a peak at its respective angle. The closer it is to the unit circle the sharper the peak in the frequency response.
A single zero causes a dip at its respective angle. The closer it is to the unit circle the sharper the dip.

I am not sure how to deal with zero pole diagrams when they come in bundles as shown below where it is asked to match each zero-pole diagram to its frequency response.

My trials:

I matched:

Frequency Response 3 to Zero-Pole Diagram 2

Frequency Response 6 to Zero-Pole Diagram 3

These are all frequency selective filters. Note the zeros on the unit circle, which correspond to zeros of the frequency response. The filter's passband is determined by range of the angles of the poles, and the stopband is determined by the range of the angles of the zeros.

As an example, let's consider PZ-diagram $\#1$: since the stopband is around DC (i.e., angle zero), and there are no more zeros, just poles, at angles greater than approximately $\pi/3$ (and, by symmetry, smaller than approximately $-\pi/3$), this must be a high pass filter. The cut-off frequency must be close to $\pi/3$, so all frequency responses can be eliminated except for $\#4$ and $\#5$. However, frequency response $\#5$ has a zero at DC, which is not present in the PZ-diagram. Consequently, we can match PZ-diagram $\#1$ to frequency response $\#4$. I'm sure you can match the other diagrams yourself.


Power spectral density: where does the difference from decibel and regular representation come from? I am generating a simple DC signal with noise. As expected the fft shows that there is one peak at 0Hz and some noise. I am also trying to get the spectral power density representation of that signal, and I came to a problem.

When i am using a logarithmic scale, to show power to frequency relation in dB/Hz, I am observing a gradually decreasing curve (second curve in the picture below). It shows that the lower the frequencies the less power they have in relation to same value. I don't understand it  - there is supposed to be only one peak at 0Hz, all other frequencies should have the same power.

This phenomenon can be observed on a regular graph, on non-logarithmic scale. This other diagram (first curve) shows clearly that there is just one frequency with "energy".
Where does this difference come from?
Does it come from leakage? If so, then why it is not relevant for the non-logarithmic spectrum?

The matlab code below generates the signal and the graphs.

clear
home
close all

signal_multiplier = 3;
noise_bits =8;
signal_bits = 11;

fs = 1;
T = 1/fs;
t = 1:T:(2^signal_bits)-T;
noise_vec2 = randi([-(2^noise_bits) 2^noise_bits],1,length(t));

signal = signal_multiplier*t + noise_vec2;

NFFT = length(signal);
[P, F] = periodogram(signal,[],NFFT,fs,'power');
PdBW = 10*log10(P);

figure(1)
subplot(2,1,1);
plot(F,P)
grid on
title('PSD')
xlabel('Frequency (Hz)')
ylabel('Power/Frequency')
subplot(2,1,2)
plot(F,PdBW)
grid on
title('PSD using logarithmic values()')
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB/Hz)')

If you zoom in on the linear curve, you would see that it's also "gradually decreasing", just much steeper, than the logarithmic curve. The logarithm simply makes it more visible.
And, yes, what you see is indeed spectral leakage. It can be reduced by using a non-trivial window function as opposed to the rectangular window being employed by default in the periodogram function. Try using a Hann window, that should reduce spectral leakage.


What is the physics behind the width of a main lobe? We know that the square window gives the lowest main lobe width possible, and that other windows after that trade main lobe width for side lobe height. I also understand that the main lobe width is inversely proportional to the length of the DFT aperature.

Given this, is there then a physical limit to how narrow a main-lobe can be, for a fixed aperture length? For example, can we ever hope to have main lobes that are narrower than that given by a square window, (for a fixed aperture length again), via windowing or any other method?

A second related question I have here has to do with energy conservation - is the interpretation that windows trade side lobe height for main lobe width, because energy must be conserved? Thus any reduction in side lobe height means the energy "has to go somewhere" and ends up in the main lobe? Would this be a correct interpretation?

A non-rectangular bell-shaped window reduces information magnitude at the edges of the window.  This either is, or acts like, information loss.  Thus the main lobe gets fatter to indicate a greater uncertainty with respect to its peak due to information loss. (relative to the noise floor, etc. etc.)

The concepts of information theory seem to be used in both physics and signal analysis.

The side lobes get lower because the high frequency energy related to producing any discontinuities between the FFT's 2 sides using circular basis vectors is reduced by a bell shaped window.

The two effects are similarly opposite in term of energy trade-off, but not exactly opposite, as different windows have differing areas under the curve, thus producing another degree-of-freedom in the trade-off between main lobe and side lobe energy, even given Parseval's theorem.

The width of the main lobe is inversely related to the length of the information and S/N of the information fed to the FFT.  In absolutely zero noise, just 3 or 4 non-aliased points can be used to calculate the exact frequency of a single pure sinusoid in the time domain (e.g. equivalent to an interpolated main lobe width of zero in the frequency domain).  Additional information outside the DFT data can also be used to estimate a narrower main lobe: For instance, many FFT post-processing frequency estimation methods assume a-priori that the data is stationary outside the widow and thus make use of the results from other offset windows, etc.


Non-symmetrical lemniscate curve parameterization I'm trying to fit a function to data points. The data generally resembles a butterfly/lemniscate shape, see drawing.

The problem is that the shape in my data can be rotated, skewed and/or non-symmetrical.

I've been looking at Bernoulli's, Devil's curve, Watt's curve, however, these are, as far as I can see, symmetrical.

Does anyone know of a plane curve that is able to represent the example shapes? Preferably in Cartesian coordinates.

There is a polar asymmetric modelling for the lemniscate:

which you can rotate and scale more easily that in Cartesian coordinates (but easy to convert). Such parameterizations are used very often as they can be less troublesome to fit.

Similar curves in astronomy are also called analemmas:

.

From an image processing fit point of view, you can consult: A Unified Scheme for Detecting Fundamental Curves in Binary Edge Images, Asano, Tetsuo et al., Computational Geometry, 2001


Recover Frame Compression Lossiness from Many Similar Frames With many very similar frames of video that have been compressed with an unavoidably lossy encoding and available now as jpeg stills, are there available good mathematical techniques for recovering detail by considering multiple frames?  Certain areas of the image are of great interest and the artifacts there change from frame to frame.  Can they be cross-referenced to solve for higher detail?

The method to improve resolution by few frames is called Multi Image Super Resolution (As opposed to Single Image Super Resolution).
For compression artifacts you can look for JPEG Deblocking algorithms.


How to use Octave ellip() coefficients for 3-band C++ crossover I have been working on making a 3-band audio crossover in C++, with the goal of activating/deactivating a different device depending on whether the low, mid, or high range of the input is above a threshold. I got as far as using GNU Octave to generate and graph out the low, band, and high pass filter coefficients using the Elliptic IIR method. I'm pleased with the coefficient characteristics, but this is my first DSP software project so I don't know how to use them to actually implement the filter program.
I looked at a bunch of articles on the topic online, and I saw people say using the Biquad Direct type 2 form is the most efficient approach for this situation. When I looked at the math on Wikipedia I understood the concepts somewhat but became confused as to how to implement it, partially because the order of my filters is greater than 2, and partially because comparing the formulas to other sources it seems like there might even be a mistake in the Wiki source??
I feel like I am incredibly close to understanding what to do, but I need that final nudge in the right direction.
I'm not sure if it is helpful at this stage, but here is what I have made in Octave so far:
function iir_elliptic_version
  clf
  close all

  pkg load signal;
  sample_freq = 192000;

  [lp_nums, lp_denoms] = Gen_Low_Pass_Filter_Coeffs(159, sample_freq)
  [mp_nums, mp_denoms] = Gen_Mid_Pass_Filter_Coeffs(160, 1999, sample_freq)
  [hp_nums, hp_denoms] = Gen_High_Pass_Filter_Coeffs(2000, sample_freq)

  f = 0:sample_freq/2;
  W = f*(2*pi/sample_freq);
  H1 = freqz(lp_nums, lp_denoms, W);
  H2 = freqz(mp_nums, mp_denoms, W);
  H3 = freqz(hp_nums, hp_denoms, W);
  plot(f, 20*log10(abs(H1)), f, 20*log10(abs(H2)), f, 20*log10(abs(H3)))
  axis([0 4000 -60 20])
  title('Filter Frequency Response of Elliptic IIR filters')
  xlabel('Frequency (Hz)')
  ylabel('Magnitude (dB)')
  grid on
endfunction

function [nums, denoms] = Gen_Low_Pass_Filter_Coeffs( passband_end,
                                              Fs)
  %This function generates the coefficients for a low pass filter for the
  %low range given the passband end frequency, and the sampling rate.
  stopband_start = passband_end + 25;
  pass_ripple_dB = 5;
  atten_dB = 20;
  norm_pass_freq = passband_end/(Fs/2);
  norm_stop_freq = stopband_start/(Fs/2);
  [order, cutoff_radians] = ellipord(norm_pass_freq, norm_stop_freq, pass_ripple_dB, atten_dB);
  printf("The order of the low-pass is %d\n", order);
  [nums, denoms] = ellip(order, pass_ripple_dB, atten_dB, cutoff_radians);
endfunction

function [nums, denoms] = Gen_Mid_Pass_Filter_Coeffs( passband_start,
                                                      passband_end,
                                                      Fs)
  %This function generates the coefficients for a band pass filter for the
  %mid range given the passband start and end frequencies, and the sampling rate.
  pass_freqs = [passband_start, passband_end];
  stop_freqs = [passband_start - 160, passband_end + 170];
  pass_ripple_dB = 5;
  atten_dB = 20;
  norm_pass_freqs = pass_freqs/(Fs/2);
  norm_stop_freqs = stop_freqs/(Fs/2);
  [order, cutoff_radians] = ellipord(norm_pass_freqs, norm_stop_freqs, pass_ripple_dB, atten_dB);
  printf("The order of the band-pass is %d\n", order);
  [nums, denoms] = ellip(order, pass_ripple_dB, atten_dB, cutoff_radians);
endfunction

function [nums, denoms] = Gen_High_Pass_Filter_Coeffs( passband_start,
                                              Fs)
  %This function generates the coefficients for a high pass filter for the
  %high range given the passband start frequency, and the sampling rate.
  stopband_end = passband_start - 25;
  pass_ripple_dB = 5;
  atten_dB = 20;
  norm_pass_freq = passband_start/(Fs/2);
  norm_stop_freq = stopband_end/(Fs/2);
  [order, cutoff_radians] = ellipord(norm_pass_freq, norm_stop_freq, pass_ripple_dB, atten_dB);
  printf("The order of the high-pass is %d\n", order);
  [nums, denoms] = ellip(order, pass_ripple_dB, atten_dB, cutoff_radians, 'high');
endfunction

Not sure what your application requirements are but there seem to be some issues with your implementation.

Your pass band ripple of 5 dB is way too big. Something like 0.1dB would be much more typical
Elliptic filters are poor choice for cross overs. The create significant phase distortions in the cross over region that you can easily see if you sum your three bands. The sum is nowhere near "flat". Better choice would be a Butterworth filter for odd order (using the proper sign) or Linkwitz-Riley for even orders.
It also looks like you are trying to build a 3-band equalizer, not a cross-over. If that's the case the whole approach is wrong .
Your filter is WAY to steep which creates an enormous amount of time domain ringing. Your low-pass filter rings for more than 100 milliseconds.
Try to stay away from $[num,den]$ representation of filters. Use $[z,p,k]$ or $sos$ instead.
Implement filters as cascaded biquads. It's easy to do and there are dozens of code examples on the internet
Do NOT use "Direct Form II". That's pretty much the worst for audio processing. Use Direct Form I or Transposed Form II instead.

A quick explanation for the last three bullets. In audio processing, the poles of the filters are often extremely close to the unit circle. It's easy to get the transfer function coefficients from the poles and zeros, but the other direction is numerically difficult since it requires finding the roots of a polynomial where the roots are very close together.
In evaluating different filter topology it is useful to look at the transfer function from the input (or output) to the state variables. For Direct Form II (and transposed Form I) this is given by the pole-only transfer function, which can get extremely large (in audio). Even your seemingly harmless low pass at 159 Hz has a pole-only gain of a whopping 94 dB. That means that your state variables are 10000 larger than your input/output.


About speech preprocessing in mobile phones I understand that cell phones transmit speech signals in digital form through 2G/3G or other standards. But I could not find anything about the signal processing method lying in between the analog signals captured by the microphone and the final digital signal (which I think is not the complete speech signal but a reduced form) that is modulated over the respective modulation scheme. Is there any compression algorithm or some enhancement algorithm before the modulation stage?

Yes, the cellular phones use various forms of compression to convert the captured analog audio (speech) into the digital bitstream for transmission through 2G/3G/.../.

The specific method used depends on the GSM version which might dictate its own bandwidth constraint and backward or forward compatibility issues into the audio encoding stage.

Most typically a 2G/3G system utilizes a variant of voice code based on a linear predictive coding (LPC) scheme which converts the analog sound captured by the microphone into a set of coefficients (parameters) according to the vocal tract model of the speaker. Unlike PCM, mp3, ogg, ra, etc., this is not a waveform coder but a parametric voice coder therefore.

The speech is reproduced at the recevier side by using a similar code excited LPC (CELP) which uses various excitation signals for the vocal tract model to reproduce the input signal as natural as possible.

You can find more details in the literature by the keywords vocoder, LPC, CELP, GSM speech coding etc...


what is multiresolution I have gone to wavelets course and up to my knowledge, changing the scale and translation of wavelets produces multiresolution images.

is this correct or can someone tell me what is multiresolution analysis in images?

Does multiresolution directly describe wavelet or multiresolution can be achieved other than wavelets?

The concept of multiresolution has roots on the observation that in data (eg images), objects and features can be observed and processed at different resolutions. One of the crudest version is to build a pyramid with the original image, and other sub-sampled versions, keeping 1 out of 4, 1 out of 16... pixels.
This yields a crude multiresolution pyramid, which is not a wavelet per se.

Wavelets, in a strict sense, naturally yield a certain type of multiresolution, but satisfy other constraints (esp. invertibility), so they are only a quite useful way of doing multiresolution.

There are many fast ways to compute CWT, that differ from DWT, see for instance

Should I ever pick the continuous wavelet transform over the discrete one? DWT vs CWT vs STFT
Implementing Continuous Wavelet Transform


what is the relation between bandwidth and channel taps I've studied that when the bandwidth increases, the number of taps we consider for channel modelling increases(example in WLAN or LTE). why is that happening. why signals with lower bandwidth (example : voice channel) see only a single tap/ a flat fading channel?

regards,
phani tej

Just as you sample a signal to process or model it digitally, you also sample the channel's impulse response. Each sample corresponds to one tap. When the signal's bandwidth increases, you increase the sampling rate, so your channel model now has more taps.

Usually, a channel is modelled as having a single tap (flat fading) when the channel's impulse response is different than zero only for the first sample. You can then ignore all other taps, since they're zero. This means that the largest path delay through the channel is shorter than the sampling interval. The implication is that a symbol will not interfere with its neighbors.

Note that, if you oversample, you can get a channel impulse response with more than one tap and still have flat fading. Flat fading means that, at the sample corresponding to the symbol time (and all subsequent samples), the channel's impluse response is zero.


Does the determination of the ergodicity of a signal force any changes in methodology? In mathematics, the term ergodic is used to describe a dynamical system which, broadly speaking, has the same behavior averaged over time as averaged over space.  -from wikipedia

From the perspective of engineering a signal processing system, does knowing whether a signal is ergodic change the plan of attack for analyzing the signals?  I've always found this a fascinating concept, but I don't really know what do do with it once the determination is made.  What analysis methods would be more pertinent knowing this information?

The concept of space as used in the statement ascribed to Wikipedia is presumably that of the ensemble of all possible signals, not space as in "Space: the final frontier...".  In the absence of parallel universes in which different signals from the ensemble can be observed simultaneously, it is rare for engineers to determine that a signal is ergodic by comparing time averages to ensemble averages.  More commonly, signals are assumed to be ergodic so that the values of time averages can be used as estimates of ensemble averages.  For example, an engineer takes $1000$ samples of a stochastic signal and plots a histogram of values.  This looks vaguely bell-shaped and so the signal is modeled as a Gaussian random process whose mean is the mean of the $1000$ samples etc.  Ultimately, the proof of the pudding is in the eating.  If this Gaussian model leads to predictions that differ wildly from observed results, the assumption of ergodicity might be discarded, or maybe $10,000$ samples, or $100,000$ samples might be tried (for example, predicting traffic at $5:30$ pm based on observations of traffic between $2:00$ am and $3:00$ am is rarely accurate, but
observations over several weeks might lead to a better model).


Signal sparsity: with noise or without noise? In compressive-sensing, signal should be sparse. Is this with or without noise? When I differentiate signal, it is supposed to be sparse. But when I add noise on it, it isn't sparse anymore. Should this algorithm  work correctly for me or not?

Actually, the beauty of compressive-sensing lies in fact that it works for signals which are approximately sparse. Recovery algorithms in CS are based on L1 norm minimization and have moderately tolerate noise in measurements. However noise surely reduces probability of exact recovery and degrades quality of recovered signal. If there is noise in a sparse signal you might need to take more measurements for exact recovery.
Take a look at this for further information:
http://statweb.stanford.edu/~candes/papers/spm-robustcs-v05.pdf


Removing breaths using python My wife is trying her hand at making audio books and spends a lot of time cleaning up some voice actors breaths in their "finished" audio. We have tried some filters in Adobe Audition and it helps but not completely.

This is what I need help with specifically: I would like to make a python (it's the only language I'm comfortable with) program that will look through an audio file and remove all the breaths and replace them with just background room noise. I've never done anything in the audio field and I'm completely overwhelmed by all the different packages and concepts. Could someone please point me in the right direction?

I was going to post this as a comment, but then it became so long that I thought, this constitutes really as an answer in fact.

Your question reminds me of the lab-work in my MSc course on Adaptive Filters. We used "Wiener filter" to remove some unwanted background noise (wind and tire noise) from the recorded input to make the speech clearer. The Wiener filter is quite simple to begin with.

You can look at the wiener filter from the scipy package of python - scipy.signal.wiener.

Another advantage of Wiener filter could be that you don't have to do any additional computational step of taking the Fourier transform of the input. That is, Wiener filter can be applied to the time domain signal directly.

A short example of how to use Wiener would be:

from scipy.signal import wiener
import matplotlib.pylab as plt
plt.plot(t,signal,'k')
plt.plot(t,wiener(signal,mysize=55),'r',linewidth=3)
plt.show()

Going by the assumption that the "breathing sound" in your audio inputs would be quite uncorrelated to the actually speech signal, just like the wind and tire noise that we removed to make the speech clear in our audio data, I expect the filter to provided a reasonable output.

(Note: There were then additional processing steps added sometimes to improve the audio output like making the parameters of the filter adaptive. You have to play around with the parameters a bit, as in one case, I remember the Wiener filter produced an output like it was coming from a hollow pipe, if you chose/modified the parameters in  a particularly wrong way. But as a start, just applying the Wiener filter to your signal in Python would be a good first step. And you can improve from their onward, step by step.)

And if you are curious as to how this filter works then you catch up more on the math here.


Sampling of bandpass signal with bandpass filter Trying to solve the following problem about bandpass signal sampling.
I have a signal whose Fourier transform is such that $X(w)= 0$ if $w>w_h$ or $w<w_l$.
The reconstruction is done using a bandpass filter rather than a lowpass and asks me for the minimum sampling frequency $w_s$ given that $w_l > w_h - w_l$.
I followed the steps of the sampling theorem and conjectured that if $w_s> w_h-w_l$ I should be able to reconstruct the signal with a lowpass filter. Indeed, denoting with $B=(w_h-w_l)/2$ I will have a copy of the signal centred at $0$ and with extremes $-B/2, B/2$ and then another copy centred at $w_s$ and with extremes $w_s-B/2$ and $w_s+B/2$. So, if $w_s > B$ I should have no aliasing. Is it correct? A lowpass filter should recover the signal centred at $0$.
Now, for the reconstruction with a bandpass filter. I have that the filter is such that $H(w)= T $ if $w_l\leq w \leq w_h$ or $-w_l\leq w \leq -w_h$. The question is: does this filter allow me to reconstruct the original signal $x$? I think it does not. First, it would give me two copy of the signal (right?) and also, there is no guarantee that there exist an integer $k$ such that $k w_s = (w_h+w_l)/2$.
The last question is then: assuming that $w_l>w_h-w_l$ what is the smallest sampling frequency $w_s$ and largest sampling interval $T$ that allow me to reconstruct the signal?
Honestly I do not understand, the issue to me remains the same, why would that assumption be enough to reconstruct the signal?

Assuming that we're really talking about an analytic complex-valued signal with no negative frequency components, then a sampling frequency $\omega_s>\omega_h-\omega_l$ will guarantee that the shifted spectra don't overlap, i.e., there will be no aliasing. However, it's not guaranteed that you'll have an image of the spectrum centered at DC. This is only the case if the sampling frequency satisfies
$$k\omega_s=\frac{\omega_l+\omega_h}{2},\qquad k\in\mathbb{Z}^+\tag{1}$$
Eq. $(1)$ means that the sampling frequency must be an integer multiple of the signal's center frequency.
So in the case of an analytic signal, reconstruction with a bandpass filter is possible if there is no aliasing, i.e., if $\omega_s>\omega_h-\omega_l$ holds.
But since you're given the extra condition $\omega_l>\omega_h-\omega_l$ I strongly suspect that the signal is actually real-valued with a (conjugate) symmetric spectrum. In this case, this answer contains the relevant information. Since $\omega_l>\omega_h-\omega_l$, we have
$$n_{max}=\left\lfloor{\frac{\omega_l}{\omega_h-\omega_l}}\right\rfloor\ge 1\tag{2}$$
And, from Eq. $(1)$ in the quoted answer, the lowest possible sampling frequency must satisfy
$$\frac{2\omega_h}{n_{max}+1}<\omega_s<\frac{2\omega_l}{n_{max}}\tag{3}$$


Complex Conjugate Sinusoids in Forward DFT I hope this isn't such a dumb question, but I'm finally getting to grips with the inner workings of the DFT.

What I'm having trouble understanding is why the basis complex sinusoids in the "forward" DFT are conjugated, and why the inverse DFT uses non-conjugated complex sinusoids.

I've written a script in MATLAB that attempts to implement what I understand about the DFT (not FFT, I know a fast function for that already exists). I notice that based on this conjugation issue I'm having, when I run a real sinusoid of some frequency through my DFT script, then perform the inverse part I get back the original sinusoid but at what looks like a quarter cycle out of phase.

Can anyone clarify what's going on for me?

(here's my MATLAB code snippet):

%% setup a test sinusoid
Fs = 44100;
freq = 440.0;
amp = 1.0;

% setup col signal vector
x = generate_sine(freq, 0, amp, Fs, Fs);

% set length of DFT (# of bins/frequency resolution)
N = 2^13;

% get bin resolution from sampling frequency: acts as
% a Kth test frequency multiplier for DFTs shorter than
% the sampling frequency
bin_width = Fs/N;

% truncate length of input signal to size of DFT (bins)
x(N+1:end) = [];

% plot original, truncated signal
figure(1)
subplot(3,1,1)
title('original signal')
stem(x, '.');

%% perform DFT using N test sinusoids (bins)
% This for loop represents the summation process of the DFT.
%   -i is an iterator where:
%       } i = 0 frequency = DC,
%       } and i > 0 = Kth sinusoid frequency
%   -in each loop iteration a new test complex sinusoid is
%    generated based on i iterator which is used a the
%    frequency number.
%   -in each loop iteration the dot product of the current test
%    sinusoid is taken, then stored in a vector of frequency bins
for i = 0: N-1

    % set current test frequency (Kth frequency)
    k = bin_width * i;

    % generate a complex sinusoid vector
    cp_sin = complex(generate_sine(k, 0, 1, Fs, N) / (N/2), generate_sine(k, pi/2, 1, Fs, N) / (N/2));

    % store dot product result in bin number i
    bins(i+1) = dot(cp_sin,x);
    bins = bins.';
end

% plot DFT result
subplot(3,1,2)
title('DFT')
plot(abs(bins));
xlim([0 N/2]);

%% perform inverse DFT
for i = 0:N-1

    % set current test frequency (Kth frequency)
    k = bin_width * i;

    % generate a sinusoid
    cp_sin = complex(generate_sine(k, 0, 1, Fs, N) * (N/2), generate_sine(k, pi/2, 1, Fs, N) * (N/2));

    % perform inverse transform to get original signal back
    xi(i+1) = (dot(bins,cp_sin))/N;
    xi = xi.';
end

%% plot inverse DFT result
subplot(3,1,3)
title('IDFT')
stem(real(xi), '.');
ylim([-1 1]);

Really the DFT is just the inner product of your signal with complex exponentially of various frequencies.

Taking the inner product of complex vectors requires taking the complex conjugate.

Reading your code, I do not believe this is correct:

cp_sin = complex(generate_sine(k, 0, 1, Fs, N) / (N/2), generate_sine(k, pi/2, 1, Fs, N) / (N/2));

shouldn't it make a complex conjugate:

cp_sin = complex(generate_sine(k, 0, 1, Fs, N) / (N/2), -generate_sine(k, pi/2, 1, Fs, N) / (N/2));

where I've placed a "-" in front of the second generate_sine call.

Also, can you let us know what the arguments to generate_sine do. I'm not sure what the third and fourth arguments achieve?