% mallet_wavelet.m
% 此函数用于研究Mallet算法及滤波器设计
% 此函数用于消噪处理
%角度赋值
%此处赋值使滤波器系数恰为db9
%分解的高频系数采用db9较好,即它的消失矩较大
%分解的有用信号小波高频系数基本趋于零
%对于噪声信号高频分解系数很大,便于阈值消噪处理
clc;
clear;
close all;
[l,h]=wfilters('db10','d');
low_construct=l;
L_fre=20; %滤波器长度
low_decompose=low_construct(end:-1:1); %确定h0(-n),低通分解滤波器
for i_high=1:L_fre; %确定h1(n)=(-1)^n,高通重建滤波器
if(mod(i_high,2)==0);
coefficient=-1;
else
coefficient=1;
end
high_construct(1,i_high)=low_decompose(1,i_high)*coefficient;
end
high_decompose=high_construct(end:-1:1); %高通分解滤波器h1(-n)
L_signal=100; %信号长度
n=1:L_signal; %原始信号赋值
f=10;
t=0.001;
y=10*cos(2*pi*50*n*t).*exp(-30*n*t);
zero1=zeros(1,60); %信号加噪声信号产生
zero2=zeros(1,30);
noise=[zero1,3*(randn(1,10)-0.5),zero2];
y_noise=y+noise;
figure(1);
subplot(2,1,1);
plot(y);
title('原信号');
subplot(2,1,2);
plot(y_noise);
title('受噪声污染的信号');
check1=sum(high_decompose); %h0(n),性质校验
check2=sum(low_decompose);
check3=norm(high_decompose);
check4=norm(low_decompose);
l_fre=conv(y_noise,low_decompose); %卷积
l_fre_down=dyaddown(l_fre); %抽取,得低频细节
h_fre=conv(y_noise,high_decompose);
h_fre_down=dyaddown(h_fre); %信号高频细节
figure(2);
subplot(2,1,1)
plot(l_fre_down);
title('小波分解的低频系数');
subplot(2,1,2);
plot(h_fre_down);
title('小波分解的高频系数');
% 消噪处理
for i_decrease=31:44;
if abs(h_fre_down(1,i_decrease))>=0.000001
h_fre_down(1,i_decrease)=(10^-7);
end
end
l_fre_pull=dyadup(l_fre_down); %0差值
h_fre_pull=dyadup(h_fre_down);
l_fre_denoise=conv(low_construct,l_fre_pull);
h_fre_denoise=conv(high_construct,h_fre_pull);
l_fre_keep=wkeep(l_fre_denoise,L_signal); %取结果的中心部分,消除卷积影响
h_fre_keep=wkeep(h_fre_denoise,L_signal);
sig_denoise=l_fre_keep+h_fre_keep; %消噪后信号重构
%平滑处理
for j=1:2
for i=60:70;
sig_denoise(i)=sig_denoise(i-2)+sig_denoise(i+2)/2;
end;
end;
compare=sig_denoise-y; %与原信号比较
figure(3);
subplot(3,1,1)
plot(y);
ylabel('y'); %原信号
subplot(3,1,2);
plot(sig_denoise);
ylabel('sig\_denoise'); %消噪后信号
subplot(3,1,3);
plot(compare);
ylabel('compare'); %原信号与消噪后信号的比较
% clc;clear;
% %% 1.正弦波定义
% f1=50; % 频率1
% f2=100; % 频率2
% fs=2*(f1+f2); % 采样频率
% Ts=1/fs; % 采样间隔
% N=120; % 采样点数
% n=1:N;
% y=sin(2*pi*f1*n*Ts)+sin(2*pi*f2*n*Ts); % 正弦波混合
% figure(1)
% plot(y);
% title('两个正弦信号')
% figure(2)
% stem(abs(fft(y)));
% title('两信号频谱')
% %% 2.小波滤波器谱分析
% h=wfilters('db30','l'); % 低通
% g=wfilters('db30','h'); % 高通
% h=[h,zeros(1,N-length(h))]; % 补零(圆周卷积,且增大分辨率变于观察)
% g=[g,zeros(1,N-length(g))]; % 补零(圆周卷积,且增大分辨率变于观察)
% figure(3);
% stem(abs(fft(h)));
% title('低通滤波器图')
% figure(4);
% stem(abs(fft(g)));
% title('高通滤波器图')
% %% 3.MALLET分解算法(圆周卷积的快速傅里叶变换实现)
% sig1=ifft(fft(y).*fft(h)); % 低通(低频分量)
% sig2=ifft(fft(y).*fft(g)); % 高通(高频分量)
% figure(5); % 信号图
% subplot(2,1,1)
% plot(real(sig1));
% title('分解信号1')
% subplot(2,1,2)
% plot(real(sig2));
% title('分解信号2')
% figure(6); % 频谱图
% subplot(2,1,1)
% stem(abs(fft(sig1)));
% title('分解信号1频谱')
% subplot(2,1,2)
% stem(abs(fft(sig2)));
% title('分解信号2频谱')
% %% 4.MALLET重构算法
% sig1=dyaddown(sig1); % 2抽取
% sig2=dyaddown(sig2); % 2抽取
% sig1=dyadup(sig1); % 2插值
% sig2=dyadup(sig2); % 2插值
% sig1=sig1(1,[1:N]); % 去掉最后一个零
% sig2=sig2(1,[1:N]); % 去掉最后一个零
% hr=h(end:-1:1); % 重构低通
% gr=g(end:-1:1); % 重构高通
% hr=circshift(hr',1)'; % 位置调整圆周右移一位
% gr=circshift(gr',1)'; % 位置调整圆周右移一位
% sig1=ifft(fft(hr).*fft(sig1)); % 低频
% sig2=ifft(fft(gr).*fft(sig2)); % 高频
% sig=sig1+sig2; % 源信号
% %% 5.比较
% figure(7);
% subplot(2,1,1)
% plot(real(sig1));
% title('重构低频信号');
% subplot(2,1,2)
% plot(real(sig2));
% title('重构高频信号');
% figure(8);
% subplot(2,1,1)
% stem(abs(fft(sig1)));
% title('重构低频信号频谱');
% subplot(2,1,2)
% stem(abs(fft(sig2)));
% title('重构高频信号频谱');
% figure(9)
% plot(real(sig),'r','linewidth',2);
% hold on;
% plot(y);
% legend('重构信号','原始信号')
% title('重构信号与原始信号比较')
%