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

Один из вопросов, который люди часто имеют: Что непрерывно в непрерывном вейвлет? Ведь вам, предположительно, интересно делать вейвлет в компьютере и не карандашом и бумагой, а в компьютере ничего по-настоящему непрерывного. Когда термин непрерывный вейвлет-анализ используется в научной вычислительной настройке, это означает метод вейвлет-анализа с более чем одним вейвлет на октаву или удвоение частоты, и где сдвиг между вейвлетами во времени является одной выборкой. Это обеспечивает полученное непрерывное вейвлет (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

Более подробное объяснение различий между непрерывным и дискретным вейвлетом анализом можно прочитать в Continuous и Discrete Вейвлета Transforms.

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

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

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

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

Графики показывают, что частоты центра вейвлета не разнесены линейно, как обычно в случае с другими банками фильтров. В частности, центральные частоты экспоненциально уменьшаются, так что размер шага между более высокими центральными частотами больше, чем размер шага между более низкими центральными частотами. Почему это имеет смысл для вейвлетов? Напомним, что вейвлеты являются фильтрами constant-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')

Гиперболический щебет является суммой $$\sin \left(\frac{15\pi }{0.8-t}\right)$$ и $$\sin \left(\frac{5\pi }{0.8-t}\right)$$. Первый компонент активен от 0,1 до 0,68 секунд, а второй - от 0,1 до 0,75 секунд. Частота в каждый момент времени или мгновенные частоты этих компонентов $$\frac{7.5}{(0.8-t{)}^2}$$ и $$\frac{2.5}{(0.8-t{)}^2}$$ циклов/сек соответственно. Это означает, что мгновенные частоты очень низки около $$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] pp. 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 показывают peaks, где синусоидальные компоненты включают и отключают, а также локализуют дефекты синусоидов в 222 миллисекундах и 800 миллисекундах. Поскольку коэффициенты CWT имеют такое же разрешение по времени, как и данные, часто полезно использовать мелкомасштабные величины CWT в порядок для обнаружения переходных изменений в данных. Чтобы наиболее точно локализовать эти переходные процессы, используйте коэффициенты самой мелкой шкалы (самая высокая частота).

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

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

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

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

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

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

[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 Эль-Нино и десесонализировал All Indian Redual Index с 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);

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

Заключения

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

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

Ссылки

[1] Mallat, S. G. A Wavelet Tour of Signal Processing: The Sparse Way. 3-й эд. Амстердам; Бостон: Elsevier/Академическая пресса, 2009.

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