Загрузка гиперболического щебета

Загрузите сигнал, который имеет два гиперболических щебета. Данные отбираются с частотой дискретизации 2048 Гц. Первый щебет активен между 0,1 и 0,68 секунд, а второй щебет активен между 0,1 и 0,75 секунд. Мгновенная частота (в hertz) первого щебета в момент времени $$t$$ является $$\frac{15\pi }{{(0.8-t)}^2}/2\pi$$ . Мгновенная частота второго щебета в момент времени $$t$$ является $$\frac{5\pi }{{(0.8-t)}^2}/2\pi$$. Постройте график сигнала.

load hychirp
plot(t,hychirp)
grid on
title('Signal')
axis tight
xlabel('Time (s)')
ylabel('Amplitude')

Частотно-временной анализ: преобразование Фурье

Преобразование Фурье (FT) очень хорошо идентифицирует частотные составляющие, присутствующие в сигнале. Однако FT не определяет, когда происходят частотные составляющие.

Постройте график спектра величины сигнала. Изменение масштаба области от 0 до 200 Гц.

sigLen = numel(hychirp);
fchirp = fft(hychirp);
fr = Fs*(0:1/Fs:1-1/Fs);
plot(fr(1:sigLen/2),abs(fchirp(1:sigLen/2)),'x-')
xlabel('Frequency (Hz)')
ylabel('Amplitude')
axis tight
grid on
xlim([0 200])

Частотно-временной анализ: короткое преобразование Фурье

Преобразование Фурье не предоставляет информацию о времени. Чтобы определить, когда происходят изменения частоты, кратковременный подход преобразования Фурье (STFT) сегментирует сигнал в различные фрагменты и выполняет FT на каждом фрагменте. Здесь показана мозаика STFT в частотно-временной плоскости.

STFT предоставляет некоторую информацию как о времени, так и частотах, на которых происходит событие сигнала. Однако выбор размера окна (сегмента) является ключевым. Для частотно-временного анализа с использованием STFT, выбор более короткого размера окна помогает получить хорошее разрешение времени за счет разрешения частоты. И наоборот, выбор большего окна помогает получить хорошее разрешение частоты за счет временного разрешения.

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

Мгновенные частоты двух щебетаний в их начальных точках времени составляют приблизительно 5 Гц и 15 Гц. Используйте функцию helper helperPlotSpectrogram для построения спектрограммы сигнала с размером временного окна 200 миллисекунд. Исходный код для helperPlotSpectrogram приведено в приложении. Графики функций мгновенных частот по спектрограмме как черные штриховые сегменты. Мгновенные частоты разрешаются раньше в сигнале, но не так хорошо позже.

helperPlotSpectrogram(hychirp,t,Fs,200)

Теперь используйте helperPlotSpectrogram для построения графика спектрограммы с размером временного окна 50 миллисекунд. Более высокие частоты, которые происходят позже в сигнале, теперь разрешены, но более низкие частоты в начале сигнала не являются.

helperPlotSpectrogram(hychirp,t,Fs,50)

Для нестационарных сигналов, таких как гиперболический щебет, использование STFT проблематично. Ни один размер окна не может разрешить все частотное содержимое таких сигналов.

Частотно-временной анализ: непрерывный Вейвлет преобразование

Непрерывное вейвлет (CWT) было создано, чтобы преодолеть проблемы разрешения, присущие STFT. Здесь показано плиточное размещение CWT на частотно-временной плоскости.

CWT- плиточного размещения плоскости полезна, потому что многие реальные сигналы имеют медленно колеблющиеся содержимое, которые происходят на длинных шкалах, в то время как высокие частоты события, как правило, являются резкими или переходными. Однако, если бы было естественно, чтобы высокочастотные события были длинными по длительности, то использование CWT было бы нецелесообразным. Вы бы имели более низкое разрешение частоты, не получая никакого временного разрешения. Но это довольно часто не так. Слуховая система человека работает таким образом; у нас есть гораздо лучшая локализация частот на более низких частотах и лучшая локализация времени на высоких частотах.

Постройте график скалограммы CWT. Скалограмма является абсолютным значением CWT, нанесенным как функция времени и частоты. График использует логарифмическую ось частоты, потому что частоты в CWT являются логарифмическими. Наличие двух гиперболических щебетаний в сигнале ясно из скалограммы. С помощью CWT можно точно оценить мгновенные частоты на протяжении всей длительности сигнала, не беспокоясь о выборе длины сегмента.

