Документация

Используя анализ частоты времени, чтобы идентифицировать числа в сигнале 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])

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

env = envelope(tones,80,'rms');
pulsewidth(env,Fs)

ans = 3×1

    0.1041
    0.1042
    0.1047

title('Pulse Width of RMS Envelope')

Здесь вы видите три импульса, каждый приблизительно 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]);

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

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

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

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

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

Сигнал трели состоит из train тональных импульсов. Посмотрите в то время сигнал и спектрограмма, полученная 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')

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

xlim([0.3 0.55])

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Давайте восстановим фрагмент звука "нащельной рейки" путем извлечения гребня в плоскости частоты времени. Мы используем 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)

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

fsst(yChirp,Fs,'yaxis')

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

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

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

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

Восстановление гребня удалило шум из сигнала. Проигрывайте шумное и сигналы 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)

Параметры, используемые, чтобы вычислить спектрограмму, дают ясное представление частоты времени сигнала 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])

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

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

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)

view(-45,65)
colormap bone

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

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

Считайте интерференционный узкополосный сигнал встроенным в широкополосном сигнале. Сгенерируйте щебет, выбранный на уровне 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
colormap parula
pspectrum(x,fs,'persistence','FrequencyLimits',[100 290],'TimeResolution',1)

Заключения

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

Приложение

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