Сгенерируйте 1024 выборки сигнала, который состоит из суммы синусоидов, встроенных в белый Гауссов шум. Нормированные частоты синусоидов $$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 Гц и достигает 1000 Гц в конце. Второй щебет имеет удвоенную амплитуду первого щебета.

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

Сравните синхронизированное преобразование с коротким преобразованием Фурье (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')

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

Вычислите и отобразите synchrosqueezed преобразование Фурье линейного щебета, который начинается с постоянного тока и пересекает 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 с. Используйте 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=7418Hz$$. Файл содержит запись женского голоса, говорящего слово «MATLAB ®».

load mtlb

% To hear, type sound(mtlb,Fs)

Вычислите синхронизированное преобразование сигнала. Используйте окно Ханна длиной 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)