How do I perform the sifting process in empirical mode decomposition? I am programming a voice activity detection algorithm and I found a paper that recommends the use of HHT. I have been trying to program it and understand how to calculate $m1(x)$, but the shifting process confuses me.

Could someone explain what exactly I need to do after I calculate the mean value between the two splines of the original data?

Here is what I found at this website:
https://www.clear.rice.edu/elec301/Projects02/empiricalMode/process.html

The first component h1 is computed:
$h_1=X(t)-m_1$
In the second sifting process, $h_1$ is treated as the data, and m11 is the mean of $h_1$'s upper and lower envelopes:
$h_{11}=h_1-m_{11}$
This sifting procedure is repeated $k$ times, until $h_{1k}$ is an IMF, that is:
$h_{1(k-1)}-m_{1k}=h_{1k}$

Why does it go from $h_1$ to $h_{11}$? I am really confused by this.

It's not shifting, it's sifting. After you find the $m_1$ curve, your subtract it from your signal $X$ to get $h_1= X(t) - m_1$. You're essentially taking away that signal component and are left with $h_1$, which is your new starting point of your next iteration. I assume you were confused by the terminology.

Sifting is a general term in signal processing related to separating out components of a signal one at a time. It is frequently used in context of wavelet decomposition since this process is very similar in this respect.

The reason the next function is $h_{11}$ and not $h_2$ is because you have to keep going with this process until the resulting $h_{1k}$ is an IMF function (criteria given on the top of the page). Once you've reached that point, you get your first IMF basis function $c_1$. Only at this point have you reached your first empirical mode.

Then you start over by finding $m_2$ and so on.


Question about Difference of Gaussian (DoG) algorithm I am recently learning about Computer Vision and I am having a trouble understanding Difference of Gaussian (DoG) algorithm. I get how the algorithm works in high level idea, but I am trying to implement my own and I am confused about some steps.
For instance, I am trying to create 5 blur level for each octave, and I am confused about which filter and sigma value applying to which image. Using Matlab, for the first octave, I created a filter and applied:

sigma = 0.5;
gauss = fspecial('gaussian', [5 5], sigma);
blur1 = imfilter(img, gauss, 'replicate');
dog1 = img - blur1;
%Next level
blur2 = imfilter(blur1, gauss, 'replicate');
dog2 = blur1 - blur2;

I am not so sure if this is how I need to apply? Do I apply gaussian filter to previously applied image? I also saw code using k*sigma. I am not sure what k means and how to apply? Oh and what value should I used for sigma? Is it in [0, 1] range or can be bigger than that? Could someone help me on this? Thank you very much.
Thank you so much!

Difference of gaussian is the difference in the output of two Gaussian filters with different blur amounts (sigma).

Sigma is the size of the Gaussian filter. A bigger sigma gives you a bigger amount of blurring. A good way to think about it is a Gaussian filter with variance sigma is very roughly like averaging 3 x sigma samples wide (or 3 x 3 in an image)

e.g. from wikipedia:

Very Important when making a Gaussian filter in MATLAB make sure the size of the filter is at least 6 x sigma. In your above code you have 5 x 5 which is fine for sigma = 0.5, but for sigma = 1 you would want 6 x 6 or bigger.

The k is simply a multiplier for sigma.

e.g.

sigma = 0.5;
gauss1 = fspecial('gaussian', round([10*sigma 10*sigma]), sigma);
sigma = 1;
gauss2 = fspecial('gaussian', round([10*sigma 10*sigma]), sigma);
blur1 = imfilter(img, gauss1, 'replicate', 'same');
blur2 = imfilter(img, gauss2, 'replicate', 'same');
dog2 = blur1 - blur2;

A more complete code example that allows you to set the number of octaves and the steps per octave:

%% Filter using DoG
stepsPerOctave = 5;
octaves = 4;
mult = nthroot(2,stepsPerOctave);

% Create blurry images
sigma = 0.5;
kernelSize = [10*sigma*2^(octaves),10*sigma*2^(octaves)]
for k = 1:octaves*stepsPerOctave+1
    disp(['Sigma is ' num2str(sigma)]);
    gauss = fspecial('gaussian', kernelSize, sigma);
    blur(:,:,k) = imfilter(I, gauss, 'replicate', 'same');
    imagesc(blur(:,:,k)); colorbar; title(['Gaussian ' num2str(k)]); pause;
    sigma = sigma * mult;
end

% Create DoG
for k = 1:octaves*stepsPerOctave
    dog(:,:,k) = blur(:,:,k+1) - blur(:,:,k);
    imagesc(dog(:,:,k)); colorbar; title(['DoG ' num2str(k)]); pause;
end


Zero padding in frequency domain, any tricks to speedup? there are many situations when one has to carry some operations on spectra, like splitting a spectrum Z wide in two Z/2 spectra, or vice versa merging two Z sized spectra in a 2Z sized one. The only solution without having to go back to time domain and then to freq domain again (which is stupid imho) is operating convolution in frequency domain.

One example: I have two Z sized spectra representing, respectively, a STFT signal frame and an impulse response, which were NOT zero padded to start from;  I want to convolve them non-circularly so I have to transform both spectra to double size spectra, zero padded and left justified; I have therefore to expand both spectra to 2Z sized ones by simply doubling the bin positions leaving the odd ones empty (resulting in a 2Z spectrum representing a repeating copy of the signal in the 2nd half); convolve with a Dirichlet kernel representing a binay square wave of unitary frequency (here is the point), so to zero the second half; then I can multiply the two 2Z sized spectra to obtain linear convolution; then carry the inverse operations to split the result in two spectra of size Z for subsequent overlap-save. This is a total of four convolutions in frequency domain.

I anticipate the answers: why not simply acquire zero padded frames of proper size at the origin? Well let's assume that I can't for many reasons. One is that I am working in a STFT context already adopting a standard frame format, like non-zero padded frames, where I get frames of a given format and am expected to return frames of the same size and format as result of my processing.

OTOH, if I avoided the frequency domain convolutions I would have to do IFFT on both spectra, zero pad in time domain, FFT again for multiplying the spectra, IFFT for splitting the result and then FFT again, for a total of seven fourier transforms, which is perhaps even more stupid...

So I am asking whether some math tricks exist to carry such operations with spectra as zero padding and splitting or any trick to speedup convolution in frequency domain without having to go back to time domain. Thanks

It looks like your choice is based on whether (approx) 5*NLogN (IFFTs+FFTs) is greater or lesser than N^2 (linear convolution).

You could use coarser interpolation methods (cubic, spline, etc.) to upsample the spectra but the error would likely be audible as a clear loss in quality.


Basics of leakage phenomena in DFT and its approximation with sinc function I am trying to understand leakage phenomena when using DFT. In the following I am reffering to Lyons: Understanding Digital Signal Processing (2nd ed.), p.69-71.

Having real valued cosine with $N$ samples as an input, the amplitude response of N-point DFT can be approximated by (in terms of bin index $m$)

$|Y(m)|\approx\frac{N}{2}\operatorname{sinc}(k-m)$,

where $k$ is number of cycles inside $N$ samples.

If $k$ is a whole number, there is only 1 bin non-zero.
Otherwise spectral leakage occurs. So far so good.

Now I am want to replicate this in Matlab using following code.

N = 32;                % No of samples
fs = 32000;            % sampling frequency
m = linspace(0,N-1,N);
k = 8.5;               % No of cycles within N samples
X = sin(2*pi*k/N*m);   % cosine signal
Y = fft(X);
f = linspace(0,fs*(N-1)/N,N);
stem(f,abs(Y));hold on

With $k$ being a whole number, no leakage is obtained.

However, using $k=8.5$ produces peaks at 8 kHz and 9 kHz which differ in magnitude, althought slightly ($|Y((8 \text{kHz})/f_\text{s})|=10.15$ and $|Y((9 \text{kHz})/f_\text{s})|=10.25$ ). In my reference book number 10.17 is given.

My questions are:

Shouldn't the two peaks be exactly the same? Why they are not?
Additionaly, I have plotted function $|Y(m)|$ with ''continious'' $m$. The bins of original spectra should be exactly on the resulting sinc function. However, farther away we go from 8.5 kHz, bigger the discrepancy. Aren't the DFT bins just ''uncovered'' parts of underlying continious spectra?

