Точный частотно-временной анализ

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

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

load quadchirp;
fs = 1000;
[cfs,f] = cwt(quadchirp,'bump',fs);
helperCWTTimeFreqPlot(cfs,tquad,f,'surf','CWT of Quadratic Chirp','Seconds','Hz')

Обратите внимание, что энергия квадратичного щебета размазывается во временной частотной плоскости временной концентрацией вейвлета. Например, если вы фокусируетесь на частотно-временной концентрации величин CWT около 100 Гц, то видите, что она уже, чем наблюдаемая около 500 Гц. Это не является внутренним свойством щебета. Это программный продукт измерительного прибора (CWT). Сравните частотно-временной анализ того же сигнала, полученного с синхронизированным преобразованием.

wsst(quadchirp,1000,'bump')

Синхронизированное преобразование использует информацию фазы в CWT, чтобы уточнить частотно-временной анализ щебета.

Сигнал восстановления из синхронизированного преобразования

Можно восстановить сигнал из синхронизированного преобразования. Это является преимуществом синхронизации по сравнению с другими методами частотно-временного переназначения. Преобразование не обеспечивает идеальную инверсию, но реконструкция стабильна, и результаты обычно довольно хороши. Чтобы продемонстрировать это, нагрузку, построение графика и воспроизведение записи женщины-оратора с надписью «Я видел овец».

load wavsheep;
plot(tsh,sheep)
title(' "I saw the sheep." ');
xlabel('Time (secs)'); ylabel('Amplitude');
grid on;
hplayer = audioplayer(sheep,fs);
play(hplayer);

Постройте график синхронизированного преобразования речевой выборки.

[sst,f] = wsst(sheep,fs,'bump');
contour(tsh,f,abs(sst));
title('Wavelet Synchrosqueezed Transform');
xlabel('Time (secs)'); ylabel('Hz');
grid on;

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

xrec = iwsst(sst,'bump');
plot(tsh,sheep)
hold on;
plot(tsh,sheep,'r');
xlabel('Time (secs)'); ylabel('Amplitude');
grid on;
title('Reconstruction From Synchrosqueezed Transform');
legend('Original Waveform','Synchrosqueezed Reconstruction');
sprintf('The maximum sample-by-sample difference is %1.3f', ...
    norm(abs(xrec-sheep),Inf))
hold off;

ans =

    'The maximum sample-by-sample difference is 0.004'

Воспроизведите восстановленный сигнал и сравните с оригиналом.

hplayerRecon = audioplayer(xrec,fs);
play(hplayerRecon)

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

Идентифицируйте гребни и восстановите режимы

Частотно-временной гребень определяется локальными максимумами частотно-временного преобразования. Поскольку синхронизированное преобразование концентрирует колебательные режимы в узкой области частотно-временной плоскости и является инвертируемым, можно восстановить отдельные режимы путем:

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

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

Загрузите данные и постройте график синхронизируемого преобразования.

load batsignal;
time = 0:DT:(numel(batsignal)*DT)-DT;
[sst,f] = wsst(batsignal,1/DT);
contour(time.*1e6,f./1000,abs(sst));
grid on;
xlabel('microseconds'); ylabel('kHz');
title('Bat Echolocation Pulse');

Обратите внимание, что существует два модулированных режима, которые отслеживают изогнутые пути в частотно-временной плоскости. Попытка разделить эти компоненты с помощью обычной полосно-пропускающей фильтрации не работает, потому что фильтр должен работать изменяющимся во времени образом. Например, обычный фильтр с полосой пропускания от 18 до 40 кГц будет захватывать энергию в самом раннем импульсе, но также будет захватывать энергию от последующего импульса.

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

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

[fridge,iridge] = wsstridge(sst,5,f,'NumRidges',2);
hold on;
plot(time.*1e6,fridge./1e3,'k--','linewidth',2);
hold off;

Гребни следуют изменяющемуся во времени характеру модулируемых импульсов. Восстановите режимы сигнала путем инвертирования синхронизированного преобразования.

