本帖最后由 1393107100 于 2019-5-1 11:10 编辑

clear,clc

close all

%%%%%%　小波变换　%%%%%%%%%%%%%%%%

fs=1024;

t=1/fs:1/fs:1;

f1=100;f2=200;f3=300;

s=cos(2*pi*f1*t.*(t>=0&t<0.3))+2*cos(2*pi*f2*t.*(t>=0.3&t<0.8))+3*sin(2*pi*f3*t.*(t>=0.8&t<=1));

%s=cos(2*pi*f1*t);

figure(1);

plot(t,s,'b')

title('原始单频信号')

Y = awgn(s,6);%加高斯白噪声

subplot(412);

plot(t,Y,'r')

title('加信噪比为6的白噪声')

%%%%%%%%%%%%%%%%%小波时频图绘制 %%%%%%%%%%%%%%%%%%%%%

wavename='cmor4-4';    %%选用带宽参数和中心频率均为4 的复morlet小波

totalscal=256;          %尺度序列的长度，即scal的长度

fc=centfrq(wavename);%% 小波的中心频率

%%%%%%%%% 绘制时间域小波波形 %%%%%%%%%%%%%%%%%%%

cparam=2*fc*totalscal;%为得到合适的尺度所求出的参数

a=totalscal:-1:1;

scal=cparam./a;%得到各个尺度，以使转换得到频率序列为等差序列

coefs=cwt(Y,scal,wavename,1/fs);%求连续小波系数,小波系数coefs(系数是复数时要取模)

f=scal2frq(scal,wavename,1/fs);% 将尺度转换为频率.F = scal2frq(A,'wname',DELTA),该函数能将尺度转换为实际频率，其中A为尺度，wname为小波名称，DELTA为采样周期。

figure(2);

imagesc(t,f,abs(coefs)); %绘制色谱图

colormap(jet)

colorbar

xlabel('时间 t/s')

ylabel('频率 f/Hz')

title('小波时频图')

abs_coefs=abs(coefs).^2;%abs函数：数值的绝对值和复数的幅值,

bw=abs_coefs>=0.7*mean(mean(abs_coefs));%mean数组的均值.mean(A)如果A为矩阵，那么返回包含每列均值的行向量。

bw=bwlabel(bw);%在BW数组中，0代表黑背景，1代表白

for i=1:max(max(bw))

if length(find(bw==i))>1000

bw(find(bw==i))=-1;

end

end

bw=bwlabel(bw==-1);

figure(3)

imagesc (bw); title('脊线区域')

out=zeros(size(coefs));% 输出：大小为输入的矩阵大小，各元素的值为0

for i=1:max(max(bw))

out=out+local_maxima(abs_coefs,bw,i);%小波系数模局部极大值

end

out=out.*bw;

RIDGE=zeros(size(coefs));

for i=1:max(max(out))

[x,y]=find(out==i);

x=x(1:100:end);

if length(x)<3

break;

end

cs = spline(y(1:100:end),x);

%spline三次方样条数据插值

%pp = spline(x,y) 返回一个分段多项式结构体以用于 ppval 和样条实用工具 unmkpp

x_new=round(ppval(cs,y(1):y(end)));

x_new(find(x_new<=0))=1;

for i=1:length(x_new)

RIDGE(x_new(i),y(i))=1;

end

end

figure,imagesc (RIDGE);  title('脊线')

时频图.jpg

(22.45 KB, 下载次数: 2)

2019-5-1 11:10 上传

时频图

脊线.jpg

(14.9 KB, 下载次数: 0)

2019-5-1 11:08 上传

脊线图