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

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

Гиперболический щебет является суммой $$\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')

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

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

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

В реальных сигналах высокочастотные события часто кратковременны, в то время как более низкие события частоты более длинны в длительности. Как пример, загрузите и постройте человеческую электрокардиограмму (ECG) сигнал. Данные производятся на уровне 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 Гц, а также переходные события, обозначающие heartbeat (комплексы 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), вероятно, происходящего в различной солнечной области, которая достигла космического корабля Улисса в то же время, что и более близкая ударная волна.

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

Локализуйте переходные процессы

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

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 или модулируемые окна в случае кратковременного преобразования Фурье, обязательно мажет изображение сигнала вовремя и частоты. Неопределенность в локализации энергии сигнала вовремя и частоты прибывает из распространения фильтров вовремя и частоты. Synchrosqueezing является методом, который пытается компенсировать это смазывание путем "сжатия" преобразования вдоль оси частоты. Чтобы проиллюстрировать это с вейвлетами, получите CWT сигнала, состоящего из амплитуды, и частота модулировала (AM-FM) сигнал и синусоиду на 18 Гц. Сигнал AM-FM задан уравнением

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

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

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

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

[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 Гц) в белом шуме. Синусоиды имеют немного отличающиеся поддержки времени. Обратите внимание на то, что компоненты на 50 Гц и на 10 Гц в Y изолированы 1/4 циклом относительно соответствующих компонентов в X. Постройте когерентность вейвлета двух сигналов.

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

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

Для реального примера считайте область Эль-Ниньо 3 данными и deseasonalized Весь индийский индекс Ливня от 1 871 до последних 2003. Данные производятся ежемесячно. Нино 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-to-1/2 задержка цикла между двумя временными рядами в те периоды. Это указывает, что периоды моря, нагревающегося сопоставимый с Эль-Ниньо, зарегистрированным недалеко от берегов Южной Америки, коррелируются с суммами ливня в Индии на расстоянии приблизительно в 17 000 км, но что этот эффект задерживается приблизительно 1/2 цикл (1 к 3,5 годам).

Заключения

В этом примере вы учились:

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

Ссылки

[1] Mallat, S. G. Тур Вейвлета по Обработке сигналов: Разреженный Путь. 3-й редактор Амстердам  ; Бостон: Нажатие Elsevier/Academic, 2009.

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