Непрерывный и дискретный вейвлет-анализ

Один вопрос, который часто возникает у людей: Что является непрерывным в непрерывном вейвлет-анализе? Ведь вам, предположительно, интересно делать вейвлет-анализ в компьютере, а не карандашом и бумагой, а в компьютере, ничего действительно непрерывного. Когда термин непрерывный вейвлет-анализ используется в научной вычислительной установке, он означает метод вейвлет-анализа с более чем одним вейвлетом на октаву, или удвоением частоты, и где сдвиг между вейвлетами во времени составляет одну выборку. Это обеспечивает результирующее непрерывное вейвлет-преобразование (CWT) имеет два свойства, которые очень полезны в приложениях:

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

Подробное описание обработки вейвлет-сигналов, включая непрерывный вейвлет-анализ со сложными вейвлетами, приведено в [1].

Фильтры или голоса на октаву

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

fb = cwtfilterbank

fb = 
  cwtfilterbank with properties:

      VoicesPerOctave: 10
              Wavelet: 'Morse'
    SamplingFrequency: 1
       SamplingPeriod: []
         PeriodLimits: []
         SignalLength: 1024
      FrequencyLimits: []
        TimeBandwidth: 60
    WaveletParameters: []
             Boundary: 'reflection'

При проверке свойств банка фильтров по умолчанию используется 10 вейвлет-фильтров на октаву (или VoicesPerOctave). Постройте график частотных характеристик вейвлет-фильтров.

freqz(fb)

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

psi = wavelets(fb);
F = centerFrequencies(fb);
figure
plot(real(psi(9,:)))
ylabel('Amplitude')
yyaxis right
plot(real(psi(end-16,:)))
ylabel('Amplitude')
axis tight
S1 = 'Wavelets With Center Frequencies';
S2 = [num2str(F(9),'%1.2f') ' and ' num2str(F(end-16),'%1.2f') ' cycles/sample'];
title({S1 ; S2})
xlabel('Sample')
legend('0.25 c/s Wavelet','0.01 c/s Wavelet')

На графике частотных откликов центральные частоты не превышают 0,45 циклов/выборка. Найдите, сколько вейвлет-фильтров в банке фильтров в одной октаве, или удвоение частоты. Подтвердить, что число равно VoicesPerOctave.

numFilters10 = nnz(F >= 0.2 & F<= 0.4)

numFilters10 = 10

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

fb = cwtfilterbank('VoicesPerOctave',8);
figure
freqz(fb)

Убедитесь, что в банке фильтров имеется восемь фильтров на октаву.

F = centerFrequencies(fb);
numFilters8 = nnz(F >= 0.2 & F<= 0.4)

numFilters8 = 8

Можно прочитать более подробное объяснение различий между непрерывным и дискретным вейвлет-анализом при непрерывных и дискретных вейвлет-преобразованиях.

Логарифмически разнесенные центральные частоты

Одним из аспектов вейвлет-анализа, который люди могут найти немного запутанным, является логарифмический интервал фильтров.

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

subplot(2,1,1)
plot(F)
ylabel('Cycles/Sample')
title('Wavelet Center Frequencies')
grid on
subplot(2,1,2)
semilogy(F)
grid on
ylabel('Cycles/Sample')
title('Semi-Log Scale')

Графики показывают, что частоты вейвлет-центров не расположены линейно, как обычно в случае с другими блоками фильтров. В частности, центральные частоты экспоненциально уменьшаются, так что размер шага между более высокими центральными частотами больше, чем размер шага между более низкими центральными частотами. Почему это имеет смысл для вейвлетов? Напомним, что вейвлеты являются постоянными Q-фильтрами, что означает, что их полосы пропускания пропорциональны их центральным частотам. Если у вас есть банк фильтров, где каждый фильтр имеет одинаковую полосу пропускания, как в кратковременном преобразовании Фурье, или спектрограмме, банк фильтров, то для поддержания постоянного частотного перекрытия между фильтрами необходим постоянный размер шага. Однако при использовании вейвлетов размер шага должен быть пропорционален частоте, подобной полосе пропускания. Для непрерывного вейвлет-анализа наиболее распространенным интервалом является основание 2 ^ (1/NV), где NV - количество фильтров на октаву, повышенное до целых степеней. Для сравнения, интервал, используемый исключительно в дискретном вейвлет-анализе, является основанием 2, увеличенным до целочисленных степеней.

