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

Можно разделить почти любой изменяющийся во времени сигнал на интервалы времени, достаточно короткие, чтобы сигнал был по существу неподвижным в каждой секции. Частотно-временной анализ чаще всего выполняется путем сегментирования сигнала на эти короткие периоды и оценки спектра по скользящим окнам. 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 миллисекунд. Однако вы не можете сказать, какие номера были набраны. График частотной области помогает понять это, потому что он показывает частоты, присутствующие в сигнале.

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

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.

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

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

Разрешение 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 для точения спектра шумного варианта сигнала splat, tfridge для идентификации гребня чирп-звука, и ifsst реконструировать чирп. Процесс отключает восстановленный сигнал.

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

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)

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

Приложение

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