No, the two peaks don't need to be the same. The discrete-time Fourier transform (DTFT) (samples of which are computed by the DFT) is also not symmetric with respect to the sinusoid's frequency (as explained in Robert's comment). If you compute the DTFT of a truncated sinusoid

$$x[n]=\sin(\omega_0n)(u[n]-u[n-N])$$

you get

$$X(e^{j\omega})=\frac{1}{2j}\left(e^{-j(\omega-\omega_0)\frac{N-1}{2}}\frac{\sin\left(\frac{(\omega-\omega_0)N}{2}\right)}{\sin\left(\frac{\omega-\omega_0}{2}\right)}-e^{-j(\omega+\omega_0)\frac{N-1}{2}}\frac{\sin\left(\frac{(\omega+\omega_0)N}{2}\right)}{\sin\left(\frac{\omega+\omega_0}{2}\right)}\right)$$

The first term in the parentheses has its main contribution around $\omega=\omega_0$, whereas the second term has it around $\omega=-\omega_0$. Note that the two functions are no sinc functions; only for $\omega$ very close to $\omega_0$ (or $-\omega_0$) can they be approximated by a sinc function. And that still wouldn't take into account their overlap.

So your results are correct, and the plot below shows your DFT on top of the DTFT as computed above:

Here's the Matlab/Octave code for the DTFT plot:

N = 32;
k = 8.5;
fs = 32000;
w = 2 * pi * linspace(0,1,500);
w0 = 2*pi*k/N;
X = exp(-1i*(w-w0)*(N-1)/2) .* sin(N/2*(w-w0))./sin((w-w0)/2);
X = X - exp(-1i*(w+w0)*(N-1)/2) .* sin(N/2*(w+w0))./sin((w+w0)/2);
X = X/(2*1i);
plot(w*fs/(2*pi),abs(X))


Is there a closed form expression for noise other than white noise? For white noise, where the PSD is constant, we can model that as a simple Gaussian. Is there any way to model the probability density function of other types of noise (pink noise, red noise, blue noise, etc.), and if so, is there any way to derive such a thing?

This is usually done with filtering of white noise.

White noise generator + LP filter give you a red noise generator.
White noise generator + HP filter give you a blue noise generator

Design a filter with required frequency response and apply it to a white noise generator.


Pre-Warping is not working on bilinear transform New here and I just wanna ask the question of what happens when your pre warping doesnt work on your bilinear transform?
I am assuming my isnt working as part of my research I have read that the higher the frequency is the more non linear it becomes, so I am assuming that my frequency I chose just has to many non linearities that it wouldn't ever get re-aligned to the continuous function critical frequency. It tries  but doesnt get there.
Sampling Frequency = 44.410kHz
Nyquist Frequency = 22.205kHz
Critical Frequency = 20kHz
Ratio = Fs/Nyquist frequency =  20/22.205 = 90.1%
I noticed the closer you are critical frequency is to your Nyquist frequency when applying the bilinear transform it gets ugly. And it looks like the pre warping can only do so much.
The only way I see around this if I just changed my Sampling frequency to like 96kHz which would make my Nyquist frequency = 48kHz would be then the ratio = 42% this would ease the constraint of the non-linearities and have the pre warping affect work much better, but this raises the question is it possible to have a digital filter with a sampling frequency of 44.410kHz and a Fc = 20kHz? Because the way this is going I dont see it.

CODE:
s = tf('s');

Fs = 96e3;
Ts = 1/Fs;
Fc = 20e3;

GP = -3;
GS = -20;
WP = 2*pi*Fc;
WS = 251327.412;

n = ceil(log10((10^(-GS/10)-1)/(10^(-GP/10)-1))/(2*log10(WS/WP)));

WC_1 = WP/((10^(-GP/10)-1)^(1/(2*n)));
WC_2 = WS/((10^(-GS/10)-1)^(1/(2*n)));

LP = 1/((s/WC_1)^4+2.613*(s/WC_1)^3+3.414*(s/WC_1)^2+2.613*(s/WC_1)+1);
LP_2 = 1/((s/WC_2)^4+2.613*(s/WC_2)^3+3.414*(s/WC_2)^2+2.613*(s/WC_1)+1);

LP_44_410kHz = c2d(LP,1/44.410e3,['Method','tustin','PrewarpFrequency',WP]);
LP_96kHz = c2d(LP,Ts,['Method','tustin','PrewarpFrequency',WP]);

options = bodeoptions;
options.FreqUnits = 'Hz';
options.Title.String = 'With Pre-Warping';
bode(LP_44_410kHz,LP,LP_96kHz,options);

is it possible to have a digital filter with a sampling frequency of 44.410kHz and a Fc = 20kHz?

Yes it is possible... .
I use MIM (Magnitude Invariance Method) for this. Here's plot showing 2nd order LPF responses for fc=20kHz-25kHz (1kHz step) @ 44.kHz sampling:

MIM method works well for even lower sample rate implementations but, there are also issues not sorted out yet.
Here's Matlab implementation of the MIM and PIM (Phase Invariance Method) method.


How to compute camera pose from Homography matrix? Say I use only one calibrated camera. From this camera, I get images A and B. I know the homography between A and B, computed through OpenCV's findHomography().
I know the pose (rotation matrix R and translation vector t) of image A, and I need the pose of image B. Once I get it, I suppose I'll be able to compute every further poses of following images.
Do you know an implementation of computing B's pose? I found several articles on the web, but I couldn't find an easily implementable solution...

I found a nice implementation, using OpenCV: http://nghiaho.com/?p=1298

Even if my answer comes too late for you, maybe other people find this useful. I have the codes for an openCV Pose from Homography. I found the method at this really useful website, euclideanspace.

void cameraPoseFromHomography(const Mat& H, Mat& pose)
{
    pose = Mat::eye(3, 4, CV_64FC1); //3x4 matrix
    float norm1 = (float)norm(H.col(0));
    float norm2 = (float)norm(H.col(1));
    float tnorm = (norm1 + norm2) / 2.0f;

    Mat v1 = H.col(0);
    Mat v2 = pose.col(0);

    cv::normalize(v1, v2); // Normalize the rotation

    v1 = H.col(1);
    v2 = pose.col(1);

    cv::normalize(v1, v2);

    v1 = pose.col(0);
    v2 = pose.col(1);

    Mat v3 = v1.cross(v2);  //Computes the cross-product of v1 and v2
    Mat c2 = pose.col(2);
    v3.copyTo(c2);

    pose.col(3) = H.col(2) / tnorm; //vector t [R|t]
}

////


scipy.signal.spectrogram I want to do time-frequency analysis of a signal. So i decided to use spectrogram function from scipy. This function returns, three parameters f,t,Sxx. I am using below arguements:

x : my signal
fs : 1000 Hz
window : hann, with size = 300
nperseg : 300, same as window size
noverlap : 25

My question is about returned time segments of this function.

As far as I understand the first time segment for which we calculate fft must be 300 samples, then subtract 25 from it calculate fft for next window of 300 samples,and so on. So the signal of length 1s must be divided in 4 time segments. .275,.500,.775,1. , but the function returns :

f,t,Sxx=spectrogram(sig1,fs=1000,window=get_window('hann',300),noverlap=25,scaling='spectrum',nperseg=300)
t
Out[172]: array([ 0.15 ,  0.425,  0.7  ])

can someone please suggest if I am thinking wrong or there is a bug in spectrogram function.

As the returned value shows, the first spectrum is computed at time $t=0.15$ s. With a sampling rate of 1000 Hz, that means that the first analysis window is centred at sample $n=150$. So the first analysis frame spans the samples $n \in [0,300[$.

Because your overlap is set to 25 samples, the second analysis frame spans the samples $n \in [275,575[$, thus overlapping in its first 25 samples with the first analysis frame. The second analysis frame is then centred at time $t=0.425$ s or in samples $n=425$.

The third analysis frame again overlaps 25 samples with the second one. Thus, it spans the samples $n \in [550,850[$ and is centred at time $t=0.7$ s or in samples $n=700$.

The next analysis frame would span the samples $n \in [825,1125]$ and is not computed because it would go over the range of samples in the input signal.


Resolution of Discrete Fourier Transform is 1/T - Mathematical proof? In many articles I see that the frequency resolution of the Discrete Fourier Transform (DFT) equals Fs/N where Fs is the sampling rate and N is the total number of samples. Fs/N is equivalent to 1/T where T is the total time the signal was sampled.

So the longer the time of sampling took place the better the resolution of the spectrum (allowing to separate different tones)

I'm looking for the mathematical proof of the 1/T resolution as I do not see it directly from the definition of the DFT.

$\frac{Fs}{N}$ is not equivalent to $\frac{1}{T}$, since $Fs = \frac{1}{T}$,  $\frac{1}{T N} \neq \frac{1}{T}$ for $N>1$ (as it is implied here of course).

The point about the resolution of the Discrete Fourier Transform follows from its definition where an $N$ point transform, "divides the circle" in $N$ discrete points. As you can see $k n$ (the discrete frequency and sample "instance" variables) are multiplied by the constant $\frac{i 2 \pi}{N}$.

Here, it is still possible to "sample" a non-integer harmonic (say $k = 1.267$) but that would not improve the situation if you only have a limited number of true $n$ (samples).

Hope this helps.


Justification of bilinear transform I would like to understand the "justification" for the bilinear transform. The basic idea as I understand it is that by integration rule of Laplace transform we have for continuous $y(t)$:

$$\mathcal{L}\{y\}(s)=\mathcal{L}\{\int_0^ty'(x)dx\}(s)=\frac{1}{s}\mathcal L\{y'\}(s)$$
$$\Leftrightarrow \frac{\mathcal L\{y\}(s)}{\mathcal L\{y'\}(s)}=\frac{1}{s}$$

Using the trapezoid rule gives the approximation ($h$ is the sampling period):
$$y(t)=y(t-h)+\int_{t-h}^ty'(x)dx\approx y(t-h)+\frac{h}{2}(y'(t)+y'(t-h))$$

And taking the $z$ transform using the approximation gives
$$\mathcal Z\{y\}(z)\approx \mathcal Z\{ y[n-1]+\frac{h}{2}(y'[n]+y'[n-1])\}(z)$$
$$\approx z^{-1}\mathcal{Z\{y\}(z)}+\frac{h}{2}(\mathcal Z\{y'\}(z)+z^{-1}\mathcal Z\{y'\}(z))$$
$$\Leftrightarrow\frac{\mathcal Z\{y\}(z)}{\mathcal Z\{y'\}(z)}\approx \frac{h}{2}(\frac{1+z^{-1}}{1-z^{-1}})$$

Putting the above together we finally end up with the estimate
$$\frac{\mathcal L\{y\}}{\mathcal L\{y'\}}(\frac{2}{h}(\frac{1-z^{-1}}{1+z^{-1}}))\approx \frac{\mathcal Z\{y\}(z)}{\mathcal Z\{y'\}(z)}$$

Question 1. How do we get from the above estimate to the desired estimate

$$\mathcal L\{y\}(\frac{2}{h}(\frac{1-z^{-1}}{1+z^{-1}}))\approx \mathcal Z\{y\}(z)$$

Question 2. It is possible to use trapezoid rule on any linear ODE directly (iterating the rule for higher order derivatives as needed). However, interestingly this gives a different answer than the bilinear transform. The resulting coefficients $a_i$ are the same, however the $b_i$s are different.

So it appears that the bilinear transform distorts the start of the impulse. For example, for 1st order Butterworth lowpass filter we have from bilinear approximation that $b_0=b_1$. As such the impulse response ramps up before ramping down, while the continuous response and trapezoidal estimate only ramps down (see example below). Why are the $b_i$s given by the bilinear transform the "best" (e.g. most widely used) estimates (note: we only consider $b_i$s up to the order of the filter)?

Example:
Consider the 1-pole Butterworth filter given by the transfer function:

$$H(s)=\frac{1}{1+s/\omega}$$

The differential  eq. is given by:
$$y=-y'/\omega$$

By trapezoid method and the above eq.:
$$y(t+h)\approx y(t)+\frac{h}{2}(y'(t+h)+y'(t))=y(t)-\frac{h\omega}{2}(y(t+h)+y(t))$$
Collecting the terms we get the difference equation
$$y[n+1]=\frac{(1-\frac{h\omega}{2})}{(1+\frac{h\omega}{2})}y[n]$$

Which is the same as with bilinear transform, but lacking the term $b_1$.

Let me rephrase your derivation of the bilinear transform to clarify the result. For any $t_0<t-T$ we have

$$y(t)=\int_{t_0}^tx(\tau)d\tau=y(t-T)+\int_{t-T}^tx(\tau)d\tau\tag{1}$$

Approximating the right-most integral in $(1)$ by the trapezoidal rule we get

$$y(t)\approx y(t-T)+\frac{T}{2}\big(x(t)+x(t-T)\big)\tag{2}$$

Switching over to samples of $y(t)$ and $x(t)$ and using the sloppy but common and handy notation $y[n]=y(nT)$, we define from $(2)$ the discrete-time (DT) system:

$$y[n]=y[n-1]+\frac{T}{2}\big(x[n]+x[n-1]\big)\tag{3}$$

This system is now used as the discrete-time approximation to an integrator. In the $\mathcal{Z}$-transform domain it is described by

$$H(z)=\frac{Y(z)}{X(z)}=\frac{T}{2}\frac{z+1}{z-1}\tag{4}$$

So the bilinear transform replaces the transfer function of a continuous-time (CT) integrator $H_c(s)=1/s$ by the transfer function $(4)$, i.e., the bilinear transform uses the mapping

$$s\leftrightarrow \frac{2}{T}\frac{z-1}{z+1}\tag{5}$$

So far for part $1$ of your question.

I must admit that I'm not sure I understand part $2$ correctly. Could you give a simple concrete example showing how you obtain two different DT systems from a given CT system?

As for you final question, the bilinear transform is of course not the best way to transform a CT system to a DT system, simply because there cannot be one best way to do so. However, in some applications - such as optimal filter design - it is considered the most useful transform because it preserves optimality of the magnitude. This is the case because the frequency response of the resulting DT system is just a compressed version of the frequency response of the original CT system, i.e., all properties concerning the magnitude (e.g., ripple size, number of zeros in the stopband, etc.) remain unchanged. Just the frequency axis gets warped, but we can take this  warping into account when designing the CT system such that the resulting DT system implements the desired band edges.

Furthermore, the bilinear transform is the only direct mapping from the $s$-plane to the $z$-plane that preserves the frequency response in the sense that the $j\omega$-axis is mapped to the unit circle of the $z$-plane. This is not the case with the forward and backward Euler methods. Furthermore, the left half-plane of the complex $s$-plane is mapped to the inside of the unit circle of the $z$-plane, i.e., stability is preserved in the most complete way. This is not the case with the forward Euler method, and even though stability is also preserved with the backward Euler method, the backward Euler method maps the left half-plane to a small region (a circle centered at $z=\frac12$ with radius $\frac12$) inside the unit circle. This means among other things that certain frequency characteristics cannot be achieved with the backward Euler method.

EDIT: Concerning your example of a first order Butterworth filter

$$H(s)=\frac{Y(s)}{X(s)}=\frac{1}{1+\frac{s}{\omega_0}}\tag{6}$$

you can't ignore the input of the filter in the differential equation, so you get from $(6)$

$$y'(t)=\omega_0\big(x(t)-y(t)\big)\tag{7}$$

Now if you use $(7)$ in the trapezoidal approximation

$$y(t)\approx y(t-T)+\frac{T}{2}\big(y'(t)+y'(t-T)\big)\tag{8}$$

you obtain the same difference equation as if you simply had replaced $s$ in $(6)$ by $\frac{2}{T}\frac{z-1}{z+1}$. This must be the case because the bilinear transform simply is the trapezoidal rule (also called Tustin's method).


How can resonance be added to a fractional order Butterworth low pass I am interested in how to add resonance (Q) to the magnitude response of a Butterworth low pass when it is expressed in the form:

$$
G^2(\omega)=\frac {1}{1+\left(\frac{\omega}{\omega_c}\right)^{2n}}
$$

I know how to add resonance $Q > 0$ for $n >= 2$ i.e. include it in one (and only one) of the second second order terms (if the equation is expanded as a product of second order sections):

$$
H(s) = \frac{ \omega_0^2 }{ s^2 + \frac{ \omega_0 }{Q} + \omega_0^2 }
$$

However, I wonder whether there is a general formula for resonance $Q$ that is a function of $\omega$ that can be used for fractional $n$ such that it can be added to the transfer function $G$ or $G^2$ to yield a new magnitude response?

Given the formula for second order sections:

$$
s_k = \omega_c e^{\frac{j(2k+n-1)\pi}{2n}}\qquad k = 1,2,3,\ldots, n \qquad s = i\omega
$$

$$
H(s)=\frac{G_0}{\prod_{k=1}^n (s-s_k)/\omega_c}
$$

Use k = 1 for the resonance section, add $Q$ in the normal way i.e. calculate $H_{k=1}^2 = H_{k=1}(s) \overline{H_{k=1}(s)}$ where $s = i\omega$ and the denominator will be of the form $w^4 + cw^2 + 1$ where $c$ is a constant calculated from $s_k$. We add $Q$ to the first order term (note that this is a quadratic in $\omega^2$). Using any positive real $n$ in $s_k$ gives valid values for $c$.

Let $R^2$ be a modified $H_{k=1}^2$ with $Q$ added in to the first order term of the quadratic in $w^2$ then to add resonance to the original equation in the question multiply it by:

$$
\frac {R^2}{H_{k=1}^2}
$$

This effectively removes the $k = 1$ second order section and adds a resonant section back in. It turns out once you do this, $n$ can be any positive real value when used in $c$ but note this wouldn't work for fractional roots of the sections (you would have to use fractional calculus if you wanted the roots of the fractional equation).

Without resonance:

$$
c = 2 \cos^2(\frac{\pi (1+n) } {2n}) - 2
$$

With resonance:

$$
c = \frac{2}{Q^2}\cos^2(\frac{\pi (1+n) } {2n}) - 2
$$


Discrete-time Biquad filter: Relation between peak frequency and pole frequency This question arose from a discussion in the comments to this question and its answer. The author of that question discussed a discrete-time second order filter described by the following difference equation:

$$y[n]=b_0x[n]+b_1x[n-1]+b_2x[n-2]-a_1y[n-1]-a_2y[n-2]\tag{1}$$

with $b_0=b_2=(1-r^2)/2$, $b_1=0$, $a_1=-2r\cos(2\pi f_0 T)$, and $a_2=r^2$.

The transfer function corresponding to the difference equation $(1)$ is

$$H(z)=\frac{1-r^2}{2}\frac{1-z^{-2}}{1-2r\cos(\omega_p)z^{-1}+r^2z^{-2}}\tag{2}$$

where I use the abbreviation $\omega_p=2\pi f_0 T$.

In the original question quoted above, $f_0$ was called "resonant frequency", and the discussion inspired by that question was about whether or not $f_0$ equals the peak frequency of the filter's frequency response.

So the questions I pose are:

What is the meaning of $r$ and $f_0$ as used in the coefficients of Eqs. $(1)$ and $(2)$?
Does $f_0$ equal the filter's peak frequency (i.e., the frequency where the magnitude of the filter's frequency response attains its maximum)?

According to the Audio EQ Cookbook the transfer function of a second-order band-pass filter with 0dB peak gain is given by

$$H(z)=\frac{\alpha-\alpha z^{-2}}{(1+\alpha)-2\cos(\omega_0)z^{-1}+(1-\alpha)z^{-2}}\\=\frac{\alpha}{1+\alpha}\cdot\frac{1-z^{-2}}{1-\frac{2}{1+\alpha}\cos(\omega_0)z^{-1}+\frac{1-\alpha}{1+\alpha}z^{-2}}\tag{1}$$

where $\omega_0$ is the angular frequency where the peak occurs ("peak frequency"), and $\alpha=\sin(\omega_0)/(2Q)$ is a parameter determined by $\omega_0$ and the desired quality factor $Q$.

If we express the same transfer function $H(z)$ in terms of its poles we get

$$H(z)=k\frac{1-z^{-2}}{\left(1-re^{j\omega_p}z^{-1}\right)\left(1-re^{-j\omega_p}z^{-1}\right)}=
k\frac{1-z^{-2}}{1-2r\cos(\omega_p)z^{-1}+r^2z^{-2}}\tag{2}$$

where $k$ is some gain constant to be determined, $0<r<1$ is the pole radius, and $\omega_p=2\pi f_p$ is the pole angle, corresponding to the pole frequency $f_p$. Comparing the denominator of $(2)$ to the filter coefficients given in this question (and shown under Eq. $(1)$ in the question above), it is clear that the frequency $f_0$ in that question equals the pole frequency $f_p$ in our formulation. And $r$ in the question above is obviously the pole radius.

Comparing $(1)$ and $(2)$ we can express $\alpha$ and $\omega_0$ of Eq. $(1)$ in terms of the pole radius $r$ and the pole angle $\omega_p$ used in Eq. $(2)$:

$$\begin{align}&\frac{1-\alpha}{1+\alpha}=r^2  \Longrightarrow \alpha=\frac{1-r^2}{1+r^2},\qquad
k=\frac{\alpha}{1+\alpha}=\frac{1-r^2}{2}\\
&\frac{2}{1+\alpha}\cos(\omega_0)=(1+r^2)\cos(\omega_0)=2r\cos(\omega_p)\end{align}\tag{3}$$

The second line of Eq. $(3)$ shows that the pole angle $\omega_p$ and the peak frequency $\omega_0$ are not equal, and that the following relationship holds:

$$\frac{1+r^2}{2r}\cos(\omega_0)=\cos(\omega_p)\tag{4}$$

The following Matlab/Octave example illustrates the above. We choose $\omega_0=0.2\pi$ and $Q=1$ and design the corresponding biquad band-pass filter according to Eq. $(1)$. The figure below shows the magnitude of the resulting frequency response, and the locations of the peak frequency $\omega_0$ (red line) and of the pole angle $\omega_p$ (green line), both normalized by $\pi$.

w0 = .2*pi;
Q = 1;
al = sin(w0)/(2*Q);
b = [al,0,-al]/(1+al);                % numerator coeffs
a = [1+al,-2*cos(w0),1-al]/(1+al);    % denominator coeffs
p = roots(a);                         % poles
r = abs(p(1));                        % pole radius
r2 = sqrt((1-al)/(1+al));             % pole radius according to formula
abs(r-r2)                             % ans = 0
wp = abs(angle(p(1)));                % pole angle: wp = 0.17877*pi
wp2 = acos((1+r^2)/(2*r)*cos(w0));    % pole angle according to formula
abs(wp-wp2)                           % ans = 0
[H,w]=freqz(b,a,1024);
plot(w/pi,abs(H))
hold on, plot([w0,w0]/pi,[0,1],'r',[wp,wp]/pi,[0,1],'g'), hold off

EDIT:

In Eq. $(2)$ I assumed that there are two complex conjugate poles, because in the original question, which triggered this question and answer, the filter coefficients were chosen in such a way that only complex conjugate poles were possible (cf. the coefficients under Eq. $(1)$ of the question above).

However, the general transfer function in Eq. $(1)$ can also have a double real-valued pole or two different real-valued poles. It can be shown that a double real-valued pole occurs when the quality factor $Q$ is chosen as $Q=\frac12$. This case is still covered by Eq. $(2)$ with $\omega_p=0$ (i.e., a positive real-valued pole) or $\omega_p=\pi$ (a negative real-valued pole). For $Q<\frac12$ we get two different real-valued poles, and Eq. $(2)$ is not valid anymore. Instead, the transfer function becomes

$$H(z)=k\frac{1-z^{-2}}{(1-\beta_1z^{-1})(1-\beta_2z^{-1})}=
k\frac{1-z^{-2}}{1-(\beta_1+\beta_2)z^{-1}+\beta_1\beta_2z^{-2}}\tag{a}$$

where $\beta_1$ and $\beta_2$ are the two real-valued poles. Comparing this transfer function to the transfer function given in Eq. $(1)$, we get the following relationships between $\alpha$, $\omega_0$, $k$ and $\beta_1$ and $\beta_2$:

$$\begin{align}\alpha&=\frac{1-\beta_1\beta_2}{1+\beta_1\beta_2}\\
\cos(\omega_0)&=\frac{\beta_1+\beta_2}{1+\beta_1\beta_2}\\
k&=\frac{1-\beta_1\beta_2}{2}\end{align}\tag{b}$$

Of course, the pole angles are now either zero or $\pi$, because both poles are real-valued. Note that the peak frequency $\omega_0$ is not  necessarily $0$ or $\pi$. The peak frequency $\omega_0$ can in fact take any value, and the two pole angles will always be either $0$ or $\pi$ as long as the quality factor $Q$ satisfies $Q<\frac12$.


what does musical noise sound like and is there a sample clip for it? Musical residual noise is created when a noisy signal passes through a noise reduction algorithm. For instance, please see the paper:  Esch, T; Vary, P: "Efficient Musical Noise Suppression for Speech Enhancement Systems", IEEE Intl Conf Acoustics, Speech and Signal Proc, (ICASSP): 4409-4412, 2009.

  Noise reduction techniques that are relying on spectral weighting
  rules often generate annoying musical noise artifacts in the processed
  signal.

What exactly does it sound like? Are there any sample clips available (for example, a wav file before passing through a noise reduction algorithm routine and a wav file after noise reduction with musical noise) for comparison?

Here's one example: http://www.izotope.com/tech/aes_suppr/

Googling "Musical Noise" certainly turns up a lot of noise music, so I'd suggest looking for academic papers that have web links to sound examples.


Regarding the advantages of generalized linear phase filters I know that for a linear-phase filter with frequency response given by
$$H(e^{j\omega}) = |H(e^{j\omega})|e^{j\phi(\omega)} $$
if the input of the system is
$$x[n] = s[n]\cos[\omega_0 n] $$
where $s[n]  $ is a narrowband signal with bandwidth $W\ll\omega_0$, then
the output of the system is approximately
$$y[n] \approx |H(e^{j\omega_0})|s[n-\tau_g]\cos[\omega_0 (n -\tau_p)] $$
(if $\tau_g$ and $\tau_p$ are integers) where
$$\tau_g =\left.-\frac{d\phi(\omega)}{d\omega}\right|_{\omega=\omega_o} $$
$$\tau_p =\left.-\frac{\phi(\omega)}{\omega}\right|_{\omega=\omega_o} $$
A broadband signal can be thought as a superposition of narrowband signals, this means that if we have a filter with linear phase (i.e.  $\phi(\omega)= \tau\omega$)
both group delay and phase delay are constant and equal, therefore all the narrowband signal packets will be shifted in the same amount.
But what about generalized linear phase filters? In this case $$\phi(\omega) = n_d \omega + \phi_o $$
so the output of a narrowband input centered at $\omega_0$ will now be
$$y[n]= |H(e^{j\omega_0})|s[n-n_d]\cos\left[\omega_0\left(n -n_d-\frac{\phi_0}{\omega_0}\right)\right] $$
and the term $\frac{\phi_0}{\omega_0} $ will differ from one signal packet to another as the center frequency ($\omega_0$) of each packet changes.
So why are generalized linear phase filters as important  as linear phase ones? Isn't the term $\frac{\phi_0}{\omega_0} $ important?

The only related question I could find here did not have a complete answer to my question.

Note that with the definition of generalized linear phase $\phi(\omega)$ according to

$$H(e^{j\omega})=|H(e^{j\omega})|e^{j\phi(\omega)}\tag{1}$$

and

$$\phi(\omega)=\alpha\omega+\beta\tag{2}$$

the restriction that the impulse response $h[n]=\text{IDTFT}\{H(e^{j\omega})\}$ be real-valued only allows two possible values for $\beta$:

$$\beta\in\{0,\pi\}\tag{3}$$

This is due to the required conjugate symmetry of $H(e^{j\omega})$: $H(e^{j\omega})=H^*(e^{-j\omega})$. The value $\beta=0$ obviously results in a system without phase distortion, and the value $\beta=\pi$ results in a sign reversal without any other phase distortions.

However, if you define generalized linear phase according to

$$H(e^{j\omega})=A(\omega)e^{j\phi(\omega)}\tag{4}$$

with $\phi(\omega)$ given by $(2)$, and with a real-valued but possibly bipolar amplitude function $A(\omega)$, then for real-valued filters the constant $\beta$ can take on four different values:

$$\beta\in\{0,\pi/2,\pi,3\pi/2\}\tag{5}$$

For $\beta=0$ and $\beta=\pi$ the amplitude function $A(\omega)$ is even, and we get no phase distortion apart from a sign reversal for $\beta=\pi$. However, for $\beta=\pi/2$ and $\beta=3\pi/2$, $A(\omega)$ is odd, and there will be phase distortion, i.e., the group delay does not equal the phase delay. So for narrow-band signals, the envelope experiences a different delay than the carrier. Basically, - apart from the delay - the carrier is just shifted by $\pm 90$ degrees. This is what is desired when implementing differentiators or Hilbert transformers. Consequently, systems with generalized linear phase and with $\beta=\pi/2$ or $\beta=3\pi/2$ are important for implementing causal (FIR) approximations to differentiators and Hilbert transformers, where the phase shift is implemented exactly, apart from a delay necessary to make the filter causal.

In sum, systems with generalized linear phase generally distort the phase (unless $\beta=0$ or $\beta=\pi$), but these distortions are desired when implementing differentiators and Hilbert transformers, and - apart from a delay - they implement the desired phase exactly.

As a final remark, note that the approximation of the response of an LTI system to a narrow-band input signal of bandwidth $W$ using phase delay and group delay is only valid if the magnitude of the system's frequency response is (approximately) constant in the band $[\omega_0-W,\omega_0+W]$.


Hardware kits for DSP Lab? I have been assigned responsibility as Lab instructor of DSP at undergrad Level and so i have to equip my lab,please recommend which processors should i purchase for hardware implementation of dsp concepts

For an undergrad DSP lab, there's absolutely no benefit in complicating what hardware your students need to learn: Let them do DSP on a completely normal PC.
There's plenty of offline processing and real-time capable DSP platforms for PCs. I'd personally recommend getting a bunch of Linux PCs, and having Python3 with scipy/numpy/matplotlib on them for basic algorithms exploration on that.
For advanced topics, like doing live signal processing on audio coming from a microphone or from an internet stream, or doing software defined radio with real SDR hardware, I'd go for GNU Radio, which comes with a lot of useful building blocks, and can be easily extended using Python and C++.¹

¹ Disclaimer: I'm involved with the GNU Radio project, but I do, and I know several educators, use it in DSP/communications engineering labs. Say "hello" on the mailing list.


How does summation over time samples influences DFT? In case of DFT where we have following
$$
Y_{k}=\frac{1}{N}\sum_{n=0}^{N-1}y_ne^{-j{\frac{2\pi nk}{N}}}
$$

What happens in case when we change summation, e.g., increase or decrease number of $n$ that we are summation over?
How does the rest of equation change?
For example if we decide to sum over $NM$ (if $M<1$ we have windowing and if $M>1$ we have zero padding so I do not assume adding any new info data just scaling the duration of the signal) samples would we rewrite the latter equation like this
\begin{equation}
Y_{k}=\sum_{n=0}^{MN-1}y_ne^{-j{\frac{2\pi nk}{MN}}}\tag{1}
\end{equation}
or this
\begin{equation}
Y_{k}=\sum_{n=0}^{MN-1}y_ne^{-j{\frac{2\pi nk}{N}}}\quad \tag{2}
\end{equation}
And what is the explanation behind. Also what is the appropriate scaling when summation limits change?

If the length of your data is $NM$ long (regardless of $N$), your DFT should be $NM$ long, so the correct equation is (1).

The explanation is simple: how long is your data? Use that length for the DFT.

You may want to define:
$$
y'_n = \left \{ \begin{array}{cl}
y_n, &  0 \le n \le N-1\\
0, &  n \ge N\\
\end{array}
\right .
$$
and operate on that because $y_n$ for $n\ge N$ is not defined.


How is Bayesian Estimation related to filtering? I am reading about estimation theory, including topics like Bayesian Estimation (e.g. Wiener Filtering).

It seems that we usually define a filter in terms of tis frequency response (e.g. High Pass, Low Pass). On the other hand, Wiener Filter works by filtering out noise without specific reference to frequency domain.

How are the two related? If not, are they just two different approaches to the filtering problem?

In Wiener filtering, you filter a noisy signal to more closely resemble a desired signal that you have access to.  In Bayesian estimation, you take prior knowledge into account to estimate some state given noisy measurements.  In frequency filtering, you just remove frequency content from a signal.


Carrier frequency of chirp signal I am computing the link budget for a communication-radar (CoRadar) system. The formula for maximum distance, in communication, is:

$d_{max} = \sqrt{\frac{P_t}{\tau}\frac{G_t G_{r, comm} \lambda^2}{(4\pi^2)FkT_0 B_w SNR_{comm}}}$

where,

$P_t$ is the transmitted power

$\tau$ is the duty cycle

$G_t$ is the transmission antenna gain

$G_{r, comm}$ is the communication receiver antenna gain

$\lambda$ is the wavelength

$F$ is the noise factor

$k$ is the Boltzmann's constant

$T_0$ is the noise reference temperature

$B_w$ is the bandwidth

$SNR_{comm}$ is the $SNR$ of the communication system

I am using a LFM chirp signal that sweeps from $3.5\times10^6$ to $8.5\times10^6$. Which would be the frequency, which would be associated to a wavalength $\lambda$, that should I use for the formula? $8.5\times10^6$? Or, $\frac{8.5\times10^6 + 3.5\times10^6}{2} = 6\times10^6$ (i.e. a middle point)?

And which would be the bandwidth? $8.5\times10^6 - 3.5\times10^6 = 5\times10^6$?

Many thanks.

In my experience in the radar industry, most people use the center frequency of the chirp for such calculations; in your case, that would be 6MHz.  Since your bandwidth is only 5MHz this is probably fine. You'd start to run into issues with that calculation if your bandwidth was much wider, say several GHz.


How do I take the real part of this bandpass filter system's output? I'm stuck on a final step in this problem. Essentially, there are the two systems above, which we'll call System 1 (Fig. 4.26, with ideal lowpass $H(jw)$) and System 2 (with $H_1(jw)$). The question proposes that the output y(t) of System 2 is identical to that which would be obtained by retaining $\Re\{f(t)\}$ of System 1.

I have a couple questions -

(1) Why do we need to take the real part of system 1 and not system 2? Aren't they both already bandpass filters?

(2) I solved out answers for $F(jw)$ and $Y(jw)$. $H_1(jw)=H(jw)$ is assumed (they seem to be the same):

$$F(jw)=H(j(w+w_c))X(jw)\\Y(jw)=\frac{1}{2}(H(j(w-w_c))X(jw)+\frac{1}{2}H(j(w+w_c))X(jw)$$

I'm relatively confident this answer is correct, but I'm unsure of how to take the real part of $f(t)$. Can I take the real part of $F(jw)$ and get it to look like $Y(jw)$?

Thank you!

The signal $f(t)$ is a complex bandpass signal. Note that if $H(e^{j\omega})$ is an ideal lowpass filter (and $\omega_c>\omega_0$) then $f(t)$ only has positive frequency components. For a real-valued signal, the spectrum shows conjugate symmetry, so there have to be negative frequency components as well.

By taking the real part of $f(t)$ you get $y(t)$:

$$\Re\{f(t)\}=\frac12\left[f(t)+f^*(t)\right]\Longleftrightarrow
\frac12\left[F(j\omega)+F^*(-j\omega)\right]=\\=\frac12\left[X(j\omega)H(j(\omega+\omega_c))+X^*(-j\omega)H^*(j(\omega_c-\omega))\right]\tag{1}$$

(where '$\Longleftrightarrow$' means Fourier transform pair).

Since $x(t)$ and $h(t)$ (the impulse response of the ideal lowpass filter) are real-valued, their Fourier transforms satisfy $X(j\omega)=X^*(-j\omega)$ and $H(j\omega)=H^*(-j\omega)$. Using this property, Eq. (1) can be rewritten as

$$\Re\{f(t)\}\Longleftrightarrow\frac12 X(j\omega)\left[H(j(\omega+\omega_c))+H(j(\omega-\omega_c))\right]=Y(j\omega)$$


Dynamics processor for speech - apply gain on delayed samples? A Dynamics processor (limiter, compressor, expander, and gate) is calculating the output gain at some point in time, and it is applying this gain at some delayed point in time.

This allows for Dynamics processor to "look ahead" and act on some audio events a few mili-seconds before they actually arrive. This is usefull in cases where sharp, transient signals occur. However, my use case is speech, in a bit noisy environment, but nothing "sharp", no transients, no drum kicks and similar.

Would adding lookahead in such a "speech" use case be beneficial?

It sounds better this way.

Let's say you have limiter and you can't exceed a certain amplitude value (risk of clipping, for example). Now there is an audio event (say a snare hit) that peaks at twice that amplitude. If you just reduce this sample by sample instantaneously you create a lot of audible distortion. If you have look ahead, you can ramp down the gain smoothly so that by the time the audio even arrives, the gain is low enough so that you are not clipping any more.

Rapid gain changes tend to create a lot of audible artifacts. The look ahead buffer allows you to slow down the gain changes without the risk of creating an instantaneous out of bounds condition.

It's also useful to look at a bit in the future when determining the gain for the current sample: if that's the highest point in the waveform can slow down the gain reduction, but if it's part of a slope, you probably want to maintain or increase your gain reduction speed.

Of course, you are not actually "looking into the future". Look ahead simply means that you put a delay on the processing path but NOT on the detection path of the dynamic processor.


FIR filter usage / understanding FIR Let us assume that I wouldd like to create a digital lowpass with 5 taps. I am not using a window function but just the impulse response from my Fourier integral.

I get the following coefficients: $h=[0.5\, 0.86 \,1 \,0.86\, 0.5]$.

What I don't understand: how does the filtering work? Assuming that I have a constant signal  $x$ with '1' and I want to use the filter on that.

I applied the filter and what I get is

$$y=[0.5 \,
1.36\,
2.36\,
3.22\,
3.72\,
3.22\,
2.36\,
1.36\,
0.5
]\,.$$
But now i am clueless how to draw this? The output is bigger then my input signal?

The good response by Laurent Duval provides a short overview of convolution and its properties and contains most of the information that you need. I think though that I see that you have some confusion about what is meant by filtering in general and how to apply a finite impulse response (FIR) filter in particular. I am attempting to explain some intuition behind FIR filtering here.

Something to keep in mind is that the length of an FIR filter should be shorter (preferably much shorter) than the signal that you wish to filter. This is because the result of a discrete convolution of two sequences is equal to the sum of the lengths of the sequences minus one, and part of this output is the transient response of the filter.

All filters have a transient response that applies to sudden changes in the input sequence (for example, a step function or discontinuity in the input). The beginning of a signal sequence looks like a step function to the filtering process (at least if you initialize your tap delay line with zeros as is typical). The end of the sequence has a similar transient response while zeros are being shifted in to replace the samples in the delay line.

The transient response occurs over the the first $N-1$ samples while the incoming signal is replacing the zeros in the delay line memory of the filter. This means that those first $N-1$ samples are not part of the steady state output of the filter. When we refer to the filtered signal, we are typically neglecting this transient response and are instead referring to the steady state output of the filter. This is often justified when the signal length is much greater than the filter length (the typical case), although not in your example where they are equal. If you use a filter of length $N$ and a signal of length $N$ as you have done, then the filtered output has length $2 N - 1$, and you have only $(2 N - 1) - 2 (N - 1) = 1$ sample that is not part of this (typically undesired) transient response. Most further processing will require more than $1$ sample to be effective, so this is probably not practical for you. In the typical case where the signal sequence length is much larger (say an order of magnitude or more), the steady state part of the output is also much longer than the transient part of the output.

Let's extend your example to illustrate this. If we were to apply your filter to a sequence of ones that contains $15$ samples instead of $5$, we would have an output sequence consisting of $15+5-1 = 19$ samples long. The first $4$ samples would be part of the transient and the last $4$ samples would be part of the transient (in fact, you can read this off of your $y$ output signal). If you were to drop these samples, you would see the steady state output which would be the value $3.72$ repeating $19-4-4 = 11$ times. Note that $3.72$ is the DC gain of your filter (the sum of the coefficients), so this is the expected outcome.

On the other hand, if we leave the first and last $N-1$ samples attached, then you would have
$$y=[0.5 \, 1.36\, 2.36\, 3.22\, 3.72\,
...\, 3.72\, ...\,
3.72\, 3.22\, 2.36\, 1.36\, 0.5]$$
where you can easily see the transient parts at the beginning and end of the filtered sequence.


What is the name of these simple filter algorithms? I found some snippets of code on the internet that are really simple signal filters. Unfortunately there wasn't much info to back it up. Using Google I've found similar code, but none of it backed up by research or documentation. I've tried it and it works well enough for my situation.

The low pass filter looks like this:

float RC = 1.0/(CUTOFF*2*3.14);
float dt = 1.0/SAMPLE_RATE;
float alpha = dt/(RC+dt);
float out[numSamples];
out[0] = in[0];
for(i=1; i<numSamples; i++){
    out[i] = out[i-1] + (alpha*(in[i] - out[i-1]));
}

And the high pass filter looks like this:

float RC = 1.0/(CUTOFF*2*3.14);
float dt = 1.0/SAMPLE_RATE;
float alpha = RC/(RC + dt);
float out[numSamples];
out[0] = in[0];
for (i = 1; i<numSamples; i++){
    out[i] = alpha * (out[i-1] + in[i] - in[i-1]);
}

I'd like to know if there is a name for these algorithms or who invented them so I can find more information about them.

out[i] = out[i-1] + (alpha*(in[i] - out[i-1]));

From Wikipedia's article on Low pass filters:

  That is, this discrete-time implementation of a simple RC low-pass filter is the exponentially weighted moving average

  […]

  The loop that calculates each of the n outputs can be refactored into the equivalent:

for i from 1 to n
     y[i] := y[i-1] + α * (x[i] - y[i-1])

So, yes, this is a "exponentially weighted moving average". You can call it that, a single-pole IIR low pass with a decay of $\alpha$, or just a single-tap recursive low pass.

I'm happy I could help you with this, it wasn't very hard!


From MFCC to Machine learning. What are the steps? I was able to get a dataset with MFCC coefficients. However, depending on the length of my sound file I get a different sized matrix. As in, 13 (13 MFCC coefficients) by XXX, where XXX will vary depending on the length of the sound file.
Does it make sense to 'normalize it' to keep XXX consistent? If so, how?
Like in this example, the size of the matrix always varies: https://archive.ics.uci.edu/ml/machine-learning-databases/00195/Test_Arabic_Digit.txt

Also, how will I feed this into a Machine Learning algorithm? (i.e. k-NN, HMM, etc.) I somewhat figured how to do it if it's just one line for each sound file (e.g. 1 by 13 for each sound). What are the steps if it's MFCC? I am a little lost here.

Thank you for your help.

There are two general things you can do with a machine learning algorithm:
1. Regression-- take input data and return an output numerical value. This is also known as curve-fitting, parameter estimation, etc.
2. Classification-- take input data and return a decision from a discrete set. In traditional statistical terminology, this is similar to hypothesis testing.

(I'm assuming that you are interested in classification.)

If you are doing supervised learning, you will need some training data, i.e. sets of MFCC coefficients that have been labeled with a priori knowledge. E.g. if you are identifying American English accents, you will need some audio data from New Englanders, Southerners, Midwesterners, etc., from which you extract the MFCC coefficients and label them manually. These coeffcients are known as features and the algorithm that distills down the high-dimensional dataset (i.e. the audio clip) down to a few coefficients (i.e. the MFCC coefficients) is known as feature extraction.

Once you have your training data and labels, you can apply a simple algorithm like KNN quite easily using libraries such as the Python sklearn library.


What is $H_2$ and $H_{\infty}$ control? I can create a Linear Quadratic Gaussian Integral(LQGI) controller very easy by using GNU Octave. LQGI is in the area of Optimal Control Theory.

But there is something called Robust Control Theory. And every time I'm searching for that I got these in my  search field:

$H_2$
$H_{\infty}$

Questions:

What are they?
How do I apply them to control the processes?
Are they good for control systems?
Are they better that PID and LQGI?
Are they easy?
Are they for nonlinear system?
Can they have a Extended Kalman Filter?
What's the difference between them both?
Can you show me an easy example how to apply them?
Other?

They are system norms, a metric that you can compare two different systems in terms of their generalized gain and spread. You can look these up no need for attaching physical motivation.
You don't. They are for assesing and using in the optimization programs.
Yes, but only if you know what you are doing.
Can be. But again, they are for assessment and optimization. If you hit the minimum with a PID controller then it won't make any difference. But hitting the minimum is pretty much unlikely.
Unfortunately no.
see 5 (you can pull some tricks by bounding the nonlinearity with a linear system etc. or for very limited nonlinear operators you can apply them)
Not applicable. They are different things.
They are different norms. But both of them are norms. Read about them.
See the documentation of hinfsyn of matlab
Yes


Does every 8 by 8 DCT matrix have a valid (in range) IDCT? For a given quantization level of JPEG compression, are all possible DCT matrices have valid IDCTs. By valid, I mean matrices whose elements are in the range [0,255] ?

All images blocks that are applied by DCT matrix have a valid IDCT - i.e. they can always bring back original pixels and in general inverse transfer is theoretically as well as computationally viable.
However, while your pixels have values in range 0-255 - the DCT matrix of the block never results in values which are confined between 0-255. Not only that, they are not even integers - they are real values. So while they are valid numbers, they can't be stored in the same container as pixels.
Given that they are real numbers, it would be useless to be used for compression against integers of the original source. So hence they are confined in a particular range and then divided by a specific number - this is a process of quantization. Hence, what really goes inside your JPG or MPG are quantized co-eofficients.
The quantization is not equal. Certain components are preserved more than others (with finer quantization) - hence dynamic range of each component, though finite, is not 0-255 but different for each DCT co-efficient. this is practically taken care by Huffman coding which maps these values into binary bits.
So they are valid in 'general' definition but not what you are interpreting as.


Carrier and timing offset Assume my reference is 1 ppm so if my transmit frequency is 400 MHz and
symbolrate is 100 kHz.

Then, my carrier will have +/- 400 Hz offset  and then my symbol will have +/-0.001 Hz offset. Similarly, the reciever carrier has +/- 300 Hz offset
and symbol  has +/- 0.002 Hz offset.

So is the carrier recovery is to track 400 + 300 Hz
and timing recovery is 0.001 +0.002 Hz ?

Answer
If you need to guarantee operation within these ranges, then yes. (Assuming the 0.001 Hz is a typo and meant 0.1 Hz.)
Stochastic Modelling
That "1ppm" is typically a "condensed" number, e.g. the observed maximum over 24h of oscillations, or can be found as integral of phase noise power, which is often given in an oscillator's datasheet. Typically, small frequency errors are more likely than large ones.
It might make sense to model your system probabilistically: Your frequency errors in receiver and transmitter, $X_{RX}$ and $X_{TX}$ are independent and will have a probability density function $f_{X_{RX}}(x)$ and $f_{X_{TX}}(x)$, respectively.
So, the probability density function for the difference $\Delta=X_{TX}-X_{RX}$ is given by $f_\Delta(\delta) = f_{X_{TX}}(x)*f_{X_{TX}}(-x)$, $*$ being the convolution operator.
This will allow you to define something like "in x % of cases, it should work, so I only need to account for a range of at most this and that in frequency error".
For example: This makes a lot of sense, since frequency synchronization is not the only thing that can fail in a receiver; if your channel has a 0.1% chance of simply being too noisy, then building a PLL loop filter that is wide enough to work for 99.99996% of cases doesn't make much sense, and you should probably rather reduce your loop filter bandwidth to make the receiver's phase recovery less sensitive to noise.


Difference equation with complex zero Let's assume I have the following transfer function:

$$
H(z)=\frac{z-\left(\frac{1}{\sqrt{2}}+i \cdot \frac{1}{\sqrt{2}}\right)}{z}
$$

It looks like a first order highpass-filter with a complex zero at $\frac{1}{\sqrt{2}} + i \cdot \frac{1}{\sqrt{2}}$.

If I do the inverse $z$-transform, I get the following difference equation:

$$
y(k)=x(k) - \frac{1}{\sqrt{2}} \cdot x(k-1) - i \cdot \frac{1}{\sqrt{2}} \cdot x(k-1)
$$
with $k$ as time step variable. My input signal is real valued and my output should be real valued as well.

How do I implement that filter e.g. in C?
How do I deal with that complex numbers in that case?

What you actually want to do is the following: you need a zero at $\omega=\pi/4$, but you also need to place a zero at $\omega=-\pi/4$ in order to get real-valued filter coefficients. So your transfer function becomes

$$\begin{align}H(z)&=\frac{(z-e^{j\pi/4})(z-e^{-j\pi /4})}{z^2}\\&=\frac{z^2-2z\cos(\pi/4)+1}{z^2}\\&=1-\sqrt{2}z^{-1}+z^{-2}\tag{1}\end{align}$$

The corresponding difference equation is

$$y[n]=x[n]-\sqrt{2}x[n-1]+x[n-2]\tag{2}$$


Meaning of Hop-Size in Filter Bank Interpretation of Short Term Fourier Transform The Short Term Fourier Transform (STFT) is used to obtain the time-varying spectrum from a signal. There are two ways to understand it - the overlap add interpretation, and the filter bank interpretation.
The overlap-add method is where the signal is cut up into windowed segments in time and the DFT of each segment is taken. The window is advanced forward by R samples in time, where R is called the hop-size.

The filter bank interpretation is where the window is seen as a low-pass filter, and the segment is shifted in frequency by an amount equal to the time segment. The hop-size here is interpreted as a down-sampling in frequency.

I am unable to understand this interpretation and would like to know how R is interpreted as down-sampling of the signal in frequency. This would also clarify the constraints on choosing a hop-size R, to prevent aliasing in frequency (as having a hop-size R means down-sampling in frequency).

Your message was apparently unfinished ('When reading about the DCT...'). I thus will be starting from it.

The JPEG compression format is an example of hopping. A DCT-filter bank (Discrete Cosine Transform) is applied on the whole image, and subsampled by $8\times 8$ (in both directions). Another interpretation is via a rectangular window, implemented in a fake overlap-add version, as the overlap is zero ($8\times 8$ blocks are almost disjoint). The DCT has a very simple interpretation, as it consists of an infinite continuous cosine basis, discretized and windowed with a rectangular window.

The traditional Short Term Fourier Transform (STFT) is an improvement, with the use of an overlapping window. In the continuous setting, an STFT with a window $h$ can be inverted with a lot of dual windows $g$ (under very mild assumptions on $h$ and $g$).

Going to discrete implementations, certain choices of overlap, hop  and window can lead to easy inversion. Inversion, or synthesis is the turning point:  you can "analyse" any signal with any window, or hop. But not all designs allow synthesis or inversion, although it should be permitted more easily than in the  critical and orthogonal framework, because of the potential redundancy.

There are instance of windows allowing discrete, non-redundant inversion, like the Modified DCT, the LOT (lapped orthogonal transform), the MLT (modulated lapped  transform), and a handful of others (GenLOT, GULLOT). They are used for instance in the mp3 (mpeg-1 layer 3) style of audio compression.

In the above case, the hop is equal to the number of channels, similarly to the case of the JPEG DCT: $8$-channel DCT, and a jump of $8$ pixels in each direction. Allowing only a window for all the channels is somewhat restrictive. The most efficient framework is that of oversampled filter banks:

There is no direct concept of a global window for all $M$-band filters, $M\ge N$. Instead, some windowing is embedded in each one. In the above picture, the $N$ corresponds to your $R$. My understanding is that the downsampling takes place in the time-domain (hopping), but it results in a shift in frequency (and aliasing), as filter $H_i$ outputs are sent to the base-band.

Then, given a set of analysis filters $H_i$ (including some overlapping), and if $M$ is greater than $N$, you can very often cancel all the aliasing induced by the $N$ downsampling with an infinity of inverse filter banks. One issue is to find ONE suitable inverse.

This is discussed for instance in  Optimization of Synthesis Oversampled Complex Filter Banks, J. Gauthier et al., IEEE Transactions on Signal Processing, 2009:

  An important issue with oversampled FIR analysis filter banks (FBs) is
  to determine inverse synthesis FBs, when they exist. Given any complex
  oversampled FIR analysis FB, we first provide an algorithm to
  determine whether there exists an inverse FIR synthesis system. We
  also provide a method to ensure the Hermitian symmetry property on the
  synthesis side, which is serviceable to processing real-valued
  signals. As an invertible analysis scheme corresponds to a redundant
  decomposition, there is no unique inverse FB. Given a particular
  solution, we parameterize the whole family of inverses through a null
  space projection. The resulting reduced parameter set simplifies
  design procedures, since the perfect reconstruction constrained
  optimization problem is recast as an unconstrained optimization
  problem. The design of optimized synthesis FBs based on time or
  frequency localization criteria is then investigated, using a simple
  yet efficient gradient algorithm.


Vehicle counting from a video Disclaimer: I posted this question on StackOverFlow but learned that it is better suited here.
I am a bit new to image processing so I'd like to ask you about finding the optimal solution for my problem, not help for code.
I couldn't think of a good idea yet so wanted to ask your advice.
I'm working on a project using OpenCV on counting the vehicles from a video file or a live camera. Other people working on such a project generally track the moving objects then count them but instead of that, I wanted to work with a different viewpoint; asking user to set a ROI (Region of interest) on the video window and work only for this region (for some reasons, like to not deal with the whole frame and some performance increase), as seen below. (The user can set more than one ROI and the user is asked to set the height of the ROI about two times of a normal car by sense of proportion.)

I've done some basic progress so far, like background updating, morphological filters, thresholding and getting the moving object as a binary image something like below.

After doing them, I tried to count the white pixels of the final thresholded foreground frame and estimate whether it was a car or not by checking the total white pixels number (I set a lower bound by a static calculation by knowing the height of ROI). To illustrate, I drew a sample graphic:

As you can see, it was easy to calculate the white pixels and check if it draws a curve by the time and determining whether a car or something like noise.
Problem
I was quite successful until two cars passed through my ROI together at the same time. My algorithm incorrectly counted them as one car. I tried different approaches for this problem and similar to this like long vehicles but I haven't found a satisfactory solution.
Question
Is it possible to handle this task by this approach of pixel value counting? If not, what do you suggest?

You can keep track of the last 10-15 frames and count the number of white pixel of these frames. You can differentiate if it is a vehicle or not, since the graph will be high for 10-15 frames you will be able to say that it is a vehicle otherwise it is noise.

To count multiple cars in multiple side by side lanes you can check the pixels column wise, I mean vertically. For the time being you must be checking the pixels horizontally for white pixels.


Decompose stereo signal Basically the same as this question, but it was never really answered: Stereo sound separation into three sources

Imagine I have a really trivial signal:

Left: 100Hz tone + 200Hz tone + 500Hz tone
Right: 500Hz tone + 1kHz tone + 2kHz tone

I'd like to perform some operation where I can decompose Left and Right into L', R' and C. Where:

L': 100Hz + 200Hz
R': 1kHz + 2kHz
C : 500Hz

I'd like the algorithm to work reasonably well for mixed, commercial music, which would exclude just a simple FFT analysis (I think.)

I started out looking at correlation, but I don't think that applies in this case. Now I'm thinking I need to be looking at ICA, but that seems to be focused on producing N signals from a single input, and not 3 signals from 2 inputs...

There a few commercial algorithms that do exactly that (Dolby Prologic, DTS Neo 6, Lexicon Logic 7, Bose Videostage, etc.). If your project has the funds, you can simply try to license one of those.

The inner workings of these algorithms are rather complicated and typically a mixture of time domain and frequency domain feature extraction, some steering logic, and some real time steering mechanics (time variant mixers, sub band filters, etc.)

Simple signals like your example are normally easy to tease apart, however almost all music is  partially correlated and overlapping both in time and frequency. For starters it's difficult to define what actually SHOULD be in the center and the Left and Right. The algorithm has to make a lot of "artistic" decision to make sure there are no steering artifacts and it meets the design intent over a wide variety of source material.


How to evaluate an impulse/response model in MATLAB This is my first experience in signal processing, so don't treat me badly :)
However, I have estimated an impulse-response model using the MATLAB System Identification toolbox. The output, validated on unseen data, has a good fit, approximately around 80%. However, when I export the object imp - its type is idtf - containing the estimated model, I don't know how to evaluate it on unseen data. Obviously, in the structure is contained the polynomial structure of the model, but it bothers me to write a function to estimate it... Is there any quicker procedure?

One way of evaluating the model response to arbitrary data is the lsim command. As you don't provide more information with respect to your data and/or model an example below (taken straight from the MATLAB help). See help lsim or doc lsim for more information.

% create model
%
num = [1 4];
den = [1 20 5];
G = idtf(num,den);

% evaluate response
%
t = 0:0.01:5;
u = sin(t);
y = lsim(G,u,t);

% plot
%
figure;
plot(t,y)

In addition, if desired, you can evaluate the model impulse, step and/or frequency response using the impulse, step and bode commands.


How to chose a kernel function for a FIR low pass filter in audio processing? As the answer of a previous question, I was suggested to use a FIR filter to act as a low-pass filter to reconstruct the envelope of a sound wave (edited, emphasis mine):

  As I said in a comment, you can get the envelope of a signal by
  running it through a lowpass filter

  [...]

  Implement a lowpass filter (FIR) by creating a filter kernel of appropriate length M (h(M) ). Note that your FIR filter
  (sin(x)/x) should normally be multiplied by a window function (e.g
  Hamming, Blackman, etc).
  Convolve your signal with the filter kernel.

  [...]

I'm having a hard time choosing/finding a proper filter kernel function. The quote above suggests sin(x)/x. Is this a standard kernel function? Does it have some magic properties that make it particularly suitable? What should I check if I consider another kernel function?

Yes the sinc kernel $$h[n] = \frac{\sin(\omega_c n)}{\pi n}$$ is the theoretical result for the ideal lowpass brickwall filter whose gain is one and cutoff frequency is $\omega_c$ radian per samples.

Your filter impulse response will be samples of that function. Truncate it symmetrically to some finite length, and then apply a window on those samples. Windowing will help you achieve a much smoother response.


5G NR Network Collision Avoidance I'm looking for a technical explanation or a specification that explains exactly how 5G NR collision avoidance will work. I was watching a talk from the Mobile World Conference in Shanghai two years ago (2018) and it appears that centralized transmission scheduling may be a part of the standard - here is a timestamped link to the panel member that references this:

Steve Greaves addresses 802.11ay collision avoidance

To quote:

  If you look at what Facebook are doing with AY they are trying to get coordination & synchronization across the
  whole network. Once you have that capability what you can do is you
  can manage interference. You can control the scheduling of
  transmissions such that you can deal with the dense urban deployments
  which are interference limited.

I can't find a standard or explanation for this exact thing within 5G so I was hoping someone could point me in the right direction.

Note: I am not intimately familiar with 5G and am not sure what parts of the "standard" for 5G have been defined. For a long time there was no standard at all. Also please let me know if there is a better community for this question.

IEEE 802.11 is not related to 5G NR, so consideration for WLAN networks following 802.11ay don't apply to cellular infrastructure.

I don't know what "ad-hoc" modi 5G NR brings; I don't think it's any.

So, 5G doesn't do any collision avoidance, because it's not a problem it has: it's a cellular standard that defines that the basestation / cellular infrastructure defines who sends when, where. The mobile network operators own the spectrum, so there's by definition no uncooperative users.

Generally, whereas interference mitigation is in fact a non-negligible topic for coordinated microwave networks, with these deci- to centimeter waves still taking many non-optical paths and receiver sensitivity being high enough for that to create "unwantedly far reach", you can basically ignore that for cellular mmWave applications, where you can be pretty sure that you, as the owner of the used spectrum, will have a relatively easy time making sure that your own femtocell base stations don't lead to mutual interference: The propagation is pretty much optical, so, the potential numbers of base stations a single handset might reach is pretty limited by lines of sight and a pretty small radius.


Golay complementary sequences for channel estimation - MATLAB There are set of complementary sequences known as Golay sequences that are used for channel estimation because of the nice property they have which is that the sum of the autocorrelation function of each gives dirac function.

For example, these are two complementary sequences

$$Ga=[+1, +1, +1, −1, +1, +1, −1, +1]$$
$$Gb=[+1, +1, +1, −1, −1, −1, +1, −1]$$

$$R_a(i)+R_b(i) = \delta(i)$$

Due to this nice correlation property, they are used in digital communication for as known pilot sequences for channel estimation of wireless channel. For example assume the channel we want to estimate is $h(t)$ then sending $Ga$ and $Gb$ over channel first and then performing correlation at receiver with Ga and Gb provides a good estimate of channel.

$$(h(t)*Ga)*Ga + (h(t)*Gb)*Gb =h(t)*(R_a+R_b)= h(t) \delta(t)$$

In theory, this is easy to understand. However I am trying to implement a more realistic setup in MATLAB.

If $G_a$ and $G_b$ are 10 chip in length with have a chip time of $T_c=0.57 ns$ and the channel is a 2 tap with delay spread of the channel is $T_{delay}=2.7 ns$. And assume we have Additive White Gaussian noise. My aim to estimate the channel. Below is my code.

channel =[1 0.2]';
Ga=[+1, +1, +1, -1, +1, +1, -1, +1]';
Gb=[+1, +1, +1, -1, -1, -1, +1, -1]';
y1=conv(channel,Ga);
y2=conv(channel,Gb);
r1=awgn(y1,10);
y2=awgn(y2,10);
hac1 = dsp.Autocorrelator;
correlation_output1=step(hac1,Ga);
correlation_output2=step(hac1,Gb);
channel_estimation=correlation_output1+ correlation_output2;

As you can tell, I have not included the delay spread $T_{delay}$ neither the chip rate $T_c$. Can anyone help me understand how I can incorporate such in my MATLAB code?

I would try with this code:

channel =[1 0.2]';
Ga=[+1, +1, +1, -1, +1, +1, -1, +1]';
Gb=[+1, +1, +1, -1, -1, -1, +1, -1]';
y1=conv(channel,Ga);
y2=conv(channel,Gb);
correlation_output1=xcorr(y1,Ga);
correlation_output2=xcorr(y2,Gb);
channel_estimation=correlation_output1+ correlation_output2;
stem(channel_estimation./length(Gb))

Best regards,


What Is a Weighted Local Histogram? I'm reading through a couple of academic papers, and this terms often comes up "local weighted histogram". An example of quote is the following:

  First, cumulative histograms are built for every pixel from its neighborhood, using Gaussian-neighborhood weighting

I would be able to construct a local histogram, I'm not sure what this weighting is about.

Say the central pixel is at $p_0 = (y_0,x_0)$, and we have two pixels of same gray scale value $r$ at coordinates $p_1 = (y_1,x_1)$ and $p_2 = (y_1,x_1)$, but such that $d(p_0,p_1) < d(p_0,p_2)$. A normal histogram would count that value $r$ twice. With a Gaussian weighting for example how would we count such pixel value in the histogram construction?

A weighted local histogram would mean filtering the image $I(x,y)$ with a localised filter $H(x,y)$ (gaussian in this example). The resultant image is the 2D convolution $Y(x,y) = I(x,y)*H(x,y)$. The normal histogram of $Y(x,y)$ is the localised histogram of $I(x,y)$. It is used to view certain specific characteristics of an image depending on the filter used


What is the Technique to Find Variance of Estimation Error Given an $n$-vector $y$ (responses) and a design matrix $X$, I wish to fit them with a simple linear regression model
$$y=X\beta+e,$$ or,

$y_t = x_t'\beta_0 + e_t$
where $e\sim\mathcal{N}(0, \sigma^2I)$ and $\beta_0$ is the true parameter. Then, we have
$$y\sim\mathcal{N}(X\beta, \sigma^2I).$$
Then the maximum likelihood estimations (MLE) of $\beta$ which is the Ordinary Least Square estimator is
$$\hat\beta_n=(X_n^TX)^{-1}X_n^TY_n$$. Then how to calculate the Expectation of the variance of the estimation error viz. $E[(\hat{\beta_n} - \beta_0){(\hat{\beta_n} - \beta_0)}^T]$ . I am struggling to understand what to substitute for $\beta_0$ in this calculation and how to find the expression. Any help will be extremely useful. Thank you.

I don't understand the subscript $n$ notation, however, in the least squares problem that is given by:
\begin{equation}
{\bf{y}}={\bf{H}}{\theta}+\bf{n},
\end{equation}

where ${\bf{n}}\sim\mathcal{N}(\bf{0}, \sigma^2I_N)$ is a zero mean additive white Gaussian noise and $I_N$ is the $N \times N$ identity matrix, the maximum likelihood and the least squares estimators are equivalent and are efficient estimators (i.e. unbiased and attains the Cramer Rao bound) and are given by
\begin{equation}
\widehat \theta_{ML}=\left ( \bf{H^TH} \right )^{-1} \bf{H^Ty} .
\end{equation}

Now, let's calculate the first two moments of the estimator:

\begin{equation}
E\left \{ \widehat \theta_{ML} \right \} = E\left \{ \left ( \bf{H^TH} \right )^{-1} \bf{H^Ty} \right \} = \\
\left ( \bf{H^TH} \right )^{-1} \bf{H^T} \underset{=\bf{H}\theta}{\underbrace{E\left \{  \bf{y}\right \}}} = \theta,
\end{equation}
where $\theta$ is the true parameter value.
Thus, the estimator is unbiased and its Mean Square Error(MSE) will equal its covariance.

Now for the covariance:

Remember that if $\bf{z}=\bf{Ax}$ then $cov\left ( \bf{z} \right ) = \bf{A} cov\left ( x \right ) \bf{A^T}$.
Since the mean doesn't change the covariance in our problem $\bf{y} \sim \mathcal{N}\left ( \bf{H\theta}, \sigma^2I_N \right )$, thus,

\begin{equation}
cov\left ( \widehat \theta_{ML} \right )=\left ( \bf{H^TH} \right )^{-1} \bf{H^T} cov\left ( \bf{n} \right ) \bf{H}\left ( \bf{H^TH} \right )^{-1} = \\
= \left ( \bf{H^TH} \right )^{-1} \bf{H^T} \sigma^2 \bf{I_N} \bf{H}\left ( \bf{H^TH} \right )^{-1} = \\
= \sigma^2 \left ( \bf{H^TH} \right )^{-1} \bf{H^T}  \bf{H}\left ( \bf{H^TH} \right )^{-1} = \sigma^2 \left ( \bf{H^TH} \right )^{-1}.
\end{equation}

Note that $\bf{H^TH}$ is symmetric so it equals its transpose. Same is correct for its inverse.

I believe that this answers your question completely.


How to check if the poles of transfer function are stable? $$H(z)=\frac{az}{[z-(b+ja)][z-(b-ja)]}$$

The two poles are looking like this:
$z_{\inf, 1}=b+ja, z_{\inf, 2}=b-ja$

I know they must be inside the unit circle: $|z_{\inf,i}|<1$

So I replace the term: $|b+ja|<1$ by the first pole, but how I must continue?

It's not $|z_{\inf,i}<1|$ it's $|z_{\inf,i}|<1 $ which means that, as @Mohammad Mohammadi said in the comments, your pole on a complex plane must be inside the unit circle.

Analytically, you get the magnitude with
$$ r = |z| = \sqrt{a^2+b^2} $$
If $r<1$, your pole is stable.

Note that $z_1 = b + ja$ and $z_2 = b-ja$ have the same magnitude, which means that if one is stable, the other is stable as well (in the complex plane, they are symmetrical to the horizontal axis).

Illustration

Matlab code for the illustration

%%
clear
close all
clc;

% STABLE

a = 0.73;
b = 0.24;

P1 = a + b*1i; % Pole 1
P2 = a - b*1i; % Pole 2

r = sqrt(a^2+b^2); % Magnitude

% UNSTABLE

a2 = 1.23;
b2 = 0.84;

r2 = sqrt(a2^2+b2^2);

P12 = a2 + b2*1i; % Pole 1
P22 = a2 - b2*1i; % Pole 2

zplane([0 0; 0 0], [P1 P12; P2 P22]); % Plot the poles on a complex plane
legend('', '', ...
    ['z = ' num2str(a) ' \pm ' num2str(b) 'i (r = ' num2str(r) ')'],...
    ['z = ' num2str(a2) ' \pm ' num2str(b2) 'i (r = ' num2str(r2) ')'], ...
    'Unit Circle (r=1)')


The physical meaning of comparing the phases of two signals with *different frequency*? I was examining my simulation data about vibration. I was trying to compare the displacement phase difference between two cylinders; the cylinders carry out nearly sinusoidal vibrations. However, I just figured out that at certain conditions, the dominant displacement frequencies of two cylinders are different (i.e. C2 has dominant frequencies at both equivalent and two times of the C1's frequency), although in most cases, their frequencies are essentially the same.

So my question is: is the phase difference (e.g. in phase, in anti-phase) still physically meaningful when the frequencies are different?

Appreciate any explanation. (I think some physical examples would be very helpful)

Edit: found a discussion here, but still trying to understand.

There is not a phase in the usual sense in which there is a number $\phi_0$ such that

$$x_1(t)=x_2(t-\phi_0) \ \forall t$$

where $x_1(t)$ and $x_2(t)$ are two periodic signals of different frequencies.

However, there are some phase-like analysis one could perform, depending on the application. I'll talk about two examples that come across my mind quickly.

Suppose we have the two signals $x_1(t)=\sin(t)$ and $x_2(t)=\sin(3t-\phi_0)$. As the latter is the 3rd harmonic of the former, their sum is periodic. However, the shape of the resulting wave would vary a lot depending on $\phi_0$ (some kind of "phase difference" between the signals). In the image below, I plotted three different waves ($x_1(t)+x_2(t)$ for three different values of $\phi_0$):

$\phi_0=0$ in red.
$\phi_0=\pi$ in blue.
$\phi_0=\pi/2$ in green.

As you can see, there is a strong dependence on the value of $\phi_0$. This would be very important in music, for example, as the resulting wave would have a pretty different sound depending on the "phase" between the harmonics.

Following the same example, there is another way to see the phase between these two signals. Take for example the case in which $\phi_0 = 0$:

As you can see, whenever the fundamental wave crosses the $t$-axis, its 3rd harmonic also does. On the other hand, whenever the first one reaches a maximum (a minimum), its 3rd harmonic reaches a minimum (a maximum). This could lead to one thinking that they are in anti-phase. Again, this depends on what your application is and what you are looking for doing this analysis, as well as on what you define as "phase".

As a sidenote, there is a typical example that uses the phase between two waves of different frequencies (or something like that). It's called PLL (or phase-locked loop). You can check that out if you are interested in the topic.


Connecting segmented region in image I am working on segmenting hands from monocular images and I have already come up with a pretty descent segmentation. The only issues I am having are when people have rings on their finger. Then the hand is no longer a connected component and I get the hand and the finger component. The following image shows my segmented mask, where the ring finger is detached because of a ring.

Is there any natural method to connect components like these together? Otherwise I could try to segment the rings, but not all people have rings, so I would prefer to try to tackle this in this way.

Also, if there is another SE better suited for this question, please let me know.

I guess that the missing area is due to a ring. A solution would be to use a PCA in order to compute the finger orientation. Then you apply a closing (mathematica morphology) with an segment oriented segment (same orientation as the finger).

It might modify a little bit the gap between the fingers where they are connected to the hand.


How to initialize an IIR filter in Matlab I designed an IIR filter and want to do zero, step and projection initialization.  There is mention in different reference papers but no explicit examples for Matlab.  How do you do IIR initialization, for the previously mentioned types, in Matlab?

Not easy. The Matlab filter() function implements a general purpose IIR filter and allows you to pass in an initial state. However, they don't publish the exact definition of the state and the algorithm deployed. It seems to be close to a transposed form II filter but with some subtle numerical differences, so the state doesn't quite match.

In any case it would be better to break this down into second order sections using sosfilter(). However that doesn't allow you to set an initial state. You may have to write your own filter function for that.

YOu still can use filter() to some extent. The default is zero initialization. You can also do step initialization by simply running a constant input for a while and collecting the state as an output argument from the filter, like this:

[b a] = butter(4,1000*2/44100);
% create state for a constant input of "1"
[~,state1] = filter(b,a,ones(1000,1));
% now filter for a constant of 2 but with a state for "1" input
ystep = filter(b,a,2*ones(100,1),state1);
plot(ystep);

Not sure about projection initialization, though.


Is there any trade-off between the sparseness and residual of the constructed data? I am going to test a sparsity promoting algorithm with different parameters to solve the $Ax=b$ to obtain a sparse model $x$. I have to test different parameters for which my algorithm results in different answers for $x$.

With different $x$'s there are different $Ax=\hat{b}$ 's. To have a better comparison, I forced my algorithm to iterate while the $\left\|Ax-b\right\|^2<\epsilon$.

Now that I obtained all $x$'s for which $\left\|Ax-b\right\|^2<\epsilon$ , how can I determine which model, $x$, is better? I know that I am looking for a model that is more sparse and have a better fit to $b$ but I think I have to consider somehow a trade-off between the sparseness and residual.

You can reformulate your problems as a minimization problem:

$$P1:  \min_x  \|Ax -b\|_2^2 + \lambda \|x\|_o $$

Where $\|x\|_o$ is the $\ell_o$ norm of x where it is the number of nonzero elements. Looking for a solution to p1 is essentially looking for a solution that best explains the model $\|Ax -b\|$ while having sparse coefficients. The problem is NP non convex NP hard and is extremely difficult to optimize. Another approach is to use the best convex envelope to the $\ell_o$ norm which is the $\ell_1$ norm.

$$P2:  \min_x  \|Ax -b\|_2^2 + \lambda \|x\|_1 $$

Under some conditions on A, the two problems are in fact equivalent in the sense they have the same global optimal solution. There are many efficient solvers for P2 including ADMMS, randomized methods, ISTA and FISTA.

To solve P2 in MATLAB, SPAMS is your friend here.


Audio template matching - power correlation? I'm looking at some code that matches audio templates in a longer audio file. The calculation correlates the power spectra of the template and audio file, maximizing over the possible alignments. This seems sub-optimal to me, because by going over to the power spectrum we're throwing away phase information. Yet when I try doing a full correlation I get worse results. I'm not 100% sure yet my implementation is correct; I could have made a silly mistake.

Is there any reason it could be better correlating the power spectra than doing a full correlation? E.g. is it more noise resistant? (I don't see why it would be. As it happens, I do have something resembling white noise in my test data.)

The only obvious thing I can think is that the power spectral correlation is probably better if the alignment error is >= 1/highest frequencies in signal because then the template and signal will be out of phase. I don't think this should be the case for me: I'm optimizing to within more temporal precision than that.

Other ideas?

EDIT: in view of Peter K's comment, I should clarify that I'm using the short time Fourier transform, summing over windows of size around 0.01s. That's how the alignment dependence enters.

If you're cross correlating short time power spectra and disregard phase, then you're not really disregarding phase.

Phase of the global Fourier transform encodes the temporal structure, including the position of the transients, evolution of tones or the local incoherence of noise.

By capturing the time dependence in the moving window of the time-frequency power spectrum, the really important aspects of phase are represented in the temporal evolution of the frequency power density.

So what you are discarding is merely "local phase", which contains a lot of information that is however modified by even the most subtle processes like sound propagation, speaker reproduction, microphone recording, etc. These modifications don't affect the qualitative content of the sound a lot, and a robust sound recognition algorithm should be mostly insensitive to them.

That means discarding the local phase will make your recognition or correlation algorithm more robust and avoid misclassifications due to small inaudible phase errors, while at the same time preserving the total temporal structure of the signal.


Sinc Interpolation Using DFT (FFT) Lets say I want to double the number of points in an array f. I had the idea to do this:

F=fft(f);N=length(f);
FF=[F(1:N/2) zeros(1,N) F(N/2+1:N)];
f=ifft(FF);

But the result is not exactly correct. Why is this?

Doing DFT based interpolation has to keep 3 principles:

Keep the zeros outside the data (Looking form $ -\pi $ to $ \pi $). Namely add zeros to the higher frequencies.
If it is real signal, keep the symmetry (If needed, divide the highest frequency into 2 items).
Preserve the energy (By the division by 2 if needed).

By the way, the interpolation isn't by a Sinc Kernel, it is done with its discrete form - The Dirichlet Kernel.
See The Proper Way to Do Sinc Downsampling (DFT Downsampling) for Uniformly Sampled Discrete Signals with Finite Number of Samples


What is the difference between the Point Spread Function and Sampling Aperture? I've been told that the point spread function of a pixel is its distribution of intensity, while the sampling aperture of a pixel is its distribution of sensitivity. I'm a little unsure of what the difference between intensity and sensitivity is. It sounds like the former deals with how a pixel is displayed, while the latter deals with how it is captured. Is that correct?

Sounds correct to me. The point-spread function indicates how the intensity of a point source is spread among multiple pixels in an image. In an imaging sensor, the sampling aperture of a particular pixel indicates how it responds to intensity from points in space as a function of the spatial offset from the pixel's center.


Automated accent detection How can I compare the way that word is pronounced with a sample? For example, most of dictionaries have sample word pronunciation recorded. The goal is to provide feedback to the user so that he or she can improve pronunciation.

I wonder if there are API available for such task or existing programs that are able to do it.

I found similar question here How do I compare 2 audio files, returning a percentage of difference , but I am providing a definition what actual difference is.

Regards,
Eugene

That was the same topic of my Graduation Project, I am not yet sure of the best way to do it as this matter is still a research matter, but what my team and I did was as follow:

1- Provide a database representing the features of the right pronunciation of each word. This database contains for each word around 10 different features for 10 different people reading this word.

2- Get the features of the input word ie. a .wav file representing the word

3- compare the features of the input word with the 10 features representing the right pronunciation of the word

To get the features, you will need to get the MFCC file a .wav file containing the word pronounced
The MFCC is a matrix of numbers, so to compare 2 MFCCs, you will need to use dynamic programming algorithm, in our case we've use Dynamic Time Wrapping algorithm.
You might also need to make a filtering on the input words to filter them from noise.
You might also need to check for the following programs:
Praat  >> to make a filtering if needed, and also can be used to make MFCC file
SPPAS  >> to get the feature of the .wav file (MFCC)
Matlab >> to compare 2 MFCC files

Hope this can help


Question of Digital filter; stability and order Below 1~6 are my question.

Shortcomings of IIR filter.
Mean of unstable. (Things that will happen when system is unstable.)
How to select the best order(n).
What happens when I raise order too much.
How can I know the maximum order or appropriate order.
How to determine stability of system in MATLAB.

I learn IIR filter's shortcoming is unstable. But I don't know mean of unstable. I want to know things that will happen when system is unstable and stable.
Additionally, I wonder the order of filter. I thought when I raise order, FIR and IIR filter will have more sharp transition zone(skirt?). It is true, then how much I raise the order? When I raise order too much, what happen will become? I used butterworth filter.

Some of your questions overlap so I’m not going to answer in the framework of 1 followed by 2, .....

An IIR filter ‘s stability is related to the location of its poles.  If the region of convergence contains the unit circle, the filter is stable. This definition covers both causal and anti causal filters, so you can have a stable causal filter and an unstable anti causal filter or vice versa.

Filters are implemented in terms of their coefficients so you need to root the denominator polynomial to find the poles and this is an ill conditioned problem.

see https://math.stackexchange.com/questions/70015/roots-of-bivariate-polynomials

for an example.  small perturbations in filter coefficients, like finite precision, can move the poles so that the region of convergence no longer contains the unit circle.

A good DSP book like Oppenhiem and Schaefer will discus work arounds like factorizing a polynomial into smaller stable polynomials and filter structures that are more robust to precision errors.

There are a number of functions in MATLAB that estimate oder for a given filter performance but these aren’t always reliable.

Ultimately one uses trial and error when designing a filter.  start with an oder and test the filter. if the filter meets spec, stop. if not, increase order by one and repeat.

How do you know if a filter is stable? you test it, which is what MATLAB is good at, although add ons like the fixed precision tool box might make things easier.

Overdetermined filters can be stable or unstable. An overdetermined filter can be a poor use of computation.

FIR filters are always stable. While an IIR filter requires fewer arithmetic operations than a FIR filter, some processors have execution pipelines that penalize recursion, so you can be doing more arithmetic faster with a FIR filter than the equivalent IIR filter.

A FIR filter can be linear phase, which can be desirable in many applications. An IIR filter can only be linear phase over a portion of bandwidth.

IIR filters can borrow from analog filters lke the Butterworth filter you mentioned

The answer to question 3 depends on what you mean by best.  It isn’t really an issue of order but how oder is structured


What is the difference between closed form and opened form in z transform I am new in DSP and I need to have a clear udrastanding about below example.
Consider the below $h[n]$ as the impulse response of a system:

\begin{equation}
h[n] = a^n(u[n] - u[n-N])
\end{equation}

so the opened z transform can be expressed as :
\begin{equation}
H(z) = 1 + a^1 z^{-1} + a^2 z^{-2} + \dots + a^{N-1} z^{-(N-1)}
\end{equation}
and its closed form is expressed like :
\begin{equation}
H(z) = \frac{1}{z^{N-1}} \frac{z^N-a^N}{z-a}
\end{equation}
defenetly, closed form is better than the opened form for hardware implementation because as below, the closed form can be written in the time domain like this:
\begin{equation}
y[n] = x[n] - a^N x[n-N] + a y[n-1]
\end{equation}

Now the only difference between the closed and opened form is a zero and a pole at
$z=a$. Actually, it seems to me, by adding zero and pole, we obtained a better representation. But my question is,

What is the drawback of adding zero and pole to the system at $z=a$.

How can we analyze the benefits and drawbacks between closed and opened form?

I would appreciate any help you can provide.

One possible difference occurs for a pole outside of unit circle. Note that the original system is a truncated IIR filter. As such, it's an FIR filter and have no stability issue for any finite $a$.

However, the closed form representation includes a pole, and implements a pole zero cancellation. Eventhough theoretically (with infinite precision arithmetic) this would not pose a problem, for finite length arithmetic and for large values of $N$, and for $|a| > 1$, you could end up in numerical issues.

Specifcally, the instantaneous dynamic range (the mantissa length of the floating point representation) would be large enough to be able to accurately represent the arithmetic difference indicated in the numerator. Indeed this could also be aproblem with the finite length implementation when $a^N$ is too large to be added to $ax[n-1]$, and that could happen when $|a|>1$. So which one would be more sensitive to this numerical issue? You better test by yourself. Furthermore, the closed-from version could impose an additional problem for those $ a \approx 1$ which would not be the case for the FIR filter.

As an advantage of te closed form, note that the FIR form and IIR form requires the same amount of memory (to hold past values of input up to a delay of $N$ samples), however the closed form (truncated IIR) requires much less computation MACS compared to FIR version for large N.


Response of a system to a step function (heaviside) I'd like to compute the response to a step function of a electrical/thermal system. Generally I can "easily" compute the transfer function $H$:

$$H(\omega) = \frac{V_{out}(\omega)}{V_{in}(\omega)}$$

Since the Fourier transform ($\mathcal{F}$) of the Heaviside function is (computed with WA):

$$\mathcal{F}(\theta(t)) = V_{in}(\omega) = \sqrt{\frac{\pi}{2}}\delta(\omega)+\frac{i}{\sqrt{2\pi}\omega}$$

Hence, noting $\mathcal{IF}$ the Inverse Fourier transform:

$$V_{out}(t) = \mathcal{IF} \left\{ \left( \sqrt{\frac{\pi}{2}}\delta(\omega)+\frac{i}{\sqrt{2\pi}\omega} \right) H(\omega) \right\}$$

To check my math I tried to compute the response for a simple RC system:

I should get the well known charge of the capacitor. The transfer function:

$$H(\omega) = \frac{1}{1+i\omega R C}$$

Computing the Inverse Fourier transform ($\mathcal{IF}$) with WA ($R=C=1$) I get:

This would be correct if we were going backward in time :/. So the question is... What am I doing wrong?

I did the same using Laplace Transforms and everything works fine... But I don't understand why.

P.S. I don't want another method, I just want to understand what's wrong in my approach.

P.S. the reason why I am using WA is that for my more complicated system I need to compute the Fourier transforms using WA.

Fat32's derivation of the result via the Fourier transform is correct, but I think that your original question hasn't really been answered ("what am I doing wrong?"). The real reason why WolframAlpha gives you a wrong result is the different definitions of the Fourier transform used by WolframAlpha and the one used by you to derive $H(\omega)$. By default, WolframAlpha uses this definition of the Fourier transform:

$$F(\omega)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}f(t)e^{j\omega t}dt\tag{1}$$

Apart from the factor $1/\sqrt{2\pi}$ (which makes the transform unitary), the actual problem lies in the sign of the exponent. You implicitly assume a negative sign, whereas WolramAlpha uses a positive sign.

There are two ways to solve the problem. Either you change Wolfram Alpha's default parameters using FourierParameters -> {1, -1} as explained here, or you change $H(\omega)$ accordingly, which means that you would need to use

$$H(\omega)=\frac{1}{1-j\omega RC}\tag{2}$$

As an example I tried the second method, and WolframAlpha indeed provides the correct result.

One more thing, for this kind of problems it's usually more convenient to use the Laplace transform because this avoids the use of generalized functions, such as the Dirac delta.


When was the time-frequency analysis invented? It seems there are several papers from the seventies but backtracking from the references gets quickly difficult. Who calculated for the first time a time-frequency representation of a signal?

According to the preface of Foundations of Time-Frequency Analysis, a rough timeline is as follows:

1930 - Early development of quantum mechanics by H. Weyl, E.Wigner, and J. von Neumann.
1946 - Theoretical foundation of information theory and signal analysis by D. Gabor (cf. "Theory of communication").
1980 - Time-frequency analysis established as an independent mathematical field (apart from engineering) by Guido Janssen.
1990 - Development of wavelet theory. Overview of the mutual influence given in "The wavelet transform, time-frequency localization and signal analysis" by Ingrid Daubechies.