Подробное описание дискретного вейвлет-анализа см. в [2]. Следующая таблица суммирует основные сходства и различия между дискретными и непрерывными вейвлет-методами. Для дискретных методик имена репрезентативных алгоритмов в MATLAB приведены в скобках.

Частотно-временной анализ

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

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

load hyperbolicChirp
figure
plot(t,hyperbolchirp)
axis tight
xlabel('Time (s)')
ylabel('Amplitude')
title('Hyperbolic Chirp')

Гиперболическая чирпа - это сумма $$\mathrm{греха}\frac{\mathrm{(}\mathrm{}}{\mathrm{}\mathrm{}}15security0,8-t$$) и $$\mathrm{}греха\frac{\mathrm{}(}{\mathrm{}\mathrm{}}$$5security0,8-t). Первый компонент активен от 0,1 до 0,68 секунды, а второй компонент от 0,1 до 0,75 секунды. Частота в каждый момент времени или мгновенные частоты этих компонентов $$\frac{\mathrm{составляют}\mathrm{}}{\mathrm{7,5}\mathrm{(}{}^{\mathrm{}}}$$0,8-t) 2 $$\frac{\mathrm{}\mathrm{и}}{\mathrm{}\mathrm{2,5}({}^{\mathrm{}}}$$0,8-t) 2 цикла/с соответственно. Это означает, что мгновенные частоты очень низки $$\mathrm{}$$вблизи t = 0 и быстро увеличиваются по мере приближения времени к 0,8 с.

Получить CWT сигнала частотной модуляции и построить график CWT вместе с истинными мгновенными частотами, показанными пунктирными линиями.

Fs = 1/mean(diff(t));
[cfs,f] = cwt(hyperbolchirp,Fs);
helperHyperbolicChirpPlot(cfs,f,t)

Мы видим, что CWT показывает частотно-временное представление, которое точно фиксирует мгновенные частоты гиперболической чирпы. Теперь проанализируйте тот же сигнал с помощью кратковременного преобразования Фурье, которое использует фильтры постоянной полосы пропускания. Используйте размер окна по умолчанию и перекрытие.

pspectrum(hyperbolchirp,2048,'spectrogram')

Следует отметить, что спектрограмма не способна различать две различные мгновенные частоты, когда частоты низкие, или эквивалентно, когда две мгновенные частоты близки друг к другу. На основе входной длины 2048 спектрограмма вычисляется по умолчанию с использованием длины окна 128 выборок. Это эквивалентно временному разрешению 62,5 мс при частоте дискретизации 2048 Гц. На основе окна Кайзера по умолчанию это приводит к частотному разрешению 41 Гц. Это слишком велико, чтобы дифференцировать мгновенные частоты вблизи t = 0. Что произойдет, если мы уменьшим разрешение частоты приблизительно на 1/2 в надежде на лучшее разрешение более низких частот?

pspectrum(hyperbolchirp,2048,'spectrogram','FrequencyResolution',20, ...
    'OverlapPercent',90)

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

В реальных сигналах высокочастотные события часто имеют короткую продолжительность, в то время как низкочастотные события имеют большую продолжительность. В качестве примера, нагрузка и график сигнала электрокардиограммы человека (ЭКГ). Данные дискретизируют на частоте 180 Гц.

load wecg
tm = 0:1/180:numel(wecg)*1/180-1/180;
plot(tm,wecg)
grid on
axis tight
title('Human ECG')
xlabel('Time (s)')
ylabel('Amplitude')

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

