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‘)