Частотно-временной анализ модулированных сигналов

Загрузите квадратичный щебет и отобразите график его спектрограммы. Частота сигнала начинается приблизительно с 500 Гц при t = 0, уменьшается до 100 Гц при t = 2 и увеличивается до 500 Гц при t = 4. Частота дискретизации составляет 1 кГц.

load quadchirp;
fs = 1000;
[S,F,T] = spectrogram(quadchirp,100,98,128,fs);
helperCWTTimeFreqPlot(S,T,F,'surf','STFT of Quadratic Chirp','Seconds','Hz')

Получите частотно-временной график этого сигнала с помощью CWT.

[cfs,f] = cwt(quadchirp,fs,'WaveletParameters',[14,200]);
helperCWTTimeFreqPlot(cfs,tquad,f,'surf','CWT of Quadratic Chirp','Seconds','Hz')

CWT с вейвлет производит частотно-временной анализ, очень похожий на STFT.

Частотная и амплитудная модуляция часто происходит в естественных сигналах. Используйте CWT, чтобы получить частотно-временной анализ эхолокационного импульса, излучаемого большой коричневой битой (Eptesicus Fuscus). Интервал дискретизации составляет 7 микросекунд. Используйте bump вейвлет с 32 голосами на октаву. Спасибо Кертису Кондону, Кену Уайту и Аль Фенгу из Центра Бекмана при Университете Иллинойса за данные летучих мышей и разрешение использовать его в этом примере.

load batsignal
t = 0:DT:(numel(batsignal)*DT)-DT;
[cfs,f] = cwt(batsignal,'bump',1/DT,'VoicesPerOctave',32);
helperCWTTimeFreqPlot(cfs,t*1e6,f./1000,'surf','CWT of Bat Echolocation',...
    'Microseconds','kHz')

Получите и постройте график STFT данных о летучих мышах.

[S,F,T] = spectrogram(batsignal,50,48,128,1/DT);
helperCWTTimeFreqPlot(S,T.*1e6,F./1e3,'surf','Bat Echolocation (STFT)',...
    'Microseconds','kHz')

Как для моделируемых, так и для естественных модулированных сигналов, CWT предоставляет результаты, аналогичные STFT.

Обнаружение переходных процессов в колебаниях с помощью CWT

Существуют определенные ситуации во частотно-временном анализе, где CWT может обеспечить более информативное преобразование частота-время, чем короткое преобразование Фурье. Одна из таких ситуаций возникает, когда сигнал повреждается переходными процессами. Внешние виды и исчезновение этих переходных процессов часто имеет физическое значение. Поэтому важно иметь возможность локализовать эти переходные процессы в дополнение к характеристике колебательных компонентов в сигнале. Чтобы симулировать это, создайте сигнал, состоящий из двух синусоид с частотами 150 и 200 Гц. Частота дискретизации составляет 1 кГц. Синусоиды имеют несвязанные временные поддержки. 150-Hz синусоида происходит между 100 и 300 миллисекундами. 200-Hz синусоида происходит от 700 миллисекунд до 1 секунды. Кроме того, существует два переходных процесса на 222 и 800 миллисекундах. Сигнал поврежден шумом.

rng default;
dt = 0.001;
t = 0:dt:1-dt;
addNoise = 0.025*randn(size(t));
x = cos(2*pi*150*t).*(t>=0.1 & t<0.3)+sin(2*pi*200*t).*(t>0.7);
x = x+addNoise;
x([222 800]) = x([222 800 ])+[-2 2];
figure;
plot(t.*1000,x);
xlabel('Milliseconds'); ylabel('Amplitude');

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

subplot(2,1,1)
plot(t(184:264).*1000,x(184:264));
ylabel('Amplitude')
title('Transients')
axis tight;
subplot(2,1,2)
plot(t(760:840).*1000,x(760:840));
ylabel('Amplitude')
axis tight;
xlabel('Milliseconds')

