Создайте набор фильтров непрерывного вейвлет-преобразования.

fb = cwtfilterbank

fb = 
  cwtfilterbank with properties:

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

Постройте график частотной характеристики величины.

freqz(fb)

Создайте две синусоидальные волны с частотами 16 и 64 Гц. Данные дискретизируют при частоте 1000 Гц. Постройте график сигнала.

Fs = 1e3;
t = 0:1/Fs:1-1/Fs;
x = cos(2*pi*64*t).*(t>=0.1 & t<0.3)+sin(2*pi*16*t).*(t>=0.5 & t<0.9);
plot(t,x)
title('Signal')

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

fb = cwtfilterbank('SignalLength',numel(t),'SamplingFrequency',Fs);
freqz(fb)
title('Frequency Responses — Morse (3,60) Wavelet')

Аналитический вейвлет Морса (3,60) является вейвлетом по умолчанию в наборе фильтров. Вейвлет имеет произведение временной полосы, равное 60. Создайте второй набор фильтров, идентичный первому набору фильтров, но вместо этого используйте аналитический вейвлет Морса (3,5). Постройте график частотных откликов вейвлетов во втором наборе фильтров.

fb3x5 = cwtfilterbank('SignalLength',numel(t),'SamplingFrequency',Fs,...
    'TimeBandwidth',5);
figure
freqz(fb3x5)
title('Frequency Responses — Morse (3,5) Wavelet')

Обратите внимание, что частотные характеристики шире, чем в первом наборе фильтров. Импульс Морса (3,60) лучше локализован по частоте, чем импульс Морса (3,5). Примените каждый набор фильтров к сигналу и постройте график результирующих скалограмм. Обратите внимание, что импульс Морса (3,60) имеет лучшее частотное разрешение, чем импульс Морса (3,5).

figure
cwt(x,'FilterBank',fb)
title('Magnitude Scalogram — Morse (3,60)')

figure
cwt(x,'FilterBank',fb3x5)
title('Magnitude Scalogram — Morse (3,5)')

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

Создайте сигнал, состоящий из двух синусоид с непересекающейся поддержкой во времени. Одна синусоида имеет частоту 32 Гц и амплитуду, равную 1. Другая синусоида имеет частоту 64 Гц и амплитуду, равную 2. Сигнал дискретизируется в течение одной секунды при частоте 1000 Гц. Постройте график сигнала.

frq1 = 32;
amp1 = 1;
frq2 = 64;
amp2 = 2;

Fs = 1e3;
t = 0:1/Fs:1;
x = amp1*sin(2*pi*frq1*t).*(t>=0.1 & t<0.3)+amp2*sin(2*pi*frq2*t).*(t>0.6 & t<0.9);

plot(t,x)
grid on
xlabel('Time (sec)')
ylabel('Amplitude')
title('Signal')

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

fb = cwtfilterbank('SignalLength',numel(x),'SamplingFrequency',Fs,...
    'FrequencyLimits',[20 100]);
figure
freqz(fb)

Использовать cwt и банк фильтров для построения графика скалограммы сигнала.

figure
cwt(x,'FilterBank',fb)

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

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

Обобщённые вейвлеты Морса - семейство точно аналитических вейвлетов. Вейвлеты Морзе имеют два параметра, симметрию и произведение временной полосы. Эти параметры можно изменять для получения аналитических вейвлетов с различными свойствами и поведением. Дополнительные сведения см. в разделе Вейвлеты Морзе и ссылки на них.

Загрузить сейсмографические данные, зарегистрированные во время землетрясения в Кобе в 1995 году. Данные являются сейсмографическими (вертикальное ускорение, нм/кв.сек) измерениями, зарегистрированными в Университете Тасмании, Хобарт, Австралия, 16 января 1995 начиная с 20:56:51 (GMT) и продолжающимися в течение 51 минут 1 секундными интервалами. Создайте банк фильтров CWT с настройками по умолчанию, которые можно применить к данным. Используйте банк фильтров для создания скалограммы.

