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

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)

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

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

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

Загрузка данных сейсмографа, зарегистрированных во время землетрясения в Кобе в 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) вейвлет Морса. Создайте банк фильтров с помощью аналитического вейвлета Морле с теми же пределами частоты. Сгенерируйте скалограмму и сравните с скалограммой, сгенерированной (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])

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

Создайте два банка фильтров. Один банк фильтров имеет значение по умолчанию TimeBandwidth значение 60. Вторая группа фильтров имеет TimeBandwidth значение 10. The 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 filter bank, проверьте 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 нескольких временных рядов.

Загрузка данных сейсмографа, зарегистрированных во время землетрясения в Кобе в 1995 году. Данные представляют собой сейсмографические (вертикальное ускорение, нм/кв.с) измерения, зарегистрированные в Университете Тасмании, Хобарт, Австралия, 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 в подграфике рисунка.

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

load mtlb
cwt(mtlb,Fs)

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

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

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