Преобразование Гильберта оценивает мгновенную частоту сигнала только для монокомпонентных сигналов. Монокомпонентный сигнал описывается в частотно-временной плоскости одним «гребнем». Набор монокомпонентных сигналов включает одинарные синусоиды и такие сигналы, как щебет.
Сгенерируйте щебет, дискретизированный с частотой дискретизации 1 кГц в течение двух секунд. Укажите щебет так, чтобы его частота первоначально составляла 100 Гц и увеличивалась до 200 Гц через одну секунду.
fs = 1000; t = 0:1/fs:2-1/fs; y = chirp(t,100,1,200);
Оцените спектрограмму щебета с помощью краткосрочного преобразования Фурье, реализованного в pspectrum функция. Сигнал хорошо описан одной пиковой частотой в каждую точку времени.
pspectrum(y,fs,'spectrogram')
Вычислите аналитический сигнал и дифференцируйте его фазу, чтобы измерить мгновенную частоту. Масштабированная производная дает значимую оценку.
z = hilbert(y); instfrq = fs/(2*pi)*diff(unwrap(angle(z))); clf plot(t(2:end),instfrq) ylim([0 fs/2])
The instfreq функция вычисляет и отображает текущую частоту за один шаг.
instfreq(y,fs,'Method','hilbert')
Метод прекращает работать, когда сигнал не является монокомпонентным.
Сгенерируйте сумму двух синусоидов частот 60 Гц и 90 Гц, дискретизированных при частоте 1023 Гц в течение двух секунд. Вычислите и постройте график спектрограммы. Каждая временная точка показывает присутствие двух компонентов.
fs = 1023;
t = 0:1/fs:2-1/fs;
x = sin(2*pi*60*t)+sin(2*pi*90*t);
pspectrum(x,fs,'spectrogram')
yticks([60 90])
Вычислите аналитический сигнал и дифференцируйте его фазу. Увеличьте изображение области, охватывающей частоты синусоидов. Аналитический сигнал предсказывает мгновенную частоту, которая является средним значением синусоидальных частот.
z = hilbert(x); instfrq = fs/(2*pi)*diff(unwrap(angle(z))); plot(t(2:end),instfrq) ylim([60 90]) xlabel('Time (s)') ylabel('Frequency (Hz)')
The instfreq функция также оценивает среднее значение.
instfreq(x,fs,'Method','hilbert')
Чтобы оценить обе частоты как функции времени, используйте spectrogram для нахождения спектральной плотности и tfridge степени для отслеживания двух гребней. В tfridge, укажите штраф за изменение частоты 0,1.
[s,f,tt] = pspectrum(x,fs,'spectrogram'); numcomp = 2; [fridge,~,lr] = tfridge(s,f,0.1,'NumRidges',numcomp); pspectrum(x,fs,'spectrogram') hold on plot3(tt,fridge,abs(s(lr)),'LineWidth',4) hold off yticks([60 90])