load kobe
fb = cwtfilterbank('SignalLength',numel(kobe),'SamplingFrequency',1);
cwt(kobe,'FilterBank',fb)

Величина вейвлет-коэффициентов велика в диапазоне частот от 10 мГц до 100 мГц. Создайте новый банк фильтров с ограничениями по частоте, установленными для этих значений. Создайте скалограмму.

fb2 = cwtfilterbank('SignalLength',numel(kobe),'SamplingFrequency',1,...
    'FrequencyLimits',[1e-2 1e-1]);
cwt(kobe,'FilterBank',fb2)
title('Default (3,60) Morse')

По умолчанию cwtfilterbank использует (3,60) импульс Морса. Создайте набор фильтров, используя аналитический вейвлет Morlet с теми же пределами частоты. Создайте скалограмму и сравните со скалограммой, сгенерированной (3,60) вейвлетом Морса.

fbMorlet = cwtfilterbank('SignalLength',numel(kobe),'SamplingFrequency',1,...
    'FrequencyLimits',[1e-2 1e-1],...
    'Wavelet','amor');
cwt(kobe,'FilterBank',fbMorlet)
title('Analytic Morlet')

Вейвлет Морле не так хорошо локализован по частоте, как (3,60) вейвлет Морса. Однако, изменяя произведение временной полосы пропускания, можно создать вейвлет Морса со свойствами, аналогичными вейвлету Морле.

Создайте набор фильтров, используя вейвлет Морса со значением временной полосы 30 [2], с пределами частоты, как описано выше. Создайте скалограмму данных сейсмографа. Обратите внимание, что частота мазка почти идентична результатам Морле.

fbMorse = cwtfilterbank('SignalLength',numel(kobe),'SamplingFrequency',1,...
    'FrequencyLimits',[1e-2 1e-1],...
    'TimeBandwidth',30);
cwt(kobe,'FilterBank',fbMorse)
title('(3,30) Morse')

Теперь изучите вейвлеты, связанные с fbMorlet и fbMorse фильтровать банки. Из обоих банков фильтров получают частоты вейвлет-центров, частотные характеристики фильтра и вейвлеты временной области. Подтвердите, что центральные частоты почти идентичны.

cfMorlet = centerFrequencies(fbMorlet);
[frMorlet,fMorlet] = freqz(fbMorlet);
[wvMorlet,tMorlet] = wavelets(fbMorlet);
cfMorse = centerFrequencies(fbMorse);
[frMorse,fMorse] = freqz(fbMorse);
[wvMorse,tMorse] = wavelets(fbMorse);

disp(['Number of Center Frequencies: ',num2str(length(cfMorlet))]);

Number of Center Frequencies: 34

disp(['Maximum difference: ',num2str(max(abs(cfMorlet-cfMorse)))]);

Maximum difference: 2.7756e-17

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

wv = 13;
figure
plot(fMorlet,frMorlet(wv,:));
hold on
plot(fMorse,frMorse(wv,:));
grid on
title('Frequency Response')
xlabel('Frequency')
ylabel('Amplitude')
legend('Morlet','(3,30) Morse')

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

figure
subplot(2,1,1)
plot(tMorlet,real(wvMorlet(wv,:)))
hold on
plot(tMorse,real(wvMorse(wv,:)))
grid on
title('Real')
legend('Morlet','(3,30) Morse')
xlim([-100 100])
subplot(2,1,2)
plot(tMorlet,imag(wvMorlet(wv,:)))
hold on
plot(tMorse,imag(wvMorse(wv,:)))
grid on
title('Imaginary')
legend('Morlet','(3,30) Morse')
xlim([-100 100])

Этот пример показывает, что увеличение $${}^{\mathrm{P2\; произведения\; временной\; полосы}}$$Морса создает вейвлет с большим количеством колебаний под его огибающей. Увеличение $${}^{\mathrm{P2}}$$ сужает вейвлет по частоте.

Создайте два банка фильтров. Один банк фильтров имеет значение по умолчанию TimeBandwidth значение 60. Второй банк фильтров имеет TimeBandwidth значение 10. SignalLength для обоих банков фильтров - 4096 образцов.

