Получите MODWTMRA простого сигнала временного ряда и продемонстрируйте идеальную реконструкцию.

Создание сигнала временного ряда

t = 1:10;
x = sin(2*pi*200*t);

Получение MODWT и MODWTMRA и суммирование строк MODWTMRA.

m = modwt(x);
mra = modwtmra(m);
xrec = sum(mra);

Используйте максимум абсолютных значений, чтобы показать, что разница между исходным сигналом и реконструкцией чрезвычайно мала. Наибольшее абсолютное значение - порядка $$\mathrm{}{\mathrm{}}^{\mathrm{}\mathrm{10-25}}$$, что демонстрирует совершенную реконструкцию.

max(abs(x-xrec))

ans = 5.5738e-25

Создайте MRA сигнала ЭКГ до четвертого уровня с помощью db2 вейвлет. Данные взяты из Percival & Walden (2000), стр. 125 (данные первоначально предоставлены Уильямом Константином и Пером Рейнхаллом, Вашингтонский университет). Частота дискретизации сигнала ЭКГ составляет 180 герц.

load wecg;
lev = 4;
wtecg = modwt(wecg,'db2',lev);
mra = modwtmra(wtecg,'db2');

Постройте график формы сигнала ЭКГ и MRA.

t = (0:numel(wecg)-1)/180;
subplot(6,1,1)
plot(t,wecg)
for kk = 2:lev+2
    subplot(6,1,kk)
    plot(t,mra(kk-1,:))
end
xlabel('Time (s)')
set(gcf,'Position',[0 0 500 700])

Создайте анализ множественных решений для данных индекса южного колебания. Период отбора проб составляет один день. Постройте график подробных данных уровня 8, соответствующих шкале $${\mathrm{}}^{28}$$ дней. Детали в этом масштабе фиксируют колебания в масштабе приблизительно одного года.

load soi
wtsoi = modwt(soi);
mrasoi = modwtmra(wtsoi);
plot(mrasoi(8,:))
title('Level 8 Details')

Получение MRA для данных обменного курса Deutsch Mark - доллар США с использованием минимального масштабирования полосы пропускания и вейвлет-фильтров с четырьмя коэффициентами.

load DM_USD;
Lo = [0.4801755, 0.8372545, 0.2269312, -0.1301477];
Hi = qmf(Lo);
wdm = modwt(DM_USD,Lo,Hi);
mra = modwtmra(wdm,Lo,Hi);

Получение MRA для сигнала ЭКГ с помощью 'reflection' обработка границ. Данные взяты из Percival & Walden (2000), стр. 125 (данные первоначально предоставлены Уильямом Константином и Пером Рейнхаллом, Вашингтонский университет).

load wecg;
wtecg = modwt(wecg,'reflection');
mra = modwtmra(wtecg,'reflection');

Показать, что количество столбцов в MRA равно количеству элементов в исходном сигнале.

isequal(size(mra,2),numel(wecg))

ans = logical
   1

Загрузка 23-канальных данных ЭЭГ Espiga3 [3]. Каналы расположены столбчато. Данные дискретизируются на частоте 200 Гц.

load Espiga3

Получение MRA мультисигнала.

w = modwt(Espiga3);
mra = modwtmra(w);

В этом примере показаны различия между функциями MODWT и MODWTMRA. MODWT разделяет энергию сигнала по коэффициентам детализации и коэффициентам масштабирования. MODWTMRA проецирует сигнал на вейвлет-подпространства и масштабирование подпространства.

Выберите sym6 вейвлет. Загрузить и построить график сигнала электрокардиограммы (ЭКГ). Частота дискретизации сигнала ЭКГ составляет 180 герц. Данные взяты из Percival and Walden (2000), стр. 125 (данные первоначально предоставлены Уильямом Константином и Пером Рейнхаллом, Вашингтонский университет).

load wecg
t = (0:numel(wecg)-1)/180;
wv = 'sym6';
plot(t,wecg)
grid on
title(['Signal Length = ',num2str(numel(wecg))])
xlabel('Time (s)')
ylabel('Amplitude')

Возьмите MODWT сигнала.

wtecg = modwt(wecg,wv);

Входные данные являются выборками функции $$f(x$$), оцениваемой в $$$$N-кратных временных точках. Функция может быть выражена в виде линейной комбинации масштабной функции ($$x)$$и $$\mathrm{вейвлет\; (}\mathrm{x)}$$при различных масштабах и $$трансляциях:f\underset{(x\mathrm{)}}{\overset{\mathrm{}}{}}{}_{}{\mathrm{}}^{{}_{\mathrm{}}\mathrm{}}{\mathrm{}}^{{}_{\mathrm{}}}\mathrm{=\sum k=0N-1ck2-J0/2\varphi }(\underset{\mathrm{}}{\overset{{\mathrm{2-J0x-k}}_{\mathrm{)}}}{}}{}_{}$$+∑j=1J0fj ($$x_{\mathrm{),}}\underset{\mathrm{где}\mathrm{}}{\overset{\mathrm{fj}\mathrm{(}}{}}{x}_{)}{\mathrm{}}^{\mathrm{}}{\mathrm{}}^{}$$=∑k=0N-1dj,k2-j/2ψ ($${}_{\mathrm{}}$$2-jx-k) и J0 - количество уровней вейвлет-разложения. Первая сумма является грубой масштабной аппроксимацией $${}_{\mathrm{сигнала}},$$а fj (x) являются деталями на последовательных шкалах. MODWT $$$$возвращает $${}_{}$$коэффициенты числа N $${\{}_{\mathrm{}}ck\}$$и коэффициенты $${}_{\mathrm{детализации\; (J0\; \times \; N)}}$${dj, k} расширения. Каждая строка вwtecg содержит коэффициенты с другой шкалой.

