%------------------小波分解重构算法---------------
%------------------------------------------------
clear;clc;
% ----------选择DB4小波，确定系数（参考《实用小波分析》）
c0=0.4829629131445341;
c1=0.8365163037378079;
c2=0.2241438680420134;
c3=-0.1294095225512604;
%--------------产生模拟的振动信号，包含瞬态信号
p=0.0002;%采样频率为5000
t2=2;t1=0;
t=t1:p:t2;
N=(t2-t1)/p;
x1=sin(2*pi*65*t)+2*sin(2*pi*97*t);%振动信号
tt1=1.2;tt2=1.4;
tt=tt1:p:tt2;
x2=0.1*sin(2*pi*500*tt+0.98);
x2=[zeros(1,tt1/p+1),x2,zeros(1,(t2-tt2)/p)];%瞬态信号
s=x1+x2;
figure(1)
subplot(311);plot(s);title('原始振动信号');
axis([0,N,-5,5])
grid on
%---------------------------------------------
% 一尺度小波正变换
j=1;
n=size(s,2);
for i=1:2:n-3
tras(j)=c0*s(i)+c1*s(i+1)+c2*s(i+2)+c3*s(i+3);
tras(j+1)=c3*s(i)-c2*s(i+1)+c1*s(i+2)-c0*s(i+3);
j=j+2;
end
tras(j)=c0*s(n-1)+c1*s(n)+c2*s(1)+c3*s(2);
tras(j+1)=c3*s(n-1)-c2*s(n)+c1*s(1)-c0*s(2);
% 排序，系数A1在前D1在后
j=0;
for i=1:2:n-1
j=j+1;
sort(j)=tras(i);
end
for i=2:2:n-1
j=j+1;
sort(j)=tras(i);
end
subplot(323);plot(sort(1:n/2));title('近似系数A1');grid on
subplot(324);plot(sort(n/2+1:n-1));title('细节系数D1');grid on
% 对一尺度上的近似系数A1再进行小波分解
tras=zeros;
s=sort(1:n/2);
j=1;
n=size(s,2);
for i=1:2:n-3
tras(j)=c0*s(i)+c1*s(i+1)+c2*s(i+2)+c3*s(i+3);
tras(j+1)=c3*s(i)-c2*s(i+1)+c1*s(i+2)-c0*s(i+3);
j=j+2;
end
tras(j)=c0*s(n-1)+c1*s(n)+c2*s(1)+c3*s(2);
tras(j+1)=c3*s(n-1)-c2*s(n)+c1*s(1)-c0*s(2);
j=0;
for i=1:2:n
j=j+1;
sort(j)=tras(i);
end
for i=2:2:n
j=j+1;
sort(j)=tras(i);
end
subplot(325);plot(sort(1:n/2));title('近似系数cA2');grid on
subplot(326);plot(sort(n/2+1:n));title('细节系数cD2');grid on
% 对原信号进行重构
% 逆排序
figure(2)
subplot(321);plot(sort(1:n/2));title('近似系数cA2');grid on
subplot(322);plot(sort(n/2+1:n));title('细节系数cD2');grid on
j=-1;
for i=1:1:n/2
j=j+2;
sort_r(j)=sort(i);
end
j=0;
for i=n/2+1:1:n
j=j+2;
sort_r(j)=sort(i);
end
% 小波逆变换
tras=zeros;
tras(1)=c2*sort_r(n-1)+c1*sort_r(n)+c0*sort_r(1)+c3*sort_r(1);
tras(2)=c3*sort_r(n-1)-c0*sort_r(n)+c1*sort_r(1)-c2*sort_r(1);
j=1;
for i=1:2:n-3
j=j+2;
tras(j)=c2*sort_r(i)+c1*sort_r(i+1)+c0*sort_r(i+2)+c3*sort_r(i+3);
tras(j+1)=c3*sort_r(i)-c0*sort_r(i+1)+c1*sort_r(i+2)-c2*sort_r(i+3);
end
sort(1:n)=tras(1:n);
% 逆排序
n=size(sort,2);
subplot(323);plot(sort(1:n/2));title('近似系数cA1');grid on
subplot(324);plot(sort(n/2+1:n));title('细节系数cD1');grid on
j=-1;
for i=1:1:n/2
j=j+2;
sort_r(j)=sort(i);
end
j=0;
for i=n/2+1:1:n
j=j+2;
sort_r(j)=sort(i);
end
% 第二次小波逆变换
tras=zeros;
tras(1)=c2*sort_r(n-1)+c1*sort_r(n)+c0*sort_r(1)+c3*sort_r(1);
tras(2)=c3*sort_r(n-1)-c0*sort_r(n)+c1*sort_r(1)-c2*sort_r(1);
j=1;
for i=1:2:n-3
j=j+2;
tras(j)=c2*sort_r(i)+c1*sort_r(i+1)+c0*sort_r(i+2)+c3*sort_r(i+3);
tras(j+1)=c3*sort_r(i)-c0*sort_r(i+1)+c1*sort_r(i+2)-c2*sort_r(i+3);
end
subplot(3,2,[5,6]);plot(tras);title('恢复的原始振动信号');grid on
评论0