[cfs,f] = cwt(wecg,180);
imagesc(tm,f,abs(cfs))
xlabel('Time (s)')
ylabel('Frequency (Hz)')
axis xy
caxis([0.025 0.25])
title('CWT of ECG Data')

CWT показывает почти стационарные колебания около 30 Гц и 37 Гц, а также переходные события, обозначающие сердцебиение (комплексы QRS). Этот анализ позволяет анализировать оба явления одновременно в одном и том же частотно-временном представлении.

В качестве другого примера высокочастотных переходных процессов, пунктуирующих медленно изменяющиеся компоненты, рассмотрим временной ряд величин солнечного магнитного поля, регистрируемых ежечасно над южным полюсом Солнца космическим аппаратом «Улисс» с 21:00 UT 4 декабря 1993 года по 12:00 UT 24 мая 1994 года. Полное описание этих данных приведено в [2] стр. 218-220.

load solarMFmagnitudes
plot(sm_dates,sm)
axis tight
title('Solar Magnetic Field Magnitudes')
xlabel('Date')
ylabel('Magnitude')

Необработанные временные данные показывают общий примерно линейный тренд с несколькими импульсными событиями, или ударными волнами. Ударные волны возникают, когда быстрый солнечный ветер или быстрый выброс корональной массы настигает медленный солнечный ветер. Во временных рядах мы видим значительные структуры ударной волны, происходящие примерно в следующие даты: 11 февраля, 26 февраля, 10 марта, 3 апреля и 23 апреля.

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

helperSolarMFDataPlot(sm,sm_dates)

CWT фиксирует импульсные события в то же время, когда они происходят во временных рядах. Однако CWT также обнаруживает более низкие частотные характеристики данных, скрытых во временных рядах. Например, начиная до и выходя за пределы структуры ударной волны 11 февраля 1994 года происходит событие низкочастотного устойчивого состояния (около 0,04 циклов/день). Это является результатом события коронального выброса массы (CME), которое, вероятно, происходит в другой солнечной области, которая достигла космического аппарата «Улисс» одновременно с более близкой ударной волной.

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

Локализация переходных процессов

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

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

Сигнал состоит из двух синусоидальных компонентов, которые резко включаются и выключаются с дополнительными «дефектами» через 222 миллисекунды и 800 миллисекунд. Получить CWT сигнала и построить график только с самыми мелкими или самыми высокими коэффициентами центральной частоты. Создайте вторую ось и постройте график исходных временных данных.

cfs = cwt(x,1000,'amor');
plot(t,abs(cfs(1,:)),'linewidth',2)
ylabel('Magnitude')
set(gca,'XTick',[0.1 0.222 0.5 0.7 0.8 1.2])
yyaxis right
plot(t,x,'--')
ylabel('Amplitude')
xlabel('Seconds')
hold off

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

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

Любое частотно-временное преобразование, которое использует фильтры, такие как вейвлеты в случае CWT, или модулированные окна в случае кратковременного преобразования Фурье, обязательно размазывает картинку сигнала по времени и частоте. Неопределенность в локализации энергии сигнала во времени и частоте происходит от разброса фильтров во времени и частоте. Синхроквизирование - это техника, которая пытается компенсировать это размазывание «сдавливанием» преобразования вдоль частотной оси. Чтобы проиллюстрировать это вейвлетами, получаем CWT сигнала, состоящего из сигнала с амплитудной и частотной модуляцией (AM-FM) и синусоиды 18-Hz. Сигнал AM-FM определяется уравнением

load multicompsig
sig = sig1+sig2;
cwt(sig,sampfreq,'amor')

В то время как CWT показывает обоим волну синуса на 18 Гц, а также AM-FM компонент, мы видим, что волна синуса на 18 Гц, конечно, кажется распространенной в частоте. Теперь используйте синхроимпульсирование для резкости оценок частоты.

