Резкость частотно-временного анализа

Синхроквизирование может компенсировать расширение во времени и частоте, вызванное линейными преобразованиями, такими как кратковременные преобразования Фурье и непрерывные вейвлет-преобразования (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);

Вы видите, что гребень скачет между АМ-ЧМ сигналом и синусоидальной волной, так как область наибольшей энергии во временной частотной плоскости изменяется между двумя сигналами. Добавить штрафной срок 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]);

Заключения

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

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

Ссылки

Daubechies, I., Lu, J. и Wu, H-T. «Синхронизированные вейвлет-преобразования: эмпирический инструмент, похожий на декомпозицию». Апл. Компут. Гармонический Анал., 30 (2): 243-261, 2011.

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

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

Приложение

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