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

Загрузите квадратичный чирп-сигнал и покажите график его спектрограммы. Частота сигнала начинается приблизительно с 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 микросекунд. Используйте импульсный импульс с 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 Гц происходит от 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) испускаются улиткой (внутренним ухом), и их присутствие свидетельствует о нормальном слухе. Загрузите и постройте график данных 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')

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

Другим способом изолирования излучения является использование обратного 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.

Справка: Маллат, С. «Вейвлет-тур обработки сигналов: разреженный путь», Академическая пресса, 2009.

Приложение

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