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

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)

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