一维小波分解与重构_多尺度一维小波分解

于 2021-01-08 16:35:01 发布
阅读量744 收藏 12
点赞数 3
文章标签： 一维小波分解与重构
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/weixin_30014265/article/details/113451183
版权
本文通过code1.m和DWT_new.m, DWT.m程序探讨了一维小波分解的方法，并讲解了如何计算EMD各IMF分量的能量。" 90695718,7957678,使用Flask与CKEditor5实现富文本编辑及图片上传,"['Flask框架', '前端开发', '后端开发']
摘要由CSDN通过智能技术生成

6855d0667eb8b805c19f6fbea6eccf19.png

code1.m

fc=100;%信号频率%调频信号频率
fz=5;%采样频率
Fs=1000;%采样频率
N=1024;%采样频率
n=0:N-1;%时间轴离散,N个点
t=(0:N-1)/Fs;                      %时间
s=5*sin(2*pi*150*t)+ 5*sin(2*pi*50*t);%产生频率调制信号
%db10小波进行4层分解
[c,l]=wavedec(s, 4, 'db10');%多尺度一维小波分解函数
%重构第1~4层细节信号
d4=wrcoef('d',c,l,'db10',4);
d3=wrcoef('d',c, l,'db10',3);
d2=wrcoef('d',c, l,'db10', 2);
d1=wrcoef('d',c,l,'db10',1);
%显示细节重构信号及原始信号(略),如图5.19(a)所示
wavename='cmor3-3';%复 Morlet小波
totalscal=256;%尺度
Fc=centfrq(wavename);%小波的中心频率
c=2*Fc*totalscal;
scals=c./(1:totalscal);%得到各个尺度,以使转换得到频率序列为等差序列
f=scal2frq (scals, wavename, 1/Fs);%将尺度转换为频率
coefs=cwt(s, scals, wavename);%求连续小波系数
figure
imagesc(t, f, abs(coefs))%绘制色谱图,如图5.19(a)所示colorbar%绘制伪彩色条
colorbar

6aeb8be73c4434e3907247481fb5df7f.png

小波包分解的程序 DWT_new.m

clc,clear
a=textread('1500.txt','%.10f','headerlines',2); 
n=length(a)/3 ;
for i=1:1:n
    b0(i)=a(i)/2^14*4;%转矩
    b1(i)=a(i+n)/2^14*4;%电流 
   %b1(i)=a(i+n)/1000;%m法无滤波 
    b2(i)=a(i+2*n)/1000;%转速
end
for i=0:(length(b2)-1)/4
    x_input(i+1)=b2(4*i+1);
end
fs=500;%采样频率
t=2:1/fs:(4-1/fs);
N=length(t);
x=x_input(4001:5000);
x=x-mean(x);
y=fft(x,N);
y1=abs(y);
ayy=y1/(N/2);
ayy(1)=ayy(1)/2;
f=([1:N]-1)*fs/N;
plot(f(1:N/2),ayy(1:N/2));
hold on
xlim([0 500])
%ylim([0 0.02])
set(gca,'fontsize',16);
grid on
xlabel('频率(Hz)','fontsize',16);
ylabel('幅值(r/min)','fontsize',16);
[wpt]=wpdec(x,3,'dmey');  %小波包分解，3代表分解3层，'dmey'使用meyr小波
plot(wpt)               %画小波包树图
wpviewcf(wpt,1);        %画出时间频率图
nodes=[7;8;9;10;11;12;13;14];   %第3层的节点号
ord=wpfrqord(nodes); %小波包系数重排，ord是重排后小波包系数索引构成的矩阵　如3层分解的[1;2;4;3;7;8;6;5]
nodes_ord=nodes(ord); %重排后的小波系数
for i=1:8
    rex3(:,i)=wprcoef(wpt,nodes_ord(i));  %实现对节点小波节点进行重构
end
for i=0:2
    figure(i+4);  %绘制第3层各个节点分别重构后信号的频谱
    x_sign=rex3(:,i+1);
    N=length(x_sign); %采样点个数
    signalFFT=abs(fft(x_sign,N));%真实的幅值
    Y=2*signalFFT/N;
    f=(0:N/2)*(fs/N);
    plot(f,Y(1:N/2+1));
    ylabel('amp'); xlabel('frequency');
    grid on
    axis([0 250 0 2]);
    title(['小波包第3层',num2str(i),'节点信号频谱']);