При взятии MODWT сигнала длиной $$N$$существуют уровни $${\mathrm{}}_{\mathrm{}}$$разложения floor (log2 (N)) -many (по умолчанию). Коэффициенты детализации формируются на каждом уровне. Коэффициенты масштабирования возвращаются только для конечного уровня. В этом примере, $$поскольку\mathrm{}\mathrm{N}\mathrm{}\mathrm{=}$$$${}_{\mathrm{}}2048\text{, J0}=\mathrm{}\mathrm{}\mathrm{floor}\mathrm{}\mathrm{(}\mathrm{}log2\mathrm{}\mathrm{(}$$2048)) = 11 и количество строк вwtecg является $${}_{\mathrm{J0}}\mathrm{+}1\mathrm{}=\mathrm{}11\mathrm{}$$+ 1 = 12.

MODWT разделяет энергию по различным масштабам и коэффициентам масштабирования: $${}^{\mathrm{}}\underset{\mathrm{}}{\overset{{}_{\mathrm{}}}{}}{{}_{}}^{\mathrm{}}{{}_{{}_{\mathrm{}}}}^{\mathrm{||X||2=\sum j=1J0||Wj||2+||VJ0||2}}$$, где $$X$$ - входные данные, $${}_{\mathrm{Wj}}$$ - коэффициенты детализации в масштабе $$j$$, и $${}_{{}_{\mathrm{VJ0}}}$$ - коэффициенты масштабирования конечного уровня.

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

energy_by_scales = sum(wtecg.^2,2);
Levels = {'D1';'D2';'D3';'D4';'D5';'D6';'D7';'D8';'D9';'D10';'D11';'A11'};
energy_table = table(Levels,energy_by_scales);
disp(energy_table)

    Levels     energy_by_scales
    _______    ________________

    {'D1' }         14.063     
    {'D2' }         20.612     
    {'D3' }         37.716     
    {'D4' }         25.123     
    {'D5' }         17.437     
    {'D6' }         8.9852     
    {'D7' }         1.2906     
    {'D8' }         4.7278     
    {'D9' }         12.205     
    {'D10'}         76.428     
    {'D11'}         76.268     
    {'A11'}         3.4192     

energy_total = varfun(@sum,energy_table(:,2))

energy_total=table
    sum_energy_by_scales
    ____________________

           298.28       

Подтвердите, что MODWT сохраняет энергию, вычисляя энергию сигнала и сравнивая ее с суммой энергий по всем масштабам.

energy_ecg = sum(wecg.^2);
max(abs(energy_total.sum_energy_by_scales-energy_ecg))

ans = 7.4402e-10

Возьмите MODWTMRA сигнала.

mraecg = modwtmra(wtecg,wv);

MODWTMRA возвращает проекции функции $$f(x$$) на различные вейвлет-подпространства и конечное пространство масштабирования. То есть MODWTMRA возвращает $$\underset{\mathrm{}}{\overset{\mathrm{}}{}}{}_{}{\mathrm{}}^{{}_{\mathrm{}}\mathrm{}}\sum k=0N-1ck2-J0/2\varphi {\mathrm{}}^{({}_{\mathrm{}}}$$2-J0x-k) и $${}_{\mathrm{}}$$$${\mathrm{J0-many}}_{}\{fj$$(x)}, $$\mathrm{оцененные}$$в N-кратных временных точках. Каждая строка вmraecg - проекция $$f(x$$) на другое подпространство. Это означает, что исходный сигнал может быть восстановлен путем сложения всех проекций. Это неверно в случае MODWT. Добавление коэффициентов вwtecg не восстановит исходный сигнал.

Выберите момент времени, добавьте проекции $$f(x$$), вычисленные в этот момент времени, и сравните с исходным сигналом.

time_point = 1000;
abs(sum(mraecg(:,time_point))-wecg(time_point))

ans = 3.0846e-13

Подтвердите, что, в отличие от MODWT, MODWTMRA не является энергосберегающим преобразованием.

energy_ecg = sum(wecg.^2);
energy_mra_scales = sum(mraecg.^2,2);
energy_mra = sum(energy_mra_scales);
max(abs(energy_mra-energy_ecg))

ans = 115.7053

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

plot(t,wecg,'b')
hold on
plot(t,mraecg(4,:),'-')
hold off
grid on
xlim([4 8])
legend('Signal','Projection','Location','northwest')
xlabel('Time (s)')
ylabel('Amplitude')

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

plot(t,wecg,'b')
hold on
plot(t,wtecg(4,:),'-')
hold off
grid on
xlim([4 8])
legend('Signal','Coefficients','Location','northwest')
xlabel('Time (s)')
ylabel('Amplitude')