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

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 Гц. Данные производятся на уровне 1 000 Гц. Постройте сигнал.

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) вейвлет. Примените каждый набор фильтров к сигналу и постройте получившийся scalograms. Заметьте, что у Азбуки Морзе (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. Сигнал производится в течение одной секунды на уровне 1 000 Гц. Постройте сигнал.

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 и набор фильтров, чтобы построить scalogram сигнала.

figure
cwt(x,'FilterBank',fb)

Выполните этот скрипт и используйте Data Cursor, чтобы подтвердить, что амплитуды коэффициентов вейвлета чрезвычайно равны амплитудам синусоидальных компонентов.

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

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

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

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

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

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 с теми же пределами частоты. Сгенерируйте scalogram и сравните с scalogram, сгенерированным (3,60) вейвлет Морзе.

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

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

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

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

Этот пример показывает что, увеличивая продукт пропускной способности времени $$P^2$$ из Морзе вейвлет создает вейвлет с большим количеством колебаний под его конвертом. Увеличение $$P^2$$ сужает вейвлет в частоте.

Создайте два набора фильтров. Один набор фильтров имеет TimeBandwidth по умолчанию значение 60. Второй набор фильтров имеет TimeBandwidth значение 10. SignalLength поскольку оба набора фильтров являются 4 096 выборками.

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

Выровняйте peaks 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 нескольких временных рядов.

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

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

Используйте cwt функционируйте и возьмите CWT данных 250 раз. Отобразите используемое прошедшее время.

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 scalogram в подграфике фигуры.

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

load mtlb
cwt(mtlb,Fs)

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

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

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

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

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

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