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

Загрузите квадратичный щебет, сигнализируют и показывают график его спектрограммы. Частота сигнала начинается на уровне приблизительно 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 речью на октаву. Благодаря Кертису Кондону, Кену Вайту и Аль Фэну из Центра Бекмана в Университете Иллинойса для bat данных и разрешения использовать его в этом примере.

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 bat данных.

[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 Гц находится между 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 с помощью аналитического вейвлета Morlet.

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

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

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 Гц - один сосредоточенный в 0,2 секунды и один сосредоточенный в 0,5 секунды. Существует два компонента на 70 Гц - один сосредоточенный в 0,2 и один сосредоточенный в 0,8 секунды. Первые компоненты на 70 Гц и на 25 Гц co-occur вовремя.

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 Гц удален.

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 Гц, сосредоточенного в 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 на самом деле кодируют частоту.

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

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

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

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

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

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

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

Ссылка: Mallat, S. "Тур вейвлета по обработке сигналов: разреженный путь", Academic Press, 2009.

Приложение

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