matlab绘制hht边际谱,7、边际谱和HHT谱的Matlab例子.doc

7、边际谱和HHT谱的Matlab例子.doc

7、边际谱和HHT谱的Matlab例子

看到版上总有人在问边际谱和HHT谱的画法，又搜索了一下，好像没有这方面的主题帖子，就发两个以前写的小程序，权作抛砖引玉吧。% 边际谱与FFT比较clearT = 1;? ?? ?? ?? ?? ?? ?? ?? ?% 仿真时间f1 = 15.2;f2 = 40;fs = 1000;? ?? ?? ?? ?? ?? ???% 采样率N = T*fs;n = 1:N;s = sin(2*pi*f1/fs*n) + sin(2*pi*f2/fs*n);s_fft = abs(fft(s))/N;imf = emd(s); [A, fa, tt] = hhspectrum(imf);[E, tt1] = toimage(A,fa,tt,length(tt));for k=1:size(E,1)? ? bjp(k) = sum(E(k,:))*1/fs*1/T;??endf = (0:N-3)/N*(fs/2);figure(1);plot(f,bjp);xlabel('频率 / Hz');ylabel('幅值');figure(2);plot(0:fs/N:fs/2-fs/N, s_fft(1:end/2))% 实际信号的HHT谱和边际谱clearrand('seed', 0);T = 0.01;? ?? ?? ?? ?? ?? ?? ?? ?% 仿真时间R = 5000;? ?? ?? ?? ?? ???% 码速率fd = 10000;? ?? ?? ?? ???% 载波频差fc = 20000;? ?? ?? ?? ???% 载波频率fs = 200000;? ?? ?? ?? ?% 采样率samp = fs/R;? ?? ?? ?? ?% 每个码元上的采样点数N = T*fs;n = 1:N;x = randint(1, R*T, 2);y = fskmod(x, 2, fd, samp, fs);y = y .* exp(i*2*pi*fc/fs*n);y = real(y);% z = awgn(y, 20, 'measured'); z = y;imf = emd(z); [A, fa, tt] = hhspectrum(imf);if size(imf,1) > 1? ? [A,fa,tt] = hhspectrum(imf(1:end-1, :));else? ? [A,fa,tt] = hhspectrum(imf);end[E, tt1] = toimage(A,fa,tt,length(tt));disp_hhs(E, tt1);% 使用灰度图显示% colormap(gray(255))? ?? ???for k = 1:size(E,1)? ? bjp(k) = sum(E(k,:))*1/fs*1/T;??endf = (0:N-3)/N*(fs/2);figure(2); plot(f, bjp);

答：我用g.rilling的程序对自己的信号分析了下，感觉做出的hht谱有些奇怪，好像不对，那跟我的原始信号很相似，上传麻烦您给看下。后面又用你给的程序做出了边际谱，也不知道对不对，也上传给你看下，感觉我的边际谱很乱。我的数据也在附件里附上了。另外最后一张图，是我之前用plot_hht的程序做出的时频图也上传过来看下。

hht谱图

边际谱

plot_hht程序做的时频图

8、关于hilbert谱图的问题

请教：为什么我画的Hilbert谱的二维图颜色条显示都是负数呢？之前用另外一个程序画的图也是这样的.程序在附件中给出

答：原因很簡單啊…help disp_hhs就知道答案了…因為圖中的數值皆除以最大值，而後取dB值數值變小，當然dB值會是負值。