一段hht程序，为什么频率图出不来，求救~~
我输入的原始信号是t=0:1/(128*6400):0.35;
x=2*sin(50*pi*t)+sin(100*pi*t)+sin(170*pi*t)+sin(400*pi*t);
其中emd分解的程序为
function imf = emd(x)
x = transpose(x(:));
imf = [];
while ~ismonotonic(x)
x1 = x;
sd = Inf;
while (sd > 0.1) | ~isimf(x1)
s1 = getspline(x1);
s2 = -getspline(-x1);
x2 = x1-(s1+s2)/2;
sd = sum((x1-x2).^2)/sum(x1.^2);
x1 = x2;
end
imf{end+1} = x1;
x = x-x1;
end
imf{end+1} = x;
% FUNCTIONS
function u = ismonotonic(x)
u1 = length(findpeaks(x))*length(findpeaks(-x));
if u1 > 0, u = 0;
else, u = 1; end
function u = isimf(x)
N= length(x);
u1 = sum(x(1:N-1).*x(2:N) < 0);
u2 = length(findpeaks(x))+length(findpeaks(-x));
if abs(u1-u2) > 1, u = 0;
else, u = 1; end
function s = getspline(x)
N = length(x);
p = findpeaks(x);
s = spline(,,1:N);
最后变换和画图的程序为
function plot_hht(x,Ts)
% Plot the HHT.
% plot_hht(x,Ts)
%
% :: Syntax
% The array x is the input signal and Ts is the sampling period.
% Example on use: = wavread('Hum.wav');
% plot_hht(x(1:6000),1/Fs);
% Func : emd
% Get HHT.
imf = emd(x);
for k = 1:length(imf)
b(k) = sum(imf{k}.*imf{k});
th = angle(hilbert(imf{k}));
d{k} = diff(th)/(1/819200)/(2*pi);
end
= sort(-b);
b = 1-b/max(b);
% Set time-frequency plots.
N = length(x);
c = linspace(0,(N-2)*(1/819200),N-1);
for k = v(1:2)
figure, plot(c,d{k},'k.','Color',b(),'MarkerSize',3);
set(gca,'FontSize',8,'XLim',,'YLim',);
xlabel('Time'), ylabel('Frequency');
end
% Set IMF plots.
M = length(imf);
N = length(x);
c = linspace(0,(N-1)*Ts,N);
for k1 = 0:4:M-1
figure
for k2 = 1:min(4,M-k1), subplot(4,1,k2), plot(c,imf{k1+k2}); set(gca,'FontSize',8,'XLim',); end
xlabel('Time');
end
我用matlab仿真出来，后面能分解出8个IMf，但是Frequency的那两个图没有波形，求哪位大神帮忙看看呗~~~ 这些程序有什么用，实际中根本用不上，如果是做软件开发的可以，不是软件开发学这些真的多少用，实际中只要你能通过编程找出问题原因就OK。个人意见。
页:
[1]