matlab信號小波分解與重構入門

 

%% 1. 利用MATLAB產生分解與重構濾波器組
% [Ld, Hd, Lr, Hr] = wfilters(wn);
% wfname:小波名
% Ld:分解低通濾波器h0[-n];
% Hd:分解高通濾波器h1[-n];
% Lr:分解低通濾波器h0[-n];
% Hr:分解高通濾波器h1[-n];

% eg1:計算db2小波的四個濾波器，並畫出其時域波形。
wn='db2'; % wfname:小波名
[Ld, Hd, Lr, Hr] = wfilters(wn);
k=0:3;

figure
subplot(221);
stem(k, Ld);
title('低通分解濾波器Ld');

subplot(222);
stem(k, Lr);
title('低通重建濾波器Lr');

subplot(223);
stem(k, Hd);
title('高通分解濾波器Hd');

subplot(224);
stem(k, Hd);
title('高通重建濾波器Hr');


%% 2. 利用MATLAB計算小波函數
% [phi,psi,t]=wavefun('wname',lter)
% wname:小波名
% lter:計算過程的迭代次數
% phi：尺度函數ϕ(t)
% psi:小波函數ψ(t)
% t:尺度函數與小波函數的抽樣點

% eg2:利用MATLAB計算db2小波的尺度函數與小波函數。
wname='db2';
lter=10;

[phi, psi, t]=wavefun(wname, lter);

figure
subplot(211);
plot(t, phi);
title('尺度函數')

subplot(212);
plot(t, psi);
title('小波函數')

%% 3. 利用MATLAB計算一維DWT和IDWT
% [C, L] = wavedec(x, N, 'wname');
% 
% x = waverec(C, L, 'wname');

% wname： 小波名; 
% x: 時域信號;
% N： 小波變換/分解的級數;
% C = [cAN, cDN, cDN-1, ..., cD1]； 
% cAN為近似分量， cDN, cDN-1, ..., cD1為細節分量
% L(1)= cAN的長度， L(2)= cDN的長度， L(N+1)= cD1的長度， L(N+2)= x的長度

% 小波分解
wname='haar';
N = 5;
[C, L] = wavedec(data(:, 1), N, wname);

figure
plot(C)  % [1024, 1]

% L = 32, [32, 64, 128, 256, 512], 1024

% 信號重構
x = waverec(C, L, wname);

figure
plot(data(:, 1), 'k')
hold on
plot(x, 'r')
xlim([0, 1024])
legend('Original', 'Reconstructed')


%% 4 利用部分小波系數重建信號

% x=wrcoef('type', C, L, 'wname', N)
% type='a' 由第N級近似分量重建信號
% type='d' 由第N級細節分量重建信號
% wname： 小波名
% 小波分解是樹狀二叉分解，[cA3 cD3] == [cA2]
% x -> [cA1, cD1] -> [[cA2, cD2], cD1] -> [[[cA3, cD3], cD2], cD1]
% 若 C = [cA3 cD3 cD2 cD1]
% x=wrcoef('a',C,L, 'wname',3)
% x=IDWT{[cA3 0 0 0]}
% x=wrcoef('a',C,L, , 'wname',2)
% x=IDWT{[cA3 cD3 0 0] }=IDWT{[cA2 0 0]}


% eg3 已知某信號的波形，試計算其5級小波變換系數，
% 分別由第5、3、1級小波近似系數重建信號。

% 小波參數
wname ='db1';
% Change DWT Extension Mode 
dwtmode('per') % DWT Extension Mode: Periodization

% 信號
% t=linspace(0, 1, 1024);
% x=20*t.^2.*(1-t).^4.*cos(12*pi*t);
fs = 128;
dt = 1/fs;
t = 0:dt:(1024/128)-dt;
x = data(:, 1);

figure
subplot(511);
plot(t, x);
%axis([0 1 -0.5 0.5]);
title('Signal');

% 5級分解
[C, L] = wavedec(x, 5, wname);
subplot(512);
plot(t, C); 
%axis([0 1 -3 3]);

% 由第5級小波近似系數重建信號
A5 = wrcoef('a', C, L, wname, 5);
subplot(513);
plot(t, A5);
%axis([0 1 -0.5 0.5]);

% 由第3級小波近似系數重建信號
A3 = wrcoef('a', C, L, wname, 3);
subplot(514);
plot(t, A3);
%axis([0 1 -0.5 0.5]);

% 由第1級小波近似系數重建信號
A1 = wrcoef('a', C, L, wname, 1);
subplot(515);
plot(t, A1);
%axis([0 1 -0.5 0.5]);


%% 5. 基於小波的信號去噪

% XD = wden(X, TPTR, SORH, SCAL, N, 'wname')
% 其中：
% XD: 對噪聲信號X去噪后得到的信號；
% X: 含噪聲信號;
% TPTR: 閾值規則，主要有'rigrsure', 'heursure', 'sqtwolog', 'minimaxi';
% SORH: 閾值方法, 's' (soft閾值)， 'h' (hard閾值)；
% SCAL: 閾值尺度的調整方法，主要有'one', ' sln', ' mln' ；
% N: 離散小波變換的級數
% wname: 小波名

% eg4 試利用函數wden對含有噪聲的blocks信號進行去噪。
snr = 5; % 噪聲方差

[x, xn] = wnoise('blocks', 11, snr);
k = 0:length(x)-1;

figure
subplot(311);
plot(k, x);
title('原信號');

subplot(312);
plot(k, xn);
title('含噪信號');

% 利用soft SURE閾值規則去噪
lev = 5;
wn = 'db1';
xd1 = wden(xn, 'heursure', 's', 'one', lev, wn);
subplot(313);
plot(k, xd1);
title('去噪信號');

%% 6. 基於小波的信號壓縮
% NC= wthcoef('d', C, L, N)
% 其中：
% 'd'： 表示對DWT系數C中細節(detail)分量進行壓縮;
% C,L： 由wavedec得到的DWT系數；
%   N： 若N=[1 2 3]表示將C中1、2和3級細節分量置零。
%  NC： 由系數C經過壓縮后得到的新系數;

% eg5試利用函數wthcoef對leleccum信號進行壓縮。

load leleccum

x = leleccum(1001:2024)*0.95/100;
k=0:length(x)-1;

[C, L] = wavedec(x, 5, 'db3');

NC1 = wthcoef('d', C, L, [1]); % 壓縮比率512/1024
x1 = waverec(NC1, L, 'db3');

NC5 = wthcoef('d', C, L, [1, 2, 3, 4, 5]); % 壓縮比率32/1024
x2 = waverec(NC5, L, 'db3');

figure
subplot(311);
plot(k, x);
title('原信號');

subplot(312);
plot(k, x1);
title('512/1024壓縮');

subplot(313);
plot(k, x2);
title('32/1024壓縮');


size(C)
size(NC1)
size(NC5)

figure
plot(C, 'k')
hold on
plot(NC1, 'r')
hold on
plot(NC5, 'b')

　　

 

參考：

https://blog.csdn.net/qq_24163555/article/details/82827187