Практическое введение в частотно-временной анализ

В этом примере показано, как выполнить и интерпретировать основной анализ сигнала частоты времени. В практических применениях много сигналов являются неустановившимися. Это означает, что их представление частотного диапазона (их спектр) изменяется в зависимости от времени. Пример обсуждает преимущества использования методов частоты времени по частотному диапазону или представлениям временного интервала сигнала. Это отвечает на основные вопросы, такие как: Когда конкретная частотная составляющая существует в моем сигнале? Как я увеличиваю разрешение частоты или время? Как я могу увеличить резкость спектра компонента или извлечь конкретный режим? Как я измеряю степень в представлении частоты времени? Как я визуализирую информацию о частоте времени своего сигнала? Как я нахожу неустойчивую интерференцию в содержимом частоты сигнала интереса?

Используя частотно-временной анализ, чтобы идентифицировать числа в сигнале DTMF

Можно разделить почти любой изменяющийся во времени сигнал на временные интервалы, достаточно короткие, что сигнал является чрезвычайно стационарным в каждом разделе. Частотно-временной анализ обычно выполняется путем сегментации сигнала в те короткие периоды и оценки спектра по раздвижным окнам. pspectrum функция используется с 'spectrogram' опция вычисляет основанную на БПФ спектральную оценку по каждому раздвижному окну и позволяет вам визуализировать, как содержимое частоты сигнала изменяется в зависимости от времени.

Рассмотрите сигнальную систему цифрового телефонного набора. Сигналы, произведенные такой системой, известны как сигналы двухтонального многочастотного (DTMF). Звук, сгенерированный каждым набранным номером, состоит из суммы двух синусоид - или тоны - с частотами, взятыми из двух взаимоисключающих групп. Каждая пара тонов содержит одну частоту низкой группы (697 Гц, 770 Гц, 852 Гц или 941 Гц) и одну частоту высокой группы (1 209 Гц, 1 336 Гц или 1477 Гц) и представляет уникальный символ. Следующее является частотами, выделенными кнопкам телефонной клавиатуры:

Сгенерируйте DTMF, сигнализируют и слушают его.

[tones, Fs] = helperDTMFToneGenerator();
p = audioplayer(tones,Fs,16);
play(p)

Слушая сигнал, можно сказать, что трехзначное число было набрано. Однако вы не можете сказать, каким номером это было. Затем визуализируйте сигнал вовремя и в частотном диапазоне по полосе на 650 - 1 500 Гц. Установите 'Leakage' параметр pspectrum функционируйте к 1, чтобы использовать прямоугольное окно и улучшить разрешение частоты.

N = numel(tones);
t = (0:N-1)/Fs; 
subplot(2,1,1)
plot(1e3*t,tones)
xlabel('Time (ms)')
ylabel('Amplitude')
title('DTMF Signal')
subplot(2,1,2)
pspectrum(tones,Fs,'Leakage',1,'FrequencyLimits',[650, 1500])

Figure contains 2 axes. Axes 1 with title DTMF Signal contains an object of type line. Axes 2 with title Fres = 7.8144 Hz contains an object of type line.

График временного интервала сигнала подтверждает присутствие трех пакетов энергии, соответствуя трем нажатым кнопкам. Чтобы измерить длину пакета, можно взять ширину импульса конверта RMS.

env = envelope(tones,80,'rms');
pulsewidth(env,Fs)
ans = 3×1

    0.1041
    0.1042
    0.1047

title('Pulse Width of RMS Envelope')

Figure Pulse Width Plot contains an axes. The axes with title Pulse Width of RMS Envelope contains 10 objects of type patch, line. These objects represent pulse width, signal, mid cross, upper boundary, upper state, lower boundary, mid reference, lower state.

Здесь вы видите три импульса, каждый приблизительно 100 миллисекунд долго. Однако вы не можете сказать, какие номера были набраны. График частотного диапазона помогает вам понять это, потому что он показывает частоты, существующие в сигнале.