Получите и постройте график CWT с помощью аналитического вейвлета Морле.

figure;
[cfs,f] = cwt(x,1/dt,'amor');
contour(t*1000,f,abs(cfs))
grid on
c = colorbar;
xlabel('Milliseconds')
ylabel('Frequency')
c.Label.String = 'Magnitude';

Аналитический вейвлет Морле показывает более низкую локализацию частоты, чем вейвлет отбойника, но превосходную локализацию времени. Это делает вейвлет Морле лучшим выбором для переходной локализации. Постройте график квадратов коэффициентов мелкой шкалы, чтобы продемонстрировать локализацию переходных процессов.

figure;
wt = cwt(x,1/dt,'amor');
plot(t.*1000,abs(wt(1,:)).^2);
xlabel('Milliseconds'); ylabel('Magnitude');
grid on;
title('Analytic Morlet Wavelet -- Fine Scale Coefficients');
hold off

Вейвлет сжимается, чтобы обеспечить возможность временной локализации переходных процессов с высокой степенью точности при растягивании, чтобы обеспечить возможность частотной локализации колебаний на частотах 150 и 200 Гц.

STFT может локализовать только переходные процессы до ширины окна. Вы должны построить график STFT в децибелах (дБ), чтобы иметь возможность визуализировать переходные процессы.

[S,F,T] = spectrogram(x,50,48,128,1000);
surf(T.*1000,F,20*log10(abs(S)),'edgecolor','none'); view(0,90);
axis tight;
shading interp;
colorbar
xlabel('Time'), ylabel('Hz');
title('STFT')

Переходные процессы появляются в STFT только как широкополосное увеличение степени. Сравните краткосрочные оценки степени, полученные из STFT, до (с центром 183 мс) и после (с центром 223 мс) внешние виды первого переходного процесса.

figure;
plot(F,20*log10(abs(S(:,80))));
hold on;
plot(F,20*log10(abs(S(:,100))),'r');
legend('T = 183 msec','T = 223 msec')
xlabel('Hz');
ylabel('dB');
hold off;

STFT не обеспечивает возможность локализации переходных процессов в той же степени, что и CWT.

Удаление локализованной во времени частотной составляющей с помощью обратного CWT

Создайте сигнал, состоящий из экспоненциально взвешенных синусоид. Существует два компонента 25-Hz - один с центром 0,2 секунды и один с центром 0,5 секунды. Существует два компонента 70-Hz - один с центром 0,2 и один с центром 0,8 секунды. Первые 25-Hz и 70-Hz компоненты совместно происходят во времени.

t = 0:1/2000:1-1/2000;
dt = 1/2000;
x1 = sin(50*pi*t).*exp(-50*pi*(t-0.2).^2);
x2 = sin(50*pi*t).*exp(-100*pi*(t-0.5).^2);
x3 = 2*cos(140*pi*t).*exp(-50*pi*(t-0.2).^2);
x4 = 2*sin(140*pi*t).*exp(-80*pi*(t-0.8).^2);
x = x1+x2+x3+x4;
figure;
plot(t,x)
title('Superimposed Signal')

Получите и отобразите CWT.

cwt(x,2000);
title('Analytic CWT using Default Morse Wavelet');

Удалите компонент 25 Гц, который происходит от приблизительно 0,07 до 0,3 секунд путем обнуления коэффициентов CWT. Используйте обратный CWT (icwt) для восстановления приближения к сигналу.

[cfs,f] = cwt(x,2000);
T1 = .07;  T2 = .33;
F1 = 19;   F2 = 34;
cfs(f > F1 & f < F2, t> T1 & t < T2) = 0;
xrec = icwt(cfs);

Отображение CWT восстановленного сигнала. Компонент начального 25-Hz удаляется.

cwt(xrec,2000)

Постройте график исходного сигнала и реконструкции.

subplot(2,1,1);
plot(t,x);
title('Original Signal');
subplot(2,1,2);
plot(t,xrec)
title('Signal with first 25-Hz component removed');