cwt(hychirp,Fs)

Белая штриховая линия помечает то, что известно как конус влияния. Конус влияния показывает области в скалограмме, потенциально затронутые краевыми эффектами. Для получения дополнительной информации см. «Граничные Эффекты и конус влияния».

Чтобы понять, как быстро растет величина коэффициентов вейвлета, используйте функцию helper helperPlotScalogram3d для построения скалограммы как 3-D поверхности. Исходный код для helperPlotScalogram3d приведено в приложении.

helperPlotScalogram3d(hychirp,Fs)

Используйте функцию helper helperPlotScalogram для построения скалограммы сигнала и мгновенных частот. Исходный код для helperPlotScalogram приведено в приложении. Мгновенные частоты хорошо совпадают с функциями скалограммы.

helperPlotScalogram(hychirp,Fs)

Приложение - Вспомогательные функции

helperPlotSpectrogram

function helperPlotSpectrogram(sig,t,Fs,timeres)
% This function is only intended to support this wavelet example.
% It may change or be removed in a future release.

[px,fx,tx] = pspectrum(sig,Fs,'spectrogram','TimeResolution',timeres/1000);
hp = pcolor(tx,fx,20*log10(abs(px)));
hp.EdgeAlpha = 0;
ylims = hp.Parent.YLim;
yticks = hp.Parent.YTick;
cl = colorbar;
cl.Label.String = 'Power (dB)';
axis tight
hold on
title(['Time Resolution: ',num2str(timeres),' ms'])
xlabel('Time (s)')
ylabel('Hz');
dt  = 1/Fs;
idxbegin = round(0.1/dt);
idxend1 = round(0.68/dt);
idxend2 = round(0.75/dt);
instfreq1 = abs((15*pi)./(0.8-t).^2)./(2*pi);
instfreq2 = abs((5*pi)./(0.8-t).^2)./(2*pi);
plot(t(idxbegin:idxend1),(instfreq1(idxbegin:idxend1)),'k--');
hold on;
plot(t(idxbegin:idxend2),(instfreq2(idxbegin:idxend2)),'k--');
ylim(ylims);
hp.Parent.YTick = yticks;
hp.Parent.YTickLabels = yticks;
hold off
end

helperPlotScalogram

function helperPlotScalogram(sig,Fs)
% This function is only intended to support this wavelet example.
% It may change or be removed in a future release.
[cfs,f] = cwt(sig,Fs);

sigLen = numel(sig);
t = (0:sigLen-1)/Fs;

hp = pcolor(t,log2(f),abs(cfs));
hp.EdgeAlpha = 0;
ylims = hp.Parent.YLim;
yticks = hp.Parent.YTick;
cl = colorbar;
cl.Label.String = 'magnitude';
axis tight
hold on
title('Scalogram and Instantaneous Frequencies')
xlabel('Seconds');
ylabel('Hz');
dt  = 1/2048;
idxbegin = round(0.1/dt);
idxend1 = round(0.68/dt);
idxend2 = round(0.75/dt);
instfreq1 = abs((15*pi)./(0.8-t).^2)./(2*pi);
instfreq2 = abs((5*pi)./(0.8-t).^2)./(2*pi);
plot(t(idxbegin:idxend1),log2(instfreq1(idxbegin:idxend1)),'k--');
hold on;
plot(t(idxbegin:idxend2),log2(instfreq2(idxbegin:idxend2)),'k--');
ylim(ylims);
hp.Parent.YTick = yticks;
hp.Parent.YTickLabels = 2.^yticks;
end

helperPlotScalogram3d

function helperPlotScalogram3d(sig,Fs)
% This function is only intended to support this wavelet example.
% It may change or be removed in a future release.
figure
[cfs,f] = cwt(sig,Fs);

sigLen = numel(sig);
t = (0:sigLen-1)/Fs;
surface(t,f,abs(cfs));
xlabel('Time (s)')
ylabel('Frequency (Hz)')
zlabel('Magnitude')
title('Scalogram In 3-D')
set(gca,'yscale','log')
shading interp
view([-40 30])
end