xrec = iwsst(sst,iridge);
subplot(2,1,1)
plot(time.*1e6,batsignal); hold on;
plot(time.*1e6,xrec(:,1),'r');
grid on; ylabel('Amplitude');
title('Bat Echolocation Pulse with Reconstructed Modes');
legend('Original Signal','Mode 1','Location','SouthEast');
subplot(2,1,2);
plot(time.*1e6,batsignal); hold on;
grid on;
plot(time.*1e6,xrec(:,2),'r');
xlabel('microseconds'); ylabel('Amplitude');
legend('Original Signal','Mode 2','Location','SouthEast');

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

Срок штрафа при извлечении хребта

В предыдущем примере извлечения режима использовался термин штрафа в извлечении хребта без объяснения причин.

Когда вы извлекаете несколько гребней, или у вас есть один модулируемый компонент аддитивного шума, важно использовать штрафной термин в извлечении гребня. Штрафной термин служит, чтобы предотвратить переходы частоты, когда область наивысшей энергии в плоскости частота-перемещается.

Чтобы продемонстрировать это, рассмотрим двухкомпонентный сигнал, состоящий из амплитуды и частотно-модулированного сигнала плюс синусоиды. Дискретизация сигнала осуществляется на частоте 1000 Гц. Частота синусоид составляет 18 Гц. Сигнал AM-FM определяется:

Загрузите сигнал и получите синхронизированное преобразование.

load multicompsig;
sig = sig1+sig2;
[sst,f] = wsst(sig,sampfreq);
figure;
contour(t,f,abs(sst));
grid on;
title('Synchrosqueezed Transform of AM-FM Signal Plus Sine Wave');
xlabel('Time (secs)'); ylabel('Hz');

Первая попытка извлечь два гребня из синхроскейзированного преобразования без использования штрафа.

[fridge,iridge] = wsstridge(sst,f,'NumRidges',2);
hold on;
plot(t,fridge,'k--','linewidth',2);

Вы видите, что гребень скачет между сигналом AM-FM и синусоидой, так как область наивысшей энергии во временной частотной плоскости изменяется между этими двумя сигналами. Добавьте срок штрафа 5 к извлечению хребта. В этом случае переходы частоты штрафуются в 5 раз больше расстояния между частотами в терминах интервалов (не фактическая частота в герцах).

[fridge,iridge] = wsstridge(sst,5,f,'NumRidges',2);
figure;
contour(t,f,abs(sst));
grid on;
title('Synchrosqueezed Transform of AM-FM Signal Plus Sine Wave');
xlabel('Time (secs)'); ylabel('Hz');
hold on;
plot(t,fridge,'k--','linewidth',2);
hold off;

При штрафном сроке два режима колебаний изолированы в двух отдельных гребнях. Восстановите режимы вдоль частотно-временных гребней из синхронизированного преобразования.

xrec = iwsst(sst,iridge);

Постройте график восстановленных режимов вместе с исходными сигналами для сравнения.

subplot(2,2,1)
plot(t,xrec(:,1)); grid on;
ylabel('Amplitude');
title('Reconstructed Mode');
ylim([-1.5 1.5]);
subplot(2,2,2);
plot(t,sig2); grid on;
ylabel('Amplitude');
title('Original Component');
ylim([-1.5 1.5]);
subplot(2,2,3);
plot(t,xrec(:,2)); grid on;
xlabel('Time (secs)'); ylabel('Amplitude');
title('Reconstructed Mode');
ylim([-1.5 1.5]);
subplot(2,2,4);
plot(t,sig1); grid on;
xlabel('Time (secs)'); ylabel('Amplitude');
title('Original Component');
ylim([-1.5 1.5]);

Заключения

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

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

Ссылки

Daubechies, I., Lu, J., and Wu, H-T. Synchrosqueezed wavelet transforms: a empirical mode decomposition-like tool (неопр.) (недоступный инструмент). Апл. Компут. Гармонический анал., 30 (2): 243-261, 2011.

Thakur, G., Brevdo, E., Fuckar, N.S. и Wu H-T. «Алгоритм синхронизации для изменяющегося во времени спектрального анализа: свойства робастности и новые приложения палеоклиматов». Обработка сигналов, 93 (5), 1079-1094, 2013.

Meignen, S., Oberlin, T. and McLaughlin, S. «Новый алгоритм многокомпонентного анализа сигналов, основанный на синхронизации: с приложением к дискретизации сигналов и шумоподавления». Транзакции IEEE по обработке сигналов, том 60, № 12, стр. 5787-5798, 2012.

Приложение

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