emd分解 plot_hht 怎么出错了？

yiyezhiqiuwq

源程序
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);
一运行到
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: [x,Fs] = 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
[u,v] = 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([0 c(end)]);
    ylim([0 1/2/Ts]);
    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',[0 c(end)]);
        title(sprintf('第%d个IMF', k1+k2))
        xlabel('Time/s')
        ylabel(sprintf('IMF%d', k1+k2));
      
        subplot(4,2,2*k2)
        [yf, f] = 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',[0 c(end)]);
title('原始信号')
xlabel('Time/s')
ylabel('Origin');
subplot(212)
[Yf, f] = 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); [A,f,tt]=hhspectrum(imf); [m,n]=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',[0,N]); xlabel('采样点数/n'); if fs==1     ylabel('归一化频率'); else     ylabel('频率/Hz'); end zlabel('幅值'); title('三维联合时频图'); colormap jet; shading interp; view(ANGLE(1),ANGLE(2)); set(gca,'YLim',[0,fs/2]); end