end
%%
%wavelet packet coefficients. 求取小波包分解的各个节点的小波包系数
cfs3_0=wpcoef(wpt,nodes_ord(1));  %对重排序后第3层0节点的小波包系数0-31.25Hz
cfs3_1=wpcoef(wpt,nodes_ord(2));  %对重排序后第3层0节点的小波包系数31.25-62.5Hz
cfs3_2=wpcoef(wpt,nodes_ord(3));  %对重排序后第3层0节点的小波包系数62.5-93.75Hz
cfs3_3=wpcoef(wpt,nodes_ord(4));  %对重排序后第3层0节点的小波包系数93.75-125Hz
cfs3_4=wpcoef(wpt,nodes_ord(5));  %对重排序后第3层0节点的小波包系数125-156.25Hz
cfs3_5=wpcoef(wpt,nodes_ord(6));  %对重排序后第3层0节点的小波包系数156.25-187.5Hz
cfs3_6=wpcoef(wpt,nodes_ord(7));  %对重排序后第3层0节点的小波包系数187.5-218.75Hz
cfs3_7=wpcoef(wpt,nodes_ord(8));  %对重排序后第3层0节点的小波包系数218.75-250Hz

E_cfs3_0=norm(cfs3_0,2)^2;  %% 1-范数：就是norm(...,1)，即各元素绝对值之和；2-范数：就是norm(...,2)，即各元素平方和开根号；
E_cfs3_1=norm(cfs3_1,2)^2;
E_cfs3_2=norm(cfs3_2,2)^2;
E_cfs3_3=norm(cfs3_3,2)^2;
E_cfs3_4=norm(cfs3_4,2)^2;
E_cfs3_5=norm(cfs3_5,2)^2;
E_cfs3_6=norm(cfs3_6,2)^2;
E_cfs3_7=norm(cfs3_7,2)^2;
E_total=E_cfs3_0+E_cfs3_1+E_cfs3_2+E_cfs3_3+E_cfs3_4+E_cfs3_5+E_cfs3_6+E_cfs3_7;
p_node(1)= 100*E_cfs3_0/E_total;           % 求得每个节点的占比
p_node(2)= 100*E_cfs3_1/E_total;           % 求得每个节点的占比
p_node(3)= 100*E_cfs3_2/E_total;           % 求得每个节点的占比
p_node(4)= 100*E_cfs3_3/E_total;           % 求得每个节点的占比
p_node(5)= 100*E_cfs3_4/E_total;           % 求得每个节点的占比
p_node(6)= 100*E_cfs3_5/E_total;           % 求得每个节点的占比
p_node(7)= 100*E_cfs3_6/E_total;           % 求得每个节点的占比
p_node(8)= 100*E_cfs3_7/E_total;           % 求得每个节点的占比

figure;
x=1:8;
bar(x,p_node);
title('各个频段能量所占的比例');
xlabel('频率 Hz');
ylabel('能量百分比/%');
for j=1:8
    text(x(j),p_node(j),num2str(p_node(j),'%0.2f'),...
        'HorizontalAlignment','center',...
        'VerticalAlignment','bottom')
end
%%
%绘制重构前各个特征频段小波包系数的图形
figure(1);
subplot(4,1,1);
plot(x_input);
title('原始信号');
subplot(4,1,2);
plot(cfs3_0);
title(['层数 ',num2str(3) '  节点 0的小波0-31.25Hz',' 系数'])
subplot(4,1,3);
plot(cfs3_1);
title(['层数 ',num2str(3) '  节点 1的小波31.25-62.5Hz',' 系数'])
subplot(4,1,4);
plot(cfs3_2);
title(['层数 ',num2str(3) '  节点 2的小波62.5-125Hz',' 系数'])
%绘制重构后时域各个特征频段的图形
figure(3);
subplot(3,1,1);
plot(rex3(:,1));
title('重构后0-31.25Hz频段信号');
subplot(3,1,2);
plot(rex3(:,2));
title('重构后31.25-62.5Hz频段信号')
subplot(3,1,3);
plot(rex3(:,3));
title('重构后62.5-93.75Hz频段信号');
%%
[imf,residual,info] = emd(rex3(:,2),'Interpolation','pchip');
%[imf,residual,info] = emd(x,'Interpolation','pchip');
hht(imf,fs)