Найдите peaks частоты путем оценки средней частоты в четырех различных диапазонах частот.

f = [meanfreq(tones,Fs,[700 800]), ...
     meanfreq(tones,Fs,[800 900]), ...
     meanfreq(tones,Fs,[900 1000]), ...
     meanfreq(tones,Fs,[1300 1400])];
round(f)
ans = 1×4

         770         852         941        1336

Путем соответствия с предполагаемыми частотами к схеме телефонной клавиатуры можно сказать, что набранные кнопки равнялись '5', '8', и '0'. Однако график частотного диапазона не обеспечивает типа информации времени, которая позволила бы вам выяснять порядок, в котором они были набраны. Комбинация могла быть '580', '508', '805', '850', '085', или '058'. Чтобы решить эту проблему, используйте функцию pspectrum, чтобы вычислить спектрограмму и наблюдать, как содержимое частоты сигнала меняется в зависимости от времени.

Вычислите спектрограмму по полосе на 650 - 1 500 Гц и удалите содержимое ниже -Уровень мощности на 10 дБ, чтобы визуализировать только основные частотные составляющие. Чтобы видеть тональную длительность и их местоположения во время используют 0%-е перекрытие.

pspectrum(tones,Fs,'spectrogram','Leakage',1,'OverlapPercent',0, ...
    'MinThreshold',-10,'FrequencyLimits',[650, 1500]);

Figure Pulse Width Plot contains an axes. The axes with title Fres = 45.7145 Hz, Tres = 21.875 ms contains an object of type image.

Цвета спектрограммы кодируют уровни мощности частоты. Желтые цвета указывают на содержимое частоты с более высокой степенью; синие цвета указывают на содержимое частоты с очень малой мощностью. Сильная желтая горизонтальная линия указывает на существование тона на особой частоте. График ясно показывает присутствие тона на 1 336 Гц во всех трех набранных цифрах, говоря вам, что они - все на втором столбце клавиатуры. Из графика вы видите, что самая низкая частота, 770 Гц, была набрана сначала. Самая высокая частота, 941 Гц, была следующей. Средняя частота, 852 Гц, была последней. Следовательно, набранный номер был 508.

Обменивание разрешения времени и частоты, чтобы получить лучшее представление сигнала

pspectrum функция делит сигнал на сегменты. Более длинные сегменты обеспечивают лучшее разрешение частоты; более короткие сегменты обеспечивают лучшее разрешение времени. Длинами сегмента можно управлять с помощью 'FrequencyResolution' и 'TimeResolution' параметры. Когда никакое разрешение частоты или значения разрешения времени не заданы, pspectrum попытки найти хороший баланс между временем и разрешениями частоты на основе длины входного сигнала.

Рассмотрите следующий сигнал, произведенный на уровне 4 кГц, который состоит из фрагмента трели Тихоокеанской песни голубого кита:

load whaleTrill
p = audioplayer(whaleTrill,Fs,16);
play(p)

Сигнал трели состоит из обучения тональных импульсов. Посмотрите в то время сигнал и спектрограмма, полученная pspectrum когда никакое разрешение не задано и когда разрешение времени установлено в 10 миллисекунд. Установите 'Leakage' параметр к 1, чтобы использовать прямоугольные окна. Поскольку мы хотим локализовать положение времени импульсов, установите процент перекрытия на 0. Наконец, используйте 'MinThreshold' из -60 дБ, чтобы удалить фоновый шум из представления спектрограммы.

t = (0:length(whaleTrill)-1)/Fs;
figure
ax1 = subplot(3,1,1);
plot(t,whaleTrill)
ax2 = subplot(3,1,2);
pspectrum(whaleTrill,Fs,'spectrogram','OverlapPercent',0, ...
    'Leakage',1,'MinThreshold',-60)
colorbar(ax2,'off')
ax3 = subplot(3,1,3);
pspectrum(whaleTrill,Fs,'spectrogram','OverlapPercent',0, ...
    'Leakage',1,'MinThreshold',-60,'TimeResolution', 10e-3)
