Использование частотно-временного анализа для идентификации чисел в DTMF-сигнале

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

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

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

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

Слушая сигнал, можно сказать, что был набран трехзначный номер. Однако вы не можете сказать, какой это был номер. Затем визуализируйте сигнал во времени и в частотный диапазон с частотой от 650 до 1500 Гц полосой. Установите '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 до 1500 Гц и удалите содержимое ниже $-$10 дБ степени, чтобы визуализировать только основные частотные составляющие. Чтобы увидеть длительности тональных сигналов и их местоположения во времени, используйте перекрытие 0%.

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

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

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

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

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

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

Сигнал trill состоит из 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)

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

Частотно-временное переназначение

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

The 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 чтобы восстановить щебет. Процесс денонсирует восстановленный сигнал.

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

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 локализует энергию во временной частотной плоскости путем переназначения энергии по частоте на фиксированное время. Вычислите и постройте график синхронизированного преобразования шумного щебета.

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)

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

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';

Трехмерная визуализация водопада

С view команда позволяет визуализировать спектрограмму сигнала как трехмерный график водопада. Вы также можете изменить цвета отображения с помощью 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. Вы научились конфигурировать спектрограммный график, чтобы получить логарифмическую частотную шкалу и трехмерную визуализацию. Наконец, вы научились находить сигналы интерференции путем вычисления спектра стойкости.

Приложение

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