Нагружают и визуализируют нестационарный непрерывный сигнал, состоящий из синусоидальных волн с явным изменением частоты. Сигнал дискретизируется на частоте 250 Гц.

fs = 250;
load nonstatdistinct
t = (0:length(nonstatdistinct)-1)/fs;
plot(t,nonstatdistinct)
xlabel('Time (s)')
ylabel('Signal')
axis tight

Использовать ewt для получения анализа множественного разрешения (MRA) сигнала.

mra = ewt(nonstatdistinct);

Использовать компоненты MRA с hht и построить график гильбертова спектра.

hht(mra,fs)

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

Создать нестационарный непрерывный сигнал, состоящий из синусоидальных волн с явным изменением частоты. Сигнал дискретизируется на частоте 1000 Гц.

Fs = 1000;
t = 0:1/Fs:4;
x1 = sin(2*pi*50*t) + sin(2*pi*200*t);
x2 = sin(2*pi*25*t) + sin(2*pi*100*t) + sin(2*pi*250*t);
x = [x1 x2] + 0.1*randn(1,length(t)*2);
t1 = (0:length(x)-1)/Fs;
plot(t1,x)
xlabel('Time (s)')
ylabel('Amplitude')
title('Signal')

Использовать ewt и получают MRA сигнала. Отображение нормированных пиковых частот, идентифицированных в сигнале, и приблизительных частотных диапазонов набора фильтров. Так как частоты находятся в циклах на выборку, нормализуйте по частоте выборки. Следует отметить, что пиковые частоты соответствуют частотам синусоидальных волн.

[mra,~,wfb,info] = ewt(x);
Fs*info.PeakFrequencies

ans = 5×1

  249.9375
  200.0750
  100.1000
   50.1125
   25.1187

Fs*info.FilterBank.Passbands 

ans = 5×2

  223.6941  500.0000
  141.5896  223.6941
   70.8573  141.5896
   35.4911   70.8573
         0   35.4911

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

f = 0:Fs/length(x):Fs-1/length(x);
plot(f,wfb)
ylabel('Magnitude')
grid on
yyaxis right
plot(f,abs(fft(x)),'k--','linewidth',1.5)
ylabel('Magnitude')
xlabel('Hz')

Поскольку эмпирические вейвлеты образуют плотный кадр Парсеваля, набор фильтров анализа равен набору фильтров синтеза. Следовательно, возведение в квадрат величин на каждой частоте, суммированной по фильтрам, равно 1. Если бы сумма не была равна 1, идеальная реконструкция была бы невозможна.

Загрузите сигнал ЭКГ. Сигнал дискретизируется на частоте 180 Гц.

load wecg

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

mp = 4;
[mra,cfs] = ewt(wecg,'MaxNumPeaks',mp);

Постройте график сигнала и компонентов MRA.

fs = 180;
subplot(mp+1,1,1)
t = (0:length(wecg)-1)/fs;
plot(t,wecg)
title('MRA of Signal')
ylabel('Signal')
axis tight
for k=1:mp
    subplot(mp+1,1,k+1)
    plot(t,mra(:,k))
    ylabel(['MRA ',num2str(k)])
    axis tight
end
xlabel('Time (s)')

Убедитесь, что суммирование компонентов MRA приводит к точной реконструкции сигнала.

max(abs(wecg-sum(mra,2)))

ans = 8.8818e-16

Проверка сохранения энергии коэффициентов анализа ЭВТ.

cfsenergy = sum(sum(abs(cfs).^2));
[cfsenergy norm(wecg,2)^2]

ans = 1×2

  298.2759  298.2759