matlab信号小波分解与重构入门 - 那抹阳光1994

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