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

load mtlb;
w = cwt(mtlb);

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

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

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

cwt(mtlb,Fs);

Создайте два комплексных экспоненциала, различных амплитуд, с частотами 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.

Получите CWT данных о землетрясении Кобе. Частота дискретизации составляет 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 Гц и амплитуду, равную 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.

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