sigLen = 4096;
fb60 = cwtfilterbank('SignalLength',sigLen);
fb10 = cwtfilterbank('SignalLength',sigLen,'TimeBandwidth',10);

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

[psi60,t] = wavelets(fb60);
[psi10,~] = wavelets(fb10);

Используйте scales функция для поиска основного вейвлета для каждого набора фильтров.

sca60 = scales(fb60);
sca10 = scales(fb10);
[~,idx60] = min(abs(sca60-1));
[~,idx10] = min(abs(sca10-1));
m60 = psi60(idx60,:);
m10 = psi10(idx10,:);

Так как квота временной полосы больше для fb60 банка фильтров, проверьте m60 вейвлет имеет больше колебаний под своей оболочкой, чем m10 вейвлет.

subplot(2,1,1)
plot(t,abs(m60))
grid on
hold on
plot(t,real(m60))
plot(t,imag(m60))
xlim([-30 30])
legend('abs(m60)','real(m60)','imag(m60)')
title('TimeBandwidth = 60')
subplot(2,1,2)
plot(t,abs(m10))
grid on
hold on
plot(t,real(m10))
plot(t,imag(m10))
xlim([-30 30])
legend('abs(m10)','real(m10)','imag(m10)')
title('TimeBandwidth = 10')

Выровнять пики m60 и m10 частотные характеристики величины. Проверьте частотную характеристику m60 вейвлет уже, чем частотный отклик для m10 вейвлет.

cf60 = centerFrequencies(fb60);
cf10 = centerFrequencies(fb10);

m60cFreq = cf60(idx60);
m10cFreq = cf10(idx10);

freqShift = 2*pi*(m60cFreq-m10cFreq);
x10 = m10.*exp(1j*freqShift*(-sigLen/2:sigLen/2-1));

figure
plot([abs(fft(m60)).' abs(fft(x10)).'])
grid on
legend('Time-bandwidth = 60','Time-bandwidth = 10')
title('Magnitude Frequency Responses')

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

Загрузить сейсмографические данные, зарегистрированные во время землетрясения в Кобе в 1995 году. Данные являются сейсмографическими (вертикальное ускорение, нм/кв.сек) измерениями, зарегистрированными в Университете Тасмании, Хобарт, Австралия, 16 января 1995 начиная с 20:56:51 (GMT) и продолжающимися в течение 51 минут 1 секундными интервалами. Создайте банк фильтров CWT, который можно применить к данным.

load kobe
fb = cwtfilterbank('SignalLength',numel(kobe),'SamplingFrequency',1);

Используйте cwt и 250 раз возьмем CWT данных. Отображение затраченного времени.

num = 250;
tic;
for k=1:num
    cfs = cwt(kobe);
end
toc

Elapsed time is 6.551628 seconds.

Теперь используйте wt объектная функция банка фильтров для получения CWT данных. Подтвердите, что использование банка фильтров выполняется быстрее.

tic;
for k=1:num
    cfs = wt(fb,kobe);
end
toc

Elapsed time is 3.782376 seconds.

В этом примере показано, как нарисовать скалограмму CWT на рисунке.

Загрузите образец речи. Данные дискретизируются при 7418 Гц. Постройте график скалограммы CWT по умолчанию.

load mtlb
cwt(mtlb,Fs)

Получить непрерывное вейвлет-преобразование сигнала и частоты CWT.

[cfs,frq] = cwt(mtlb,Fs);

cwt устанавливает оси времени и частоты в скалограмме. Создайте вектор, представляющий время выборки.

tms = (0:numel(mtlb)-1)/Fs;

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

figure
subplot(2,1,1)
plot(tms,mtlb)
axis tight
title('Signal and Scalogram')
xlabel('Time (s)')
ylabel('Amplitude')
subplot(2,1,2)
surface(tms,frq,abs(cfs))
axis tight
shading flat
xlabel('Time (s)')
ylabel('Frequency (Hz)')
set(gca,'yscale','log')