colorbar(ax3,'off')
linkaxes([ax1,ax2,ax3],'x')

Figure contains 3 axes. Axes 1 with title Fres = 21.2767 Hz, Tres = 47 ms contains an object of type line. Axes 2 contains an object of type image. Axes 3 with title Fres = 100.0004 Hz, Tres = 10 ms contains an object of type image.

Эти 47 разрешений времени миллисекунд выбраны pspectrum не мал достаточно, чтобы локализовать все импульсы трели в спектрограмме. С другой стороны, разрешения времени 10 миллисекунд достаточно, чтобы локализовать каждый импульс трели вовремя. Это становится еще более ясным, если мы масштабируем в несколько импульсов:

xlim([0.3 0.55])

Figure contains 3 axes. Axes 1 with title Fres = 21.2767 Hz, Tres = 47 ms contains an object of type line. Axes 2 contains an object of type image. Axes 3 with title Fres = 100.0004 Hz, Tres = 10 ms contains an object of type image.

Теперь загрузите сигнал, который состоит из импульса эхолотирования, испускаемого большой коричневой битой (Eptesicus fuscus). Сигнал измеряется с интервалом выборки 7 микросекунд. Анализируйте спектрограмму сигнала.

load batsignal
Fs = 1/DT;
figure
pspectrum(batsignal,Fs,'spectrogram')

Figure contains an axes. The axes with title Fres = 7.3336 kHz, Tres = 350 μs contains an object of type image.

Спектрограмма со значениями параметров по умолчанию показывает четыре крупных частотно-временных гребня. Уменьшайте значение разрешения частоты до 3 кГц, чтобы получить больше деталей об изменении частоты каждого гребня.

pspectrum(batsignal,Fs,'spectrogram','FrequencyResolution',3e3)

Figure contains an axes. The axes with title Fres = 3.0075 kHz, Tres = 854 μs contains an object of type image.

Заметьте, что теперь гребни частоты лучше локализуются в частоте. Однако, поскольку частота и разрешение времени обратно пропорциональны, разрешение времени спектрограммы значительно меньше. Установите перекрытие 99% сглаживать окна времени. Используйте 'MinThreshold' из -60 дБ, чтобы удалить нежелательное фоновое содержимое.

pspectrum(batsignal,Fs,'spectrogram','FrequencyResolution',3e3, ...
    'OverlapPercent',99,'MinTHreshold',-60)

Figure contains an axes. The axes with title Fres = 3.0075 kHz, Tres = 854 μs contains an object of type image.

Новые настройки дают к спектрограмме, которая ясно показывает четыре гребня частоты сигнала эхолотирования.

Переназначение частоты времени

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

pspectrum функция способна к оценке центра энергии для каждой спектральной оценки и во время и в частоту. Если вы повторно присваиваете энергию каждой оценки к интервалу, самому близкому к новому времени и центрам частоты, можно откорректировать для части утечки окна. Можно сделать это при помощи 'Reassign' параметр. Установка этого параметра на true вычисляет переприсвоенную спектрограмму сигнала.

pspectrum(batsignal,Fs,'spectrogram','FrequencyResolution',3e3, ...
    'OverlapPercent',99,'MinTHreshold',-60,'Reassign',true)

Figure contains an axes. The axes with title Fres = 3.0075 kHz, Tres = 854 μs contains an object of type image.

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

Восстановление частотно-временного гребня

Рассмотрите следующую запись, состоя из сигнала щебета, частота которого уменьшается в зависимости от времени и итоговый звук нащельной рейки.

load splat
p = audioplayer(y,Fs,16);
play(p)
pspectrum(y,Fs,'spectrogram')

Figure contains an axes. The axes with title Fres = 133.9285 Hz, Tres = 19.165 ms contains an object of type image.

