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

Скрипт Open Live Script

Рассмотрите синусоиду, $f (x) = e^{j 2 π ν x}$ , оконный с окном Gaussian, $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} .$

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

$Ω_{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} .$

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

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

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

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

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)')

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

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

fsst(x,'yaxis')

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

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';

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

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

synchrosqueezed оценки мгновенной частоты допустимы, только если синусоиды разделяются больше, чем $2 Δ$ , где

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

для окна Gaussian и $σ$ стандартное отклонение.

Повторите предыдущее вычисление, но теперь укажите, что вторая синусоида имеет нормированную частоту $ω_{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';

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')

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

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

Смотрите также

Похожие темы

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

Поддержка