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

load mtlb;
w = cwt(mtlb);

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

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

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

cwt(mtlb,Fs)

Получите CWT данных о землетрясении Кобе. Данные являются сейсмографом (вертикальное ускорение, nm/sq.sec) измерения, зарегистрированные в Университете Тасмании, Хобарт, Австралия 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 Гц. Данные производятся на уровне 1 000 Гц. Эти две комплексных экпоненты имеют непересекающуюся поддержку вовремя.

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. Можно проверить это путем выполнения этого скрипта и исследования каждого подграфика с Data Cursor.

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

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

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

То В этом примере показано, как сгенерировать файл 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 Гц. Данные производятся на уровне 1 000 Гц. Эти две синусоиды имеют непересекающуюся поддержку вовремя. Добавьте белый Гауссов шум со стандартным отклонением 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;

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

Загрузите в сигнале и отобразите значение по умолчанию scalogram. Измените палитру в 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';

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

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

Этот пример показывает что, увеличивая продукт пропускной способности времени $$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 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')