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(-20*n*t);

 

 

figure(1);

 

plot(y);

title('原信號');

 

check1=sum(high_decompose);         %h0(n)性質校驗

check2=sum(low_decompose);

check3=norm(high_decompose);

check4=norm(low_decompose);

 

l_fre=conv(y,low_decompose);  %卷積

l_fre_down=dyaddown(l_fre);         %抽取,得低頻細節

 

h_fre=conv(y,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('小波分解的高頻系數');

 

 

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;          %信號重構

 

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');       %原信號與消噪后信號的比較

 

 

 

2．消噪



%  mallet_wavelet.m

%  此函數用于研究Mallet算法及濾波器設計

%  此函數用于消噪處理

 

%角度賦值

%此處賦值使濾波器系數恰為db9

 

%分解的高頻系數采用db9較好,即它的消失矩較大

%分解的有用信號小波高頻系數基本趨于零

%對于噪聲信號高頻分解系數很大，便于閾值消噪處理

 

[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');      %原信號與消噪后信號的比較 


小波譜分析mallat算法經典程序 
 
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('重構信號與原始信號比較')