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

load mtlb;
w = cwt(mtlb);

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

load mtlb
cwt(mtlb,'bump',Fs)

Сравните результат, полученный из CWT, используя вейвлет Морса по умолчанию.

cwt(mtlb,Fs)

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

load kobe

Постройте график данных землетрясения.

plot((1:numel(kobe))./60,kobe)
xlabel('mins')
ylabel('nm/s^2')
grid on
title('Kobe Earthquake Data')

Получить CWT, частоты и конус влияния.

[wt,f,coi] = cwt(kobe,1);

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

cwt(kobe,1)

Получение CWT, периодов времени и конуса воздействия путем указания периода выборки вместо частоты выборки.

[wt,periods,coi] = cwt(kobe,minutes(1/60));

Просмотрите те же данные, указав входной период выборки вместо частоты.

cwt(kobe,minutes(1/60))

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

Fs = 1e3;
t = 0:1/Fs:1;
z = exp(1i*2*pi*32*t).*(t>=0.1 & t<0.3)+2*exp(-1i*2*pi*64*t).*(t>0.7);

Добавьте сложный белый гауссов шум со стандартным отклонением 0,05.

wgnNoise = 0.05/sqrt(2)*randn(size(t))+1i*0.05/sqrt(2)*randn(size(t));
z = z+wgnNoise;

Получение и построение графика cwt с использованием импульса Морса.

cwt(z,Fs)

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

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

Создайте сигнал, состоящий из двух синусоид с непересекающейся поддержкой во времени. Одна синусоида имеет частоту 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)

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

В этом примере показано, как использование набора фильтров 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.

В этом примере показано, как создать файл MEX для выполнения непрерывного вейвлет-преобразования (CWT) с использованием сгенерированного кода CUDA.

cfg = coder.gpuConfig('mex');

t = 0:.001:(1e5*0.001)-0.001;
x = cos(2*pi*32*t).*(t > 10 & t<=50)+cos(2*pi*64*t).*(t >= 60 & t < 90)+ ...
    0.2*randn(size(t));

x = single(x);

codegen cwt -config cfg -args {coder.typeof(single(0),[1 1e5]),coder.typeof('c',[1 inf]),coder.typeof(0)} -report

Code generation successful: To view the report, open('codegen/mex/cwt/html/report.mldatx').

[cfs,f] = cwt_mex(x,'Morse',1e3);
image('XData',t,'YData',f,'CData',abs(cfs),'CDataMapping','scaled')
set(gca,'YScale','log')
axis tight
xlabel('Seconds')
ylabel('Hz')
title('CWT of Two-Tone Signal')

[cfs2,f2] = cwt(x,'Morse',1e3);
max(abs(cfs2(:)-cfs(:)))

ans =

  single

  6.7583e-07

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

Создайте две синусоидальные волны с частотами 32 и 64 Гц. Данные дискретизируют при частоте 1000 Гц. Две синусоидальные волны имеют несвязанную поддержку во времени. Добавьте белый гауссов шум со стандартным отклонением 0,05. Получение и построение графика CWT с использованием импульса Морса по умолчанию.

Fs = 1e3;
t = 0:1/Fs:1;
x = cos(2*pi*32*t).*(t>=0.1 & t<0.3)+sin(2*pi*64*t).*(t>0.7);
wgnNoise = 0.05*randn(size(t));
x = x+wgnNoise;
cwt(x,1000)

График использует логарифмическую частотную ось, поскольку частоты в CWT являются логарифмическими. В MATLAB логарифмические оси находятся в степенях 10 (десятилетия). Вы можете использовать cwtfreqbounds определяют, каковы минимальные и максимальные частоты полосы пропускания для заданной длины сигнала, частоты дискретизации и импульса.

[minf,maxf] = cwtfreqbounds(numel(x),1000);

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

numfreq = 10;
freq = logspace(log10(minf),log10(maxf),numfreq);

Затем переместите маркер перемещения к текущим осям и замените засечки и метки частотной оси на следующие.

AX = gca;
AX.YTickLabelMode = 'auto';
AX.YTick = freq;

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

newplot;
cwt(x,1000);
AX = gca;
freq = 2.^(round(log2(minf)):round(log2(maxf)));
AX.YTickLabelMode = 'auto';
AX.YTick = freq;

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

Загрузите сигнал и отобразите скалограмму по умолчанию. Измените карту цветов на pink(240).

load noisdopp
cwt(noisdopp)
colormap(pink(240))

Возьмем CWT сигнала и получим вейвлет-коэффициенты и частоты.

[cfs,frq] = cwt(noisdopp);

Для эффективного нахождения максимального значения коэффициентов на каждой частоте (уровне) сначала транспонируют абсолютное значение коэффициентов. Найдите минимальное значение на каждом уровне. На каждом уровне вычитайте минимальное значение уровня.

tmp1 = abs(cfs);
t1 = size(tmp1,2);
tmp1 = tmp1';
minv = min(tmp1);
tmp1 = (tmp1-minv(ones(1,t1),:));

Найти максимальное значение на каждом уровне tmp1. Для каждого уровня разделите каждое значение на максимальное значение на этом уровне. Умножьте результат на количество цветов в карте цветов. Установите равным 1 все нулевые записи. Транспонируйте результат.

maxv = max(tmp1);
maxvArray = maxv(ones(1,t1),:);
indx = maxvArray<eps;
tmp1 = 240*(tmp1./maxvArray);
tmp2 = 1+fix(tmp1);
tmp2(indx) = 1;
tmp2 = tmp2';

Просмотрите результат. Теперь значения скалограммы масштабируются на максимальное абсолютное значение на каждом уровне. Частоты отображаются в линейном масштабе.

t = 0:length(noisdopp)-1;
pcolor(t,frq,tmp2);
shading interp
ylabel('Frequency')
title('Scalogram Scaled By Level')
colormap(pink(240))

Этот пример показывает, что увеличение $${}^{\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 на рисунке.

Загрузите образец речи. Данные дискретизируются при 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')