Обнаружение близко расположенных синусоидов

Открыть Live Script

Рассмотрим синусоиду, $f (x) = e^{j 2 π ν x}$ , оконный Гауссовым окном, $g (t) = e^{- π t^{2}}$ . Кратковременное преобразование

$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} .$

При рассмотрении как функции частоты преобразование комбинирует постоянное (во времени) колебание на $ν$ с Гауссовым распадом от него. Оценка синхронизации мгновенной частоты,

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

равен значению, полученному при помощи стандартного определения, $(2 π)^{- 1} d \arg f (x) / d x$ . Для наложения синусоидов,

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

краткосрочное преобразование становится

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

Создать 1024 выборки сигнала, состоящего из двух синусоидов. Одна синусоида имеет нормированную частоту $ω_{0} = π / 5$ рад/образец. Другая синусоида имеет втрое большую частоту и втрое большую амплитуду.

N = 1024;
n = 0:N-1;

w0 = pi/5;
x = exp(1j*w0*n)+3*exp(1j*3*w0*n);

Вычислите кратковременное преобразование Фурье сигнала. Используйте 256-образное Гауссово окно с $α = 20$ , 255 выборки перекрытия между смежными секциями и 1024 точки ДПФ. Постройте график абсолютного значения преобразования.

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)
axis tight
xlabel('Samples')
ylabel('Normalized Frequency (\times\pi rad/sample)')

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

Synchrosqueezed преобразование Фурье результатов в более резкой, лучше локализованной оценке спектра.

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

fsst(x,'yaxis')

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

Синусоиды видны как постоянные колебания при ожидаемых значениях частоты. Чтобы увидеть, что распад от гребней Гауссов, постройте график мгновенного значения преобразования и наложите два образцов Гауссова. Выразите Гауссову амплитуду и стандартное отклонение в терминах $α$ и длина окна. Напомним, что стандартное отклонение окна частотного диапазона является обратным стандартному отклонению окна временной области.

rstdev = (Nw-1)/(2*alpha);
amp = rstdev*sqrt(2*pi);

instransf = abs(s(:,128));

plot(w/pi,instransf)
hold on
plot(w/pi,[1 3]*amp.*exp(-rstdev^2/2*(w-[1 3]*w0).^2),'--')
hold off
xlabel('Normalized Frequency (\times\pi rad/sample)')
lg = legend('Transform','First sinusoid','Second sinusoid');
lg.Location = 'best';

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

Synchrosqueezed преобразование Фурье концентрирует энергетическое содержимое сигнала на оцененных мгновенных частотах.

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

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

Синхронизированные оценки мгновенной частоты действительны только, если синусоиды разделены более чем $2 Δ$ , где

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

для Гауссова окна и $σ$ - стандартное отклонение.

Повторите предыдущее вычисление, но теперь задайте, что вторая синусоида имеет нормированную частоту $ω_{0} + 1.9 Δ$ рад/образец.

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)
hold on
plot(w/pi,[1 3]*amp.*exp(-rstdev^2/2*(w-[w0 w1]).^2),'--')
hold off
xlabel('Normalized Frequency (\times\pi rad/sample)')
lg = legend('Transform','First sinusoid','Second sinusoid');
lg.Location = 'best';

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

Synchrosqueezed преобразование Фурье не может хорошо разрешить синусоиды, потому что $| ω_{1} - ω_{0} | < 2 Δ$ . Спектральные оценки значительно уменьшаются в значении.

[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. The axes with title Synchrosqueezed Transform contains an object of type stem.

См. также

fsst | gausswin | ifsst | spectrogram

Документация

Обнаружение близко расположенных синусоидов

См. также

Похожие темы

Signal Processing Toolbox документация

Поддержка