Давайте восстановим фрагмент звука "нащельной рейки" путем извлечения гребня в плоскости частоты времени. Мы используем fsst увеличить резкость спектра шумной версии сигнал нащельной рейки, tfridge идентифицировать гребень звука щебета и ifsst восстановить щебет. Процесс denoises восстановленный сигнал.

Добавьте Гауссов шум во фрагмент щебета звука "нащельной рейки". Добавленный шум симулирует аудиозапись, взятую с недорогим микрофоном. Исследуйте частоту времени спектральное содержимое.

rng('default')
t = (0:length(y)-1)/Fs;
yNoise = y + 0.1*randn(size(y));
yChirp = yNoise(t<0.35);
pspectrum(yChirp,Fs,'spectrogram','MinThreshold',-70)

Figure contains an axes. The axes with title Fres = 116.8154 Hz, Tres = 21.9727 ms contains an object of type image.

Увеличьте резкость спектра с помощью synchrosqueezed преобразования Фурье, fsst. fsst локализует энергию в плоскости частоты времени путем переприсвоения энергии в частоте в течение установленного времени. Вычислите и постройте synchrosqueezed преобразование шумного щебета.

fsst(yChirp,Fs,'yaxis')

Figure contains an axes. The axes with title Fourier Synchrosqueezed Transform contains an object of type image.

Щебет появляется как локализованный гребень в плоскости частоты времени. Идентифицируйте гребень с помощью tfridge. Постройте гребень наряду с преобразованием.

[sst,f] = fsst(yChirp,Fs); 
[fridge, iridge] = tfridge(sst,f,10);
helperPlotRidge(yChirp,Fs,fridge);

Figure contains an axes. The axes with title Time-Frequency Ridge of a Chirp contains 2 objects of type image, line.

Затем восстановите сигнал щебета использование гребенчатого вектора индекса iridge. Включайте один интервал на каждой стороне гребня. Постройте спектрограмму восстановленного сигнала.

yrec = ifsst(sst,kaiser(256,10),iridge,'NumFrequencyBins',1);
pspectrum(yrec,Fs,'spectrogram','MinThreshold',-70)

Figure contains an axes. The axes with title Fres = 116.8154 Hz, Tres = 21.9727 ms contains an object of type image.

Восстановление гребня удалило шум из сигнала. Проигрывайте шумное и сигналы denoised последовательно, чтобы услышать различие.

p = audioplayer([yChirp;zeros(size(yChirp));yrec],Fs,16);
play(p);

Измерение степени

Рассмотрите комплексный импульс линейной частоты модулируется (LFM), который является общей радарной формой волны. Вычислите спектрограмму сигнала с помощью разрешения времени 1,27 микросекунд и 90%-го перекрытия.

Fs = 1e8;
bw = 60e6;
t = 0:1/Fs:10e-6;
IComp = chirp(t,-bw/2,t(end), bw/2,'linear',90)+0.15*randn(size(t));
QComp = chirp(t,-bw/2,t(end), bw/2,'linear',0) +0.15*randn(size(t));
IQData = IComp + 1i*QComp;

segmentLength = 128;
pspectrum(IQData,Fs,'spectrogram','TimeResolution',1.27e-6,'OverlapPercent',90)

Figure contains an axes. The axes with title Fres = 2.0371 MHz, Tres = 1.26 μs contains an object of type image.

Параметры использовались для расчета спектрограммы, дают ясное представление частоты времени сигнала LFM. pspectrum вычисляет спектрограмму степени, это означает, что значения цвета соответствуют истинным уровням мощности в дБ. Цветная полоса показывает, что уровень мощности сигнала вокруг -4 дБ.

Логарифмическая визуализация шкалы частоты

В определенных приложениях может быть предпочтительно визуализировать спектрограмму сигнала по логарифмической шкале частоты. Можно достигнуть этого путем изменения YScale свойство оси Y. Например, считайте логарифмический щебет произведенным на уровне 1 кГц. Частота щебета увеличивается с 10 Гц до 400 Гц за 10 секунд.

