Detect Closely Spaced Sinusoids

Open Live Script

Consider a sinusoid, $f (x) = e^{j 2 π ν x}$ , windowed with a Gaussian window, $g (t) = e^{- π t^{2}}$ . The short-time transform is

$V_{g} f (t, η) = e^{j 2 π ν t} \int_{- \infty}^{\infty} e^{- π (x - t)^{2}} e^{- j 2 π (x - t) (η - ν)} d x = e^{- π (η - ν)^{2}} e^{j 2 π ν t} .$

When viewed as a function of frequency, the transform combines a constant (in time) oscillation at $ν$ with Gaussian decay away from it. The synchrosqueezing estimate of the instantaneous frequency,

$Ω_{g} f (t, η) = \frac{1}{j 2 π} \frac{e^{- π (η - ν)^{2}} \frac{\partial}{\partial t} e^{j 2 π ν t}}{e^{- π (η - ν)^{2}} e^{j 2 π ν t}} = ν,$

equals the value obtained by using the standard definition, $(2 π)^{- 1} d \arg f (x) / d x$ . For a superposition of sinusoids,

$f (x) = \sum_{k = 1}^{K} A_{k} e^{j 2 π ν_{k} x},$

the short-time transform becomes

$V_{g} f (t, η) = \sum_{k = 1}^{K} A_{k} e^{- π (η - ν_{k})^{2}} e^{j 2 π ν_{k} t} .$

Create 1024 samples of a signal consisting of two sinusoids. One sinusoid has a normalized frequency of $ω_{0} = π / 5$ rad/sample. The other sinusoid has three times the frequency and three times the amplitude.

N = 1024; n = 0:N-1; w0 = pi/5; x = exp(1j*w0*n)+3*exp(1j*3*w0*n);

Compute the short-time Fourier transform of the signal. Use a 256-sample Gaussian window with $α = 20$ , 255 samples of overlap between adjoining sections, and 1024 DFT points. Plot the absolute value of the transform.

Nw = 256; nfft = 1024; alpha = 20; [s,w,t] = spectrogram(x,gausswin(Nw,alpha),Nw-1,nfft,'centered'); surf(t,w/pi,abs(s),'EdgeColor','none') view(2) axistightxlabel('Samples') ylabel('Normalized Frequency (\times\pi rad/sample)')

Figure contains an axes object. The axes object contains an object of type surface.

The Fourier synchrosqueezed transform results in a sharper, better localized estimate of the spectrum.

[ss,sw,st] = fsst(x,[],gausswin(Nw,alpha)); fsst(x,'yaxis')

Figure contains an axes object. The axes object with title Fourier Synchrosqueezed Transform contains an object of type image.

The sinusoids are visible as constant oscillations at the expected frequency values. To see that the decay away from the ridges is Gaussian, plot an instantaneous value of the transform and overlay two instances of a Gaussian. Express the Gaussian amplitude and standard deviation in terms of $α$ and the window length. Recall that the standard deviation of the frequency-domain window is the reciprocal of the standard deviation of the time-domain window.

rstdev = (Nw-1)/(2*alpha); amp = rstdev*sqrt(2*pi); instransf = abs(s(:,128)); plot(w/pi,instransf) holdonplot(w/pi,[1 3]*amp.*exp(-rstdev^2/2*(w-[1 3]*w0).^2),'--') holdoffxlabel('Normalized Frequency (\times\pi rad/sample)') lg = legend('Transform','First sinusoid','Second sinusoid'); lg.Location ='best';

Figure contains an axes object. The axes object contains 3 objects of type line. These objects represent Transform, First sinusoid, Second sinusoid.

傅里叶synchrosqueezed变换集中the energy content of the signal at the estimated instantaneous frequencies.

stem(sw/pi,abs(ss(:,128))) xlabel('Normalized Frequency (\times\pi rad/sample)') title('Synchrosqueezed Transform')

Figure contains an axes object. The axes object with title Synchrosqueezed Transform contains an object of type stem.

The synchrosqueezed estimates of instantaneous frequency are valid only if the sinusoids are separated by more than $2 Δ$ , where

$Δ = \frac{1}{σ} \sqrt{2 \log 2}$

for a Gaussian window and $σ$ is the standard deviation.

Repeat the previous calculation, but now specify that the second sinusoid has a normalized frequency of $ω_{0} + 1.9 Δ$ rad/sample.

D = sqrt(2*log(2))/rstdev; w1 = w0+1.9*D; x = exp(1j*w0*n)+3*exp(1j*w1*n); [s,w,t] = spectrogram(x,gausswin(Nw,alpha),Nw-1,nfft,'centered'); instransf = abs(s(:,20)); plot(w/pi,instransf) holdonplot(w/pi,[1 3]*amp.*exp(-rstdev^2/2*(w-[w0 w1]).^2),'--') holdoffxlabel('Normalized Frequency (\times\pi rad/sample)') lg = legend('Transform','First sinusoid','Second sinusoid'); lg.Location ='best';

Figure contains an axes object. The axes object contains 3 objects of type line. These objects represent Transform, First sinusoid, Second sinusoid.

The Fourier synchrosqueezed transform cannot resolve the sinusoids well because $| ω_{1} - ω_{0} | < 2 Δ$ . The spectral estimates decrease significantly in value.

[ss,sw,st] = fsst(x,[],gausswin(Nw,alpha)); stem(sw/pi,abs(ss(:,128))) xlabel('Normalized Frequency (\times\pi rad/sample)') title('Synchrosqueezed Transform')

Figure contains an axes object. The axes object with title Synchrosqueezed Transform contains an object of type stem.

Detect Closely Spaced Sinusoids

See Also

Related Topics