不知楼主解决问题了没,我这右端程序,也不是我自己写的,别人的,你可以参考下
%多脊提取
clear
SampFreq = 10;
x=1/SampFreq:1/SampFreq:20;%原理想信号
%sig=exp(-0.045*x).*cos(2*pi*1.5*x+pi)+exp(-0.098*x).*cos(2*pi*1.0*x+pi/2);
sig=exp(-0.045*x).*cos(2*pi*2*x+pi)+exp(-0.098*x).*cos(3*pi*1.0*x+pi/2);
% sig=sig';
% noise=randn(1,200);
% noise=noise';
% dB=20;
% sig=sigmerge(sig,noise,dB);
% sig=awgn(sig,20)
% sig=sig';
fmax = 0.5; % 最高分析频率(归一化频率)
fmin = 0.005; % 最低分析频率(归一化频率)
fb = 4 ; % 取cmor4-2小波进行实验,带宽参数为4
fc = 2; % 中心频率2Hz
totalscal =512 ; % 所取尺度的数目
FreqBins = linspace(fmin,fmax,totalscal); % 将频率轴在分析范围内等间隔划分
Scales = fc./ FreqBins; % 计算相应的尺度参数
RealFreqBins = FreqBins * SampFreq; % 尺度所对应的实际频率
swd=cwt(sig,Scales,'cmor4-2');
swd=abs(swd);
figure(1)
imagesc(x,Scales,swd);
[a1,b1]=max(swd(2:120,:));
[a2,b2]=max(swd(100:250,:));
figure(2)
mesh(x,RealFreqBins,swd);
colormap jet;
shading interp;
colorbar;
axis([min(x) max(x) min(RealFreqBins) max(RealFreqBins)]);
ylabel('Frequency / Hz');
xlabel('Time / sec');
figure(3)
m=0.495/512*10*(b1+2);
plot(m)
axis([1 200 0.05 5])
hold on
n=0.495/512*10*(b2+100);
plot(n);
axis([1 200 0.05 5]) |