Fs = 1e3;                    
t = 0:1/Fs:10;               
fo = 10;                     
f1 = 400;                    
y = chirp(t,fo,10,f1,'logarithmic');
pspectrum(y,Fs,'spectrogram','FrequencyResolution',1, ...
    'OverlapPercent',90,'Leakage',0.85,'FrequencyLimits',[1 Fs/2])

Figure contains an axes. The axes with title Fres = 1.0002 Hz, Tres = 1.468 s contains an object of type image.

Спектрограмма щебета становится прямой линией, когда шкала частоты является логарифмической.

ax = gca;
ax.YScale = 'log';

Figure contains an axes. The axes with title Fres = 1.0002 Hz, Tres = 1.468 s contains an object of type surface.

3D визуализация водопада

С view команда, можно визуализировать спектрограмму сигнала как 3D график водопада. Можно также изменить цвета отображения с colormap функция.

Fs = 10e3;
t = 0:1/Fs:2;
x1 = vco(sawtooth(2*pi*t,0.5),[0.1 0.4]*Fs,Fs);
pspectrum(x1,Fs,'spectrogram','Leakage',0.8)

Figure contains an axes. The axes with title Fres = 106.1965 Hz, Tres = 15.7 ms contains an object of type image.

view(-45,65)
colormap bone

Figure contains an axes. The axes with title Fres = 106.1965 Hz, Tres = 15.7 ms contains an object of type surface.

Нахождение интерференции Используя спектр персистентности

Спектр персистентности сигнала является представлением частоты времени, которое показывает процент времени, когда данная частота присутствует в сигнале. Спектр персистентности является гистограммой в пространстве частоты степени. Чем дольше особая частота сохраняется в сигнале, когда сигнал развивается, тем выше его процент времени и таким образом более яркое или "более горячее" его цвет в отображении. Используйте спектр персистентности, чтобы идентифицировать сигналы, скрытые в других сигналах.

Считайте интерференционный узкополосный сигнал встроенным в широкополосном сигнале. Сгенерируйте щебет, произведенный на уровне 1 кГц в течение 500 секунд. Частота щебета увеличивается с 180 Гц до 220 Гц во время измерения.

fs = 1000;
t = (0:1/fs:500)';

x = chirp(t,180,t(end),220) + 0.15*randn(size(t));

Сигнал также содержит интерференцию на 210 Гц с амплитудой 0,05, который присутствует только для 1/6 общей длительности сигнала.

idx = floor(length(x)/6);
x(1:idx) = x(1:idx) + 0.05*cos(2*pi*t(1:idx)*210);

Вычислите спектр мощности сигнала на интервале на 100 - 290 Гц. Слабая синусоида затенена щебетом.

pspectrum(x,fs,'FrequencyLimits',[100 290])

Figure contains an axes. The axes with title Fres = 976.801 mHz contains an object of type line.

Вычислите спектр персистентности сигнала. Теперь оба компонента сигнала ясно отображаются.

figure
colormap parula
pspectrum(x,fs,'persistence','FrequencyLimits',[100 290],'TimeResolution',1)

Figure contains an axes. The axes with title Fres = 3.9101 Hz, Tres = 1 s contains an object of type image.

Заключения

В этом примере вы изучили, как выполнить частотно-временной анализ с помощью функции pspectrum и как интерпретировать данные о спектрограмме и уровни мощности. Вы изучили, как изменить время и разрешение частоты, чтобы улучшить ваше понимание сигнала и как увеличить резкость спектров и частотно-временных гребней извлечения с помощью fsst, ifsst, и tfridge. Вы изучили, как сконфигурировать график спектрограммы получить логарифмическую шкалу частоты и 3D визуализацию. Наконец, вы изучили, как найти интерференционные сигналы путем вычисления спектра персистентности.

Приложение

Следующие функции помощника используются в этом примере.

Смотрите также

| | |

Для просмотра документации необходимо авторизоваться на сайте