Сгенерируйте 1 024 выборки сигнала, который состоит из суммы синусоид, встроенных в белый Гауссов шум. Нормированные частоты синусоид $$2\pi /5$$ рад/выборка и $$4\pi /5$$ рад/выборка. Более высокая синусоида частоты имеет 3 раза амплитуду другой синусоиды.

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

w0 = 2*pi/5;
x = sin(w0*n)+3*sin(2*w0*n);

Вычислите synchrosqueezed преобразование Фурье сигнала. Постройте результат.

[s,w,n] = fsst(x);

mesh(n,w/pi,abs(s))

axis tight
view(2)
colorbar

Вычислите обычное кратковременное преобразование Фурье сигнала для сравнения. Используйте значения по умолчанию spectrogram. Постройте результат.

[s,w,n] = spectrogram(x);
 
surf(n,w/pi,abs(s),'EdgeColor','none')

axis tight
view(2)
colorbar

Сгенерируйте сигнал, который состоит из двух щебетов. Сигнал выбирается на уровне 3 кГц в течение одной секунды. Первый щебет имеет начальную частоту 400 Гц и достигает 800 Гц в конце выборки. Второй щебет запускается на уровне 500 Гц и достигает 1 000 Гц в конце. Второй щебет имеет дважды амплитуду первого щебета.

fs = 3000;
t = 0:1/fs:1-1/fs;

x1 = chirp(t,400,t(end),800);
x2 = 2*chirp(t,500,t(end),1000);

Вычислите и постройте synchrosqueezed преобразование Фурье сигнала.

fsst(x1+x2,fs,'yaxis')

Сравните synchrosqueezed преобразовывают с кратковременным преобразованием Фурье (STFT). Вычислите STFT использование функции spectrogram. Задайте параметры по умолчанию, используемые fsst:

[stft,f,t] = spectrogram(x1+x2,kaiser(256,10),255,256,fs);

Постройте абсолютное значение STFT.

mesh(t,f,abs(stft))

xlabel('Time (s)') 
ylabel('Frequency (Hz)')
title('Short-Time Fourier Transform')
axis tight
view(2)

Вычислите и отобразите synchrosqueezed преобразование Фурье квадратичного щебета, который запускается на уровне 100 Гц и пересекает 200 Гц в t = 1 с. Задайте частоту дискретизации 1 кГц. Выразите шаг расчета как скаляр duration.

fs = 1000;
t = 0:1/fs:2;
ts = duration(0,0,1/fs);

x = chirp(t,100,1,200,'quadratic');

fsst(x,ts,'yaxis')

title('Quadratic Chirp')

synchrosqueezing алгоритм работает под предположением, что частота сигнала медленно отличается. Таким образом спектр лучше сконцентрирован в прежние времена, где скорость изменения частоты меньше.

Вычислите и отобразите synchrosqueezed преобразование Фурье линейного щебета, который запускается в DC и пересекает 150 Гц в t = 1 с. Используйте Окно Хэмминга с 256 выборками.

x = chirp(t,0,1,150);

fsst(x,ts,hamming(256),'yaxis')

title('Linear Chirp')

Вычислите и отобразите synchrosqueezed преобразование Фурье логарифмического щебета. Щебет выбирается на уровне 1 кГц, запускается на уровне 20 Гц и пересекает 60 Гц в t = 1 с. Используйте окно Kaiser с 256 выборками с β = 20.

x = chirp(t,20,1,60,'logarithmic');

[s,f,t] = fsst(x,fs,kaiser(256,20));

clf
mesh(t,f,(abs(s)))

title('Logarithmic Chirp') 
xlabel('Time (s)')
ylabel('Frequency (Hz)')
view(2)

Используйте логарифмический масштаб для оси частоты. Преобразование становится прямой линией.

ax = gca;
ax.YScale = 'log';
axis tight

Загрузите речевой сигнал, выбранный в $$F_s=\mathrm{}\mathrm{}\mathrm{}7418Hz$$. Файл содержит запись розеточной речи, говоря слово "MATLAB®".

load mtlb

% To hear, type sound(mtlb,Fs)

Вычислите synchrosqueezed преобразование сигнала. Используйте окно Hann длины 256. Отобразите время на оси X и частоту на оси Y.

fsst(mtlb,Fs,hann(256),'yaxis')

Используйте ifsst, чтобы инвертировать преобразование. Сравните исходные и восстановленные сигналы.

sst = fsst(mtlb,Fs,hann(256));

xrc = ifsst(sst,hann(256));

plot((0:length(mtlb)-1)/Fs,[mtlb xrc xrc-mtlb])
legend('Original','Reconstructed','Difference')

% To hear, type sound(xrc-mtlb,Fs)