Наконец, сравните восстановленный сигнал с исходным сигналом без компонента 25-Hz с центром в 0,2 секунды.

y = x2+x3+x4;
figure;
plot(t,xrec)
hold on
plot(t,y,'r--')
legend('Inverse CWT approximation','Original signal without 25-Hz');
hold off

Определение точной частоты через аналитический CWT

Когда вы получаете вейвлет-преобразование синусоиды с помощью аналитического вейвлета, аналитические коэффициенты CWT фактически кодируют частоту.

Чтобы проиллюстрировать это, рассмотрим отоакустическое излучение, полученное из уха человека. Отоакустические выбросы (ОАЭ) выделяются улиткой (внутренним ухом), и их присутствие указывает на нормальный слух. Загрузите и постройте график данных OAE. Данные отбираются с частотой дискретизации 20 кГц.

load dpoae
plot(t.*1000,dpoaets)
xlabel('Milliseconds')
ylabel('Amplitude')

Эмиссию вызывал стимул, начинающийся с 25 миллисекунд и заканчивающийся с 175 миллисекунд. Исходя из экспериментальных параметров, частота выбросов должна составлять 1230 Гц. Получите и постройте график CWT.

cwt(dpoaets,'bump',20000);
xlabel('Milliseconds');

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

[dpoaeCWT,f] = cwt(dpoaets,'bump',20000);
[~,idx1230] = min(abs(f-1230));
cfsOAE = dpoaeCWT(idx1230,:);
plot(t.*1000,abs(cfsOAE));
hold on
AX = gca;
plot([25 25],[AX.YLim(1) AX.YLim(2)],'r')
plot([175 175],[AX.YLim(1) AX.YLim(2)],'r')
xlabel('Milliseconds'), ylabel('Magnitude');
title('CWT Coefficient Magnitudes')

Существует некоторая задержка между началом вызывающего стимула и ОАЭ. Когда вызывающий стимул прекращается, OAE немедленно начинает разрушаться в величине.

Другой способ изоляции излучения - это использование обратной CWT для восстановления частотно-локализованного приближения во временном интервале.

Восстановите частотно-локализованное приближение излучения путем инвертирования CWT в частотной области значений [1150 1350] Гц. Постройте график исходных данных вместе с реконструкцией и маркерами, указывающими на начало и конец вызывающего стимула.

xrec = icwt(dpoaeCWT,f,[1150 1350]);
figure
plot(t.*1000,dpoaets)
hold on;
plot(t.*1000,xrec,'r')
AX = gca;
ylim = AX.YLim;
plot([25 25],ylim,'k')
plot([175 175],ylim,'k')
xlabel('Milliseconds')
ylabel('Amplitude')
title('Frequency-Localized Reconstruction of Emission')

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

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

xdft = fft(xrec);
freq = 0:2e4/numel(xrec):1e4;
xdft = xdft(1:numel(xrec)/2+1);
figure
plot(freq,abs(xdft));
xlabel('Hz'); ylabel('Magnitude')
title('Fourier Transform of CWT-Based Signal Approximation');
[~,maxidx] = max(abs(xdft));
fprintf('The frequency is %4.2f Hz\n',freq(maxidx));

The frequency is 1230.00 Hz

Выводы и дальнейшее чтение

В этом примере вы научились использовать CWT для получения частотно-временного анализа сигнала 1-D с помощью аналитического вейвлета с cwt. Вы видели примеры сигналов, где CWT предоставляет аналогичные результаты для STFT и пример, где CWT может предоставить больше интерпретируемых результатов, чем STFT. Наконец, вы научились восстанавливать локализованные приближения шкалы времени (частоты) к сигналу, используя icwt.

Для получения дополнительной информации о CWT см. документацию Wavelet Toolbox.

Ссылка: Mallat, S. «A Wavelet Tour of Signal Processing: The Sparse Way», Academic Press, 2009.

Приложение

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