边际谱和HHT谱的Matlab例子

xray

看到版上总有人在问边际谱和HHT谱的画法，又搜索了一下，好像没有这方面的主题帖子，就发两个以前写的小程序，权作抛砖引玉吧。

% 边际谱与FFT比较
clear

T = 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;  
end
f = (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谱和边际谱
clear

rand('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;  
end
f = (0:N-3)/N*(fs/2);
figure(2);
plot(f, bjp);

yanling66520

谢谢你，我昨天晚上研究了下，知道怎么用G. Rilling的程序了，原来只需要自己编程，设置一些参数，直接调用相关的函数就可以了，之前比较傻，云里雾里的搞了半天呵呵

yanling66520

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

[ 本帖最后由 yanling66520 于 2009-6-29 22:26 编辑 ]

hitdely

楼主你好 请问一下 我用了你的程序可是 报错了 我对MATLAB不熟 请帮我看看

jiangxiaolong

搜索时频工具箱，安装工具包就会有这个函数。

jiangxiaolong

楼主你好 请问一下 我用了你的程序可是还是有问题。我看了好多东西，不知道怎么解决。麻烦知道的人给小弟一提点。
??? Error: File: toimage.m Line: 59 Column: 1
This statement is not inside any function.
(It follows the END that terminates the definition of the function "toimage".)

ChaChing

原帖由 jiangxiaolong 于 2010-5-21 11:03 发表
...This statement is not inside any function.
(It follows the END that terminates the definition of the function "toimage".)

检查下end的对应情况看看!

jiangxiaolong

我弄出来了，非常感谢你的帮助。

ChaChing

建议LS说清楚些, 与大家分享你的成果!

jiangxiaolong

弄错了，我只是知道我的错误了，楼主的问题，我不会。不好意思。

501987105

请问能画出三维图吗

501987105

谢谢了
那位大哥帮帮小弟

xiaohouzi77

回复 ChaChing 的帖子

感谢，给解决大问题了   

shinhwaforever

谢谢各位了,看了这些真是受益匪浅.

a240612525

回复 xray 的帖子

请问楼主，第一个程序中边际谱和FFT谱的结果为什么会不一致？

边际谱

FFT谱