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

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;

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

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

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

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)

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