[sst,f] = wsst(sig,sampfreq);
figure
helperSSTPlot(sst,t,f)

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

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

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

xrec = iwsst(sst,iridge);
helperPlotComponents(xrec,sig1,sig2,t)

Сравнение изменяющегося во времени частотного содержания в двух сигналах

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

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

t = 0:0.001:2;
X = cos(2*pi*10*t).*(t>=0.5 & t<1.1)+ ...
        cos(2*pi*50*t).*(t>= 0.2 & t< 1.4)+0.25*randn(size(t));
Y = sin(2*pi*10*t).*(t>=0.6 & t<1.2)+...
        sin(2*pi*50*t).*(t>= 0.4 & t<1.6)+ 0.35*randn(size(t));
figure
subplot(2,1,1)
plot(t,X)
ylabel('Amplitude')
subplot(2,1,2)
plot(t,Y)
ylabel('Amplitude')
xlabel('Time (s)')

Оба сигнала состоят из двух синусоидальных волн (10 Гц и 50 Гц) в белом шуме. Синусоидальные волны имеют несколько другие опоры времени. Обратите внимание, что компоненты 10 Гц и 50 Гц в Y отстают на 1/4 цикла по отношению к соответствующим компонентам в X. Постройте график вейвлет-когерентности двух сигналов.

   figure
   wcoherence(X,Y,1000,'PhaseDisplayThreshold',0.7)

Оценка вейвлет-когерентности фиксирует, что сигналы сильно коррелируют вблизи 50 и 10 Гц и только в течение правильных временных интервалов. Темные стрелки на графике показывают фазовое соотношение между сигналами. Здесь стрелки указывают прямо вверх, что указывает на 90-градусную фазовую зависимость между двумя сигналами. Именно это соотношение продиктовано косинусными слагаемыми в сигналах.

В качестве реального примера рассмотрим данные региона Эль-Нино 3 и десезонализированный индекс всех индийских осадков с 1871 по конец 2003 года. Данные отбираются ежемесячно. Временной ряд Nino 3 является рекордом ежемесячных аномалий температуры морской поверхности в градусах Цельсия, зарегистрированных от экваториального Тихого океана на долготе от 90 градусов к западу до 150 градусов к западу, и широтах от 5 градусов к северу до 5 градусов к югу.

load ninoairdata
figure
subplot(2,1,1)
plot(datayear,nino)
title('El Nino Region 3 -- SST Anomalies')
ylabel('Degrees Celsius')
xlabel('Year')
axis tight
subplot(2,1,2)
plot(datayear,air)
axis tight
title('Deseasonalized All Indian Rainfall Index')
ylabel('Millimeters')
xlabel('Year')

Вычислите вейвлет-когерентность между двумя временными рядами. Установите порог отображения фазы равным 0,7. Для просмотра периодов в годах укажите интервал выборки как 1/12 года.

figure
wcoherence(nino,air,years(1/12),'PhaseDisplayThreshold',0.7);

На графике показаны локализованные во времени участки сильной когерентности, происходящие в периоды, которые соответствуют типичным циклам El Nino от 2 до 7 лет. График также показывает, что существует приблизительная задержка цикла 3/8-1/2 между двумя временными рядами в эти периоды. Это указывает на то, что периоды потепления моря в соответствии с Эль-Нино, зарегистрированные у берегов Южной Америки, коррелируют с количеством осадков в Индии примерно в 17 000 км, но этот эффект задерживается примерно на 1/2 цикла (1-3,5 года).

Заключения

В этом примере вы узнали:

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

Ссылки

[1] Маллат, С. Г. Вейвлет-тур обработки сигналов: разреженный путь. 3-е ред. Амстердам; Бостон: Elsevier/Академическая пресса, 2009.

[2] Персиваль, Дональд Б. и Эндрю Т. Уолден. Вейвлет-методы для анализа временных рядов. Кембридж: Cambridge University Press, 2000. https://doi.org/10.1017/CBO9780511841040.