4-1 Mallet小波分解重構程序

Mallet小波在小波屆的地位相似fft在傅立葉變化中的地位，在分解過程當中先濾波後抽取，重構過程當中先插值後濾波，能夠操做正交小波變換和雙正交小波變換。

本文中的程序是對構造的信號進行高低通濾波，以後再進行高低頻重構，實現matlab中Mallet小波的基本操做。

%% 本程序採用Mallat算法對信號進行小波分析

clc;clear;
close all;

% 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)
subplot(2,1,1)
plot(y);
title('原始信號');
subplot(2,1,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(2);
subplot(2,1,1)
stem(abs(fft(h)));
title('分解低通濾波器頻譜');

subplot(2,1,2)
stem(abs(fft(g)));
title('分解高通濾波器頻譜');

% 3 Mallet分解算法（圓周卷積的快速傅立葉變換實現）

sig1= ifft(fft(y).*fft(h));	%低通
sig2= ifft(fft(y).*fft(g));	%高通

figure(3)	% 信號圖
subplot(2,2,1)
plot(real(sig1));
title('低頻份量')

subplot(2,2,2)
plot(real(sig2));
title('高頻份量')

subplot(2,2,3)
stem(abs(fft(sig1)));
title('低頻份量頻譜')

subplot(2,2,4)
stem(abs(fft(sig2)));
title('高頻份量頻譜')

% 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(4);
subplot(2,2,1);
plot(real(sig1));
title('重構低頻信號');
subplot(2,2,2);
plot(real(sig2));
title('重構高頻信號');
subplot(2,2,3);
stem(abs(fft(sig1)));
title('重構低頻信號頻譜');
subplot(2,2,4);
stem(abs(fft(sig2)));
title('重構高頻信號頻譜');

figure(5)
plot(real(sig),'r','linewidth',2);
hold on
plot(y)
legend('重構信號','原始信號');
title('重構信號與原始信號比較')

程序運行結果

  

 

問題：

1 wfilters函數

    功能：生成四種小波濾波器，高低通濾波器，高低頻重構濾波器

    格式：[LO_D,HI_D,LO_R,HI_R] = WFILTERS('wname')

              [F1,F2] = WFILTERS('wname','type')

    說明：LO_D - 低通濾波器

              HI_D - 高通濾波器

              LO_R - 低通重構濾波器

              HI_R - 高通重構濾波器

              type：

               LO_D and HI_D if 'type' = 'd' (分解)

               LO_R and HI_R if 'type' = 'r' (重構)

               LO_D and LO_R if 'type' = 'l' (低通濾波器)

               HI_D and HI_R if 'type' = 'h' (高通濾波器)    

               wname：

               wname說明

     示例：使用中若是信號長度和濾波器的長度不一致，須要將濾波器補零

h= wfilters('db30', 'l');	%	低通
g= wfilters('db30', 'h');	%	高通

h= [h,zeros(1,N-length(h))];	%補零（圓周卷積，且增大分辨率便於觀察）
g= [g,zeros(1,N-length(g))];	%補零（圓周卷積，且增大分辨率便於觀察）

figure(2);
subplot(2,1,1)
stem(abs(fft(h)));
title('分解低通濾波器頻譜');

subplot(2,1,2)
stem(abs(fft(g)));
title('分解高通濾波器頻譜');

2 圓周卷積

線性卷積就是多項式係數乘法：設a的長度是M，b的長度是N，則a卷積b的長度是M+N-1，運算參見多項式乘法。
「L點的圓周卷積」就是把先作線性卷積，再把結果的前L點保留不動，後面的點截下來，加到結果的頭上去。若是L>M+N-1，則線性卷積和圓周卷積相同。----來自百度知道

圓周卷積 線性卷積 週期卷積

3 抽樣函數

dyaddown

功能：對時間序列進行二元採樣，每隔一個元素提取一個元素，獲得一個降採樣時間序列。

格式：

1.y = dyaddown(x, EVENODD)

當EVENODD=0時，從x中第二個元素開始採樣（偶採樣）；當EVENODD=1時，從x中第一個元素開始採樣（奇採樣）。

2.y = dyaddown(x)

EVENODD缺省，按EVENODD=0

 

dyadup

功能：對時間序列進行二元插值，每隔一個元素插入一個0元素，獲得一個時間序列。

格式：

1.y = dyadup(x, EVENODD)

當EVENODD=0時，從x中第二個元素開始採樣（偶採樣）；當EVENODD=1時，從x中第一個元素開始採樣（奇採樣）。

2.y = dyadup(x)

EVENODD缺省，按EVENODD=0