emd分解 plot_hht 怎么出错了？
源程序clc;
clear all;
close all;
% = 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);
一运行到
plot_hht(x,imf,1/Fs);
??? Cell contents reference from a non-cell array object.
Error in ==> plot_hht at 12
b(k) = sum(imf{k}.*imf{k});
怎么解决？
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: = wavread('Hum.wav');
% plot_hht(x(1:6000),1/Fs);
% Func : emd
% imf = emd(x);
for k = 1:length(imf)
b(k) = sum(imf{k}.*imf{k});
th = unwrap(angle(hilbert(imf{k})));% 相位
d{k} = diff(th)/Ts/(2*pi); % 瞬时频率
end
= sort(-b);
b = 1-b/max(b); % 后面绘图的亮度控制
% Hilbert瞬时频率图
N = length(x);
c = linspace(0,(N-2)*Ts,N-1); % 0:Ts:Ts*(N-2)
for k = v(1:2) % 显示能量最大的两个IMF的瞬时频率
figure
plot(c,d{k});
xlim();
ylim();
xlabel('Time/s')
ylabel('Frequency/Hz');
title(sprintf('IMF%d', k))
end
% 显示各IMF
M = length(imf);
N = length(x);
c = linspace(0,(N-1)*Ts,N); % 0:Ts:Ts*(N-1)
for k1 = 0:4:M-1
figure
for k2 = 1:min(4,M-k1)
subplot(4,2,2*k2-1)
plot(c,imf{k1+k2})
set(gca,'FontSize',8,'XLim',);
title(sprintf('第%d个IMF', k1+k2))
xlabel('Time/s')
ylabel(sprintf('IMF%d', k1+k2));
subplot(4,2,2*k2)
= FFTAnalysis(imf{k1+k2}, Ts);
plot(f, yf)
title(sprintf('第%d个IMF的频谱', k1+k2))
xlabel('f/Hz')
ylabel('|IMF(f)|');
end
end
figure
subplot(211)
plot(c,x)
set(gca,'FontSize',8,'XLim',);
title('原始信号')
xlabel('Time/s')
ylabel('Origin');
subplot(212)
= FFTAnalysis(x, Ts);
plot(f, Yf)
title('原始信号的频谱')
xlabel('f/Hz')
ylabel('|Y(f)|');
function plot_hht_3d(imf,numfreq,fs,ANGLE) if nargin<3 fs=1; ANGLE=[-37.5,30]; end if nargin<4 if size(fs,2)>1 ANGLE=fs; fs=1; else ANGLE=[-37.5,30]; end end N=size(imf,2); =hhspectrum(imf); =size(f); MaxFreq=max(max(f)); MaxFreq=ceil(MaxFreq/0.5)*0.5; if nargin<2 numfreq=512; end df=linspace(0,MaxFreq,numfreq); Spectrum=zeros(numfreq,n); Temp=f; Temp=min(round((Temp/MaxFreq)*numfreq)+1,numfreq); for k=1:m for u=1:n Spectrum(Temp(k,u),u)=Spectrum(Temp(k,u),u)+A(k,u); end end df=df*fs; figure clf mesh(tt,df,Spectrum); set(gca,'XLim',); xlabel('采样点数/n'); if fs==1 ylabel('归一化频率'); else ylabel('频率/Hz'); end zlabel('幅值'); title('三维联合时频图'); colormap jet; shading interp; view(ANGLE(1),ANGLE(2)); set(gca,'YLim',); end
页:
[1]