Создайте матрицу, которая напоминает частотно-временную матрицу с острым гребнем. Визуализируйте матрицу в трёх размерностях.

t = 0:0.05:10;
f = 0:0.2:8;
rv = 1;

[F,T] = ndgrid(f,t);

S = zeros(size(T));
S(abs((F-4)-cos((T-6).^2))<0.1) = rv;

mesh(T,F,S)
view(-30,60)

Добавьте шум в матрицу и повторно отобразите график.

S = S+rand(size(S))/10;

mesh(T,F,S)
view(-30,60)
xlabel('Time')
ylabel('Frequency')

Извлеките гребень и постройте график результата.

[fridge,~,lridge] = tfridge(S,f);

rvals = S(lridge);

hold on
plot3(t,fridge,rvals,'k','linewidth',4)
hold off

Сгенерируйте сигнал, который дискретизируется с частотой 3 кГц в течение одной секунды. Сигнал состоит из двух тонов и квадратичного щебета.

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

x1 = 6*chirp(t,fs/6,t(end),fs/4,'quadratic');
x2 = sin(2*pi*fs/3*t);
x3 = sin(2*pi*fs/2.5*t);

x = x1+x2+x3;

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

[sst,f] = fsst(x,fs);
mx = max(abs(sst(:)))*ones(size(t));

mesh(t,f,abs(sst))
view(2)

Извлеките и постройте график двух компонентов сигнала самой высокой энергии. Не устанавливайте штраф за изменение частоты.

penval = 0;

fridge = tfridge(sst,f,penval,'NumRidges',2);

hold on
plot3(t,fridge,mx,'w','linewidth',5)
hold off

Два тона имеют одинаковую амплитуду, и алгоритм скачет между ними. Установите штраф за изменение частоты равной 1.

penval = 1;

fridge = tfridge(sst,f,penval,'NumRidges',2);

mesh(t,f,abs(sst))
view(2)
xlabel('Time (s)')
ylabel('Frequency (Hz)')

hold on
plot3(t,fridge,mx,'w','linewidth',5)
hold off

Установите штраф в высокое значение для сравнения. Щебет наказывается, потому что его частота не постоянна.

penval = 1000;

fridge = tfridge(sst,f,penval,'NumRidges',2);

mesh(t,f,abs(sst))
view(2)
xlabel('Time (s)')
ylabel('Frequency (Hz)')

hold on
plot3(t,fridge,mx,'w','linewidth',5)
hold off

Сгенерируйте сигнал, состоящий из двух квадратичных щебетаний. Дискретизация сигнала производится на частоте 1 кГц в течение 3 секунд. Щебеты таковы, что мгновенная частота симметрична относительно полутоновой точки интервала дискретизации. Один щебет вогнутый, а другой щебет выпуклый. Вогнутый щебет имеет удвоенную амплитуду выпуклого щебета.

fs = 1e3;
t = 0:1/fs:3;

x =   chirp(t-1.5,100,1.1,200,'quadratic',[],'convex');
y = 2*chirp(t-1.5,300,1.1,400,'quadratic',[],'concave');

% To hear, type soundsc(x+y,fs)

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

sig = x+y;

[sst,f,t] = fsst(sig,fs);

fsst(sig,fs,'yaxis')

Извлеките два частотно-временных гребней, которые имеют самую высокую энергию. Укажите штраф в размере 1 за изменения частоты. Удалите 1 интервал частоты вокруг хребта с самой высокой энергией перед извлечением второго гребня. Постройте график гребней.

nfb = 1;
[fr,ir] = tfridge(sst,f,1,'NumRidges',2,'NumFrequencyBins',nfb);

plot(t,fr)

Одного интервал недостаточно: Функция находит второй хребет, который частично находится на склоне первого. Увеличьте до 50 интервалов, чтобы удалить и повторить расчет.

nfb = 50;
[fr,ir] = tfridge(sst,f,1,'NumRidges',2,'NumFrequencyBins',nfb);

plot(t,fr)

Удаление слишком большого количества интервалов искажает гребни с более низкой энергией. Уменьшите число до 15 и повторите расчет.

nfb = 15;
[fr,ir] = tfridge(sst,f,1,'NumRidges',2,'NumFrequencyBins',nfb);

plot(t,fr)

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

itr = ifsst(sst,[],ir,'NumFrequencyBins',nfb);
xrec = sum(itr');

plot(t,xrec-(x+y))
ylim([-.1 .1])

% To hear, type soundsc(xrec,fs)

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