Нагрузка гиперболического чирпа

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

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 Гц. Использовать функцию помощника 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)

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

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

helperPlotScalogram3d(hychirp,Fs)

Использовать функцию помощника 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