DWT.m

clc,clear
a=textread('2000.txt','%.10f','headerlines',2); 
n=length(a)/3 ;
for i=1:1:n
    b0(i)=a(i)/2^14*4;%转矩
    b1(i)=a(i+n)/2^14*4;%电流 
   %b1(i)=a(i+n)/1000;%m法无滤波 
    b2(i)=a(i+2*n)/1000;%转速
end
%for i=0:(length(b2)-1)/2
    %data(i+1)=b2(2*i+1);
%end
fs=2000;%采样频率
t=2:0.0005:(10-0.0005);
N=length(t);
x=b2(4001:20000);
x=x-mean(x);
x_input=x;
[wpt]=wpdec(x_input,3,'dmey');  %小波包分解，3代表分解3层，'dmey'使用meyr小波
plot(wpt)               %画小波包树图
wpviewcf(wpt,1);        %画出时间频率图
nodes=[7;8;9;10;11;12;13;14];   %第3层的节点号
ord=wpfrqord(nodes); %小波包系数重排，ord是重排后小波包系数索引构成的矩阵　如3层分解的[1;2;4;3;7;8;6;5]
nodes_ord=nodes(ord); %重排后的小波系数
for i=1:8
    rex3(:,i)=wprcoef(wpt,nodes_ord(i));  %实现对节点小波节点进行重构
end
%%
figure(3)
 x_sign=rex3(:,1);
    N=length(x_sign); %采样点个数
    signalFFT=abs(fft(x_sign,N));%真实的幅值
    Y=2*signalFFT/N;
    f=(0:N/2)*(fs/N);
    plot(f,Y(1:N/2+1));
    ylabel('amp'); xlabel('frequency');
    grid on
    xlim([0 200])
    %axis([0 128 0 2]);
    title(['小波包第3层',num2str(0),'节点信号频谱']);
%%
    for i=0:7
    figure(i+3);  %绘制第3层各个节点分别重构后信号的频谱
    x_sign=rex3(:,i+1);
    N=length(x_sign); %采样点个数
    signalFFT=abs(fft(x_sign,N));%真实的幅值
    Y=2*signalFFT/N;
    f=(0:N/2)*(fs/N);
    plot(f,Y(1:N/2+1));
    ylabel('amp'); xlabel('frequency');
    grid on
    %axis([0 128 0 0.5]);
    title(['小波包第3层',num2str(i),'节点信号频谱']);
end
%%
%wavelet packet coefficients. 求取小波包分解的各个节点的小波包系数
cfs3_0=wpcoef(wpt,nodes_ord(1));  %对重排序后第3层0节点的小波包系数0-8Hz
cfs3_1=wpcoef(wpt,nodes_ord(2));  %对重排序后第3层0节点的小波包系数8-16Hz
cfs3_2=wpcoef(wpt,nodes_ord(3));  %对重排序后第3层0节点的小波包系数16-24Hz
cfs3_3=wpcoef(wpt,nodes_ord(4));  %对重排序后第3层0节点的小波包系数24-32Hz
cfs3_4=wpcoef(wpt,nodes_ord(5));  %对重排序后第3层0节点的小波包系数32-40Hz
cfs3_5=wpcoef(wpt,nodes_ord(6));  %对重排序后第3层0节点的小波包系数40-48Hz
cfs3_6=wpcoef(wpt,nodes_ord(7));  %对重排序后第3层0节点的小波包系数48-56Hz
cfs3_7=wpcoef(wpt,nodes_ord(8));  %对重排序后第3层0节点的小波包系数56-64Hz

