EMD分解

clc

clear all

close all

% [x, Fs] = wavread('Hum.wav');

% Ts = 1/Fs;

% x = x(1:6000);

Ts = 0.001;

Fs = 1/Ts;

t=0:Ts:1;

x = sin(2*pi*10*t) + sin(2*pi*50*t) + sin(2*pi*100*t) + 0.1*randn(1, length(t));

imf = emd(x);

plot_hht(x,imf,1/Fs);

k = 4;

y = imf{k};

N = length(y);

t = 0:Ts:Ts*(N-1);

[yenvelope, yfreq, yh, yangle] = HilbertAnalysis(y, 1/Fs);

yModulate = y./yenvelope;

[YMf, f] = FFTAnalysis(yModulate, Ts);

Yf = FFTAnalysis(y, Ts);

figure

subplot(321)

plot(t, y)

title(sprintf('IMF%d', k))

xlabel('Time/s')

ylabel(sprintf('IMF%d', k));

subplot(322)

plot(f, Yf)

title(sprintf('IMF%d的频谱', k))

xlabel('f/Hz')

ylabel('|IMF(f)|');

subplot(323)

plot(t, yenvelope)

title(sprintf('IMF%d的包络', k))

xlabel('Time/s')

ylabel('envelope');

subplot(324)

plot(t(1:end-1), yfreq)

title(sprintf('IMF%d的瞬时频率', k))

xlabel('Time/s')

ylabel('Frequency/Hz');

subplot(325)

plot(t, yModulate)

title(sprintf('IMF%d的调制信号', k))

xlabel('Time/s')

ylabel('modulation');

subplot(326)

plot(f, YMf)

title(sprintf('IMF%d调制信号的频谱', k))

xlabel('f/Hz')

ylabel('|YMf(f)|');

findpeaks.m文件

function n = findpeaks(x)

% Find peaks. 找极大值点，返回对应极大值点的坐标

n = find(diff(diff(x) > 0) < 0); % 相当于找二阶导小于0的点u = find(x(n+1) > x(n));

n(u) = n(u)+1; % 加1才真正对应极大值点

% 图形解释上述过程

% figure

% subplot(611)

% x = x(1:100);

% plot(x, '-o')

% grid on

%

% subplot(612)

% plot(1.5:length(x), diff(x) > 0, '-o')

% grid on

% axis([1,length(x),-0.5,1.5])

%

% subplot(613)

% plot(2:length(x)-1, diff(diff(x) > 0), '-o')

% grid on

% axis([1,length(x),-1.5,1.5])

%

% subplot(614)

% plot(2:length(x)-1, diff(diff(x) > 0)<0, '-o')

% grid on

% axis([1,length(x),-1.5,1.5])

%

% n = find(diff(diff(x) > 0) < 0);

% subplot(615)

% plot(n, ones(size(n)), 'o')

% grid on

% axis([1,length(x),0,2])

%

% u = find(x(n+1) > x(n));

% n(u) = n(u)+1;

% subplot(616)

% plot(n, ones(size(n)), 'o')

% grid on

% axis([1,length(x),0,2])

plot_hht.m文件

function plot_hht(x,imf,Ts)

% Plot the HHT.

% :: Syntax

% The array x is the input signal and Ts is the sampling period. % Example on use: [x,Fs] = wavread('Hum.wav');

% plot_hht(x(1:6000),1/Fs);

% Func : emd