zhangxiskd 发表于 2013-1-23 00:05

hht工具箱已经安装上了为什么还出错??求大侠解答

??? Error using ==> instfreq at 66
T can not be equal to 1 neither to the last element of X
Error in ==> hhspectrum at 45
f(i,:)=instfreq(an(i,:)',tt,l)';
Error in ==> bianjipu at 25
=hhspectrum(imf(1:end-1, :));%对IMF分量求取瞬时频率与振幅:A:是每个IMF的振幅向量,f:每个IMF对应的瞬时频率,t:时间序列号

源代码
clear all;
clc;
%给定参数
fs=1000;
N=1000;
n=0:N-1;
t=n/fs;
%构造仿真信号
x=2*sin(2*pi*120*t)+(1+sin(2*pi*5*t)).*sin(2*pi*50*t+0.5*cos(2*pi*5*t));
%emd分解
imf=emd(x);
%绘制时域波形
figure,plot(t,x);xlabel('时间 t/s');ylabel('幅值');title('仿真信号时域波形');
%绘制FFT频谱
nfft= 2^nextpow2(length(t));%找出大于y的个数的最大的2的指数值
ff=fs*(0:nfft/2-1)/nfft;%�T变换后对应的频率的序列
fftx=fft(x,nfft);%求FFT变换
ps=fftx.*conj(fftx)/nfft;%conj()函数是求y函数的共轭复数,实数的共轭复数是他本身
figure;
subplot(211),plot(ff,abs(fftx(1:nfft/2))*2/nfft);ylabel('幅值');xlabel('频率');title('FFT频谱');
subplot(212),plot(ff,ps(1:nfft/2));ylabel('功率谱密度');xlabel('频率');title('信号功率谱');
%绘制emd分解结果
%plot_imf(x,t,imf);
%求时频谱
=hhspectrum(imf(1:end-1, :));%对IMF分量求取瞬时频率与振幅:A:是每个IMF的振幅向量,f:每个IMF对应的瞬时频率,t:时间序列号
%绘制瞬时包络图和瞬时频率图
figure;
subplot(221),plot(t/N,A(1,:));xlabel('时间 t/s');ylabel('幅值');title('imf1分量瞬时包络');
subplot(222),plot(t/N,f(1,:)*fs);xlabel('时间 t/s');ylabel('频率');title('imf1分量瞬时频率');
subplot(223),plot(t/N,A(2,:));xlabel('时间 t/s');ylabel('幅值');title('imf2分量瞬时包络');
subplot(224),plot(t/N,f(2,:)*fs);xlabel('时间 t/s');ylabel('频率');title('imf2分量瞬时频率');
%即时频图(用颜色表示第三维值的大小)和三维图(三维坐标系:时间,中心频率,振幅)
=toimage(A,f,t,length(t));%将每个IMF信号合成求取Hilbert谱,E:对应的振幅值,Cenf:每个网格对应的中心频率这里横轴为时间,纵轴为频率      
%绘制Hilbert-Huang谱
figure;
set(gcf,'Color','w');
imagesc(t/N,,E);
set(gca,'YDir','normal')
colormap('jet')
colorbar;
xlabel('时间 t/s');
ylabel('频率 f/Hz');axis()
title('Hilbert-Huang Spectrum') ;
%画出边际谱
%N=length(Cenf);%设置频率点数   %完全从理论公式出发。网格化后中心频率很重要,大家从连续数据变为离散的角度去思考,相信应该很容易理解
for k=1:size(E,1)
    bjp(k)=sum(E(k,:))*1/fs;
end
figure;
plot(Cenf(1,:)*fs,bjp);% 作边际谱图   进行求取Hilbert谱时频率已经被抽样成具有一定窗长的离散频率,所以此时的频率轴已经是中心频率
xlabel('频率 f/Hz');
ylabel('幅值');
title('边际谱');

页: [1]
查看完整版本: hht工具箱已经安装上了为什么还出错??求大侠解答