E_cfs3_0=norm(cfs3_0,2)^2;  %% 1-范数：就是norm(...,1)，即各元素绝对值之和；2-范数：就是norm(...,2)，即各元素平方和开根号；
E_cfs3_1=norm(cfs3_1,2)^2;
E_cfs3_2=norm(cfs3_2,2)^2;
E_cfs3_3=norm(cfs3_3,2)^2;
E_cfs3_4=norm(cfs3_4,2)^2;
E_cfs3_5=norm(cfs3_5,2)^2;
E_cfs3_6=norm(cfs3_6,2)^2;
E_cfs3_7=norm(cfs3_7,2)^2;
E_total=E_cfs3_0+E_cfs3_1+E_cfs3_2+E_cfs3_3+E_cfs3_4+E_cfs3_5+E_cfs3_6+E_cfs3_7;
p_node(1)= 100*E_cfs3_0/E_total;           % 求得每个节点的占比
p_node(2)= 100*E_cfs3_1/E_total;           % 求得每个节点的占比
p_node(3)= 100*E_cfs3_2/E_total;           % 求得每个节点的占比
p_node(4)= 100*E_cfs3_3/E_total;           % 求得每个节点的占比
p_node(5)= 100*E_cfs3_4/E_total;           % 求得每个节点的占比
p_node(6)= 100*E_cfs3_5/E_total;           % 求得每个节点的占比
p_node(7)= 100*E_cfs3_6/E_total;           % 求得每个节点的占比
p_node(8)= 100*E_cfs3_7/E_total;           % 求得每个节点的占比

figure;
x=1:8;
bar(x,p_node);
title('各个频段能量所占的比例');
xlabel('频率 Hz');
ylabel('能量百分比/%');
for j=1:8
    text(x(j),p_node(j),num2str(p_node(j),'%0.2f'),...
        'HorizontalAlignment','center',...
        'VerticalAlignment','bottom')
end
%% 绘制重构前各个特征频段小波包系数的图形
figure(1);
subplot(4,1,1);
plot(x_input);
title('原始信号');
subplot(4,1,2);
plot(cfs3_0);
title(['层数 ',num2str(3) '  节点 0的小波0-125Hz',' 系数'])
subplot(4,1,3);
plot(cfs3_1);
title(['层数 ',num2str(3) '  节点 1的小波125-250Hz',' 系数'])
subplot(4,1,4);
plot(cfs3_2);
title(['层数 ',num2str(3) '  节点 2的小波250-375Hz',' 系数'])
%绘制重构后时域各个特征频段的图形
figure(3);
subplot(3,1,1);
plot(rex3(:,1));
title('重构后0-125Hz频段信号');
subplot(3,1,2);
plot(rex3(:,2));
title('重构后125-250Hz频段信号')
subplot(3,1,3);
plot(rex3(:,3));
title('重构后250-375Hz频段信号');
%%
[imf,residual,info] = emd(rex3(:,1),'Interpolation','pchip');
%[imf,residual,info] = emd(x,'Interpolation','pchip');
hht(imf,fs)
%%
emd(x,'Interpolation','pchip')
%%
x=imf(:,1);
%x=x-mean(x);
y=fft(x,N);
y1=abs(y);
ayy=y1/(N/2);
ayy(1)=ayy(1)/2;
f=([1:N]-1)*fs/N;
% figure;
figure
plot(f(1:N/2),ayy(1:N/2));
hold on
xlim([0 100])
%ylim([0 0.005])
set(gca,'fontsize',10);
grid on
xlabel('频率(Hz)','fontsize',10);
ylabel('幅值(r/min)','fontsize',10);
%%
for i=0:7
subplot(2,4,i+1);
x_sign= rex3(:,i+1); 
N=length(x_sign); %采样点个数
signalFFT=abs(fft(x_sign,N));%真实的幅值
Y=2*signalFFT/N;
f=(0:N/2)*(fs/N);
plot(f,Y(1:N/2+1));
ylabel('amp'); xlabel('frequency');grid on
%axis([0 50 0 0.03]);
title(['小波包第3层',num2str(i),'节点信号频谱']);
end

如何求EMD各IMF分量的能量？

%经验模态分解
z=A1;
c=emd(z);%M行，N列
[m1,m2]=size(c);

%求各IMF的能量
for i=1:m1
    E(i)=sum(c(i,:).*conj(c(i,:)));
end

%求信号总能量
E_total1=sum(A1.*conj(A1)); %原始信号能量  E_total1=1.2891e+007
E_total2=sum(E); %各IMF能量总和 E_total2=1.5142e+007.
Java架构月亮
关注
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值