emd_多算例_三次样条镜像延拓_时频图三维图_希尔伯特谱边际谱_IMF分量及残余信号.docVIP

主程序：main clear; clc; %示例程序 N=1000; n=1:N; fs=1000; t=n/fs; %x=cos(2*pi*11*t); y=cos(2*pi*50*t)+1*cos(2*pi*3*t); z=x+y; %x=1+0.3*sin(pi*1/200*t).*sin(pi*1/3*t+3*sin(pi*1/60*t)); y=3*cos(5*pi*1/100*t); z=x+y; z = sin(2*pi*5*t) + sin(2*pi*14*t) + sin(2*pi*50*t) ; %+ 0.1*randn(1, length(t)); data=z; m=1; % number of intrinsic mode functions imf=emd_lmd1(data); [A,f,t]=hhspectrum(imf); %对IMF分量求取瞬时频率与振幅：A：是每个IMF的振幅向量,f:每个IMF对应的瞬时频率，t：时间序列号 [E,t,Cenf]=toimage(A,f); %将每个IMF信号合成求取Hilbert谱，E：对应的振幅值，Cenf：每个网格对应的中心频率 这里横轴为时间，纵轴为频率 %即时频图(用颜色表示第三维值的大小)和三维图（三维坐标系：时间，中心频率，振幅） cemd_visu(data,1:length(data),imf); %显示每个IMF分量及残余信号-------------------------------------------- disp_hhs(E); %希尔伯特谱---------------------------------------------------------- % colormap(flipud(gray)) % %画出边际谱 % %N=length(Cenf);%设置频率点数 %完全从理论公式出发。网格化后中心频率很重要，大家从连续数据变为离散的角度去思考，相信应该很容易理解 % for k=1:size(E,1) % bjp(k)=sum(E(k,:))*1/fs; % end % figure(3); % plot(Cenf(1,:)*fs,bjp); % 作边际谱图 进行求取Hilbert谱时频率已经被抽样成具有一定窗长的离散频率，所以此时的频率轴已经是中心频率 % xlabel(频率 / Hz); % ylabel(幅值); %%%%-------- 三维能量图测试 ---------------- L=size(imf,1); % Number of components in the decomposition % loop for on each component S=[]; clear x z m p freq x=imf; % now each column is a compenent z=hilbert(x); % analytic signal m=abs(z); % module of z p=angle(z); % phase of z for i=1:L freq(:,i)=instfreq(z(:,i)); % instantaneous frequency ceilfreq(:,i)=ceil(freq(:,i)*N); % to have integer values-also do a smoothing given the number for j=1:length(x)-2 S(ceilfreq(j,i),j+1)=m(j+1,i); end end eps=0.00001; % to avoid zero values before the log E=20*log(S.^2+eps); % Hilbert energ spectrum % plot S figure(90); t=1:length(x); % time sample f=t/length(x)*25; % frequency %%%%%%%%%5 imagesc(t,f,E); contour(t(1:999),1:500,E); %画出等高线图 surf(1:999,1:500,S.^2); shading interp; colormap(gray);axis