小波降噪的函数

小波变换——压缩与去噪算法程序2009-10-04 14:27 最近忙着做方案设计，还有国庆的杂事比较多，也好久没上空间了。呵呵今天拿出来点东西晒晒，有关小波变换的一些应用，压缩和滤除噪声方面的。

＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃

ddencmp函数

% 当前延拓模式时补零
% 产生高斯白噪声
init = 2055415866; randn('seed',init);
x = randn(1,1000);
% 寻找使用小波降噪输出的3个默认值
disp('小波降噪输出的3个默认值')
[thr,sorh,keepapp] = ddencmp('den','wv',x)
% 寻找使用小波压缩输出的3个默认值
disp('小波压缩输出的3个默认值')
[thr,sorh,keepapp] = ddencmp('cmp','wv',x)
% 寻找使用小波包降噪输出的3个默认值
disp('小波包降噪输出的4个默认值')
[thr,sorh,keepapp,crit] = ddencmp('den','wp',x)
% 寻找使用小波包压缩输出的3个默认值
disp('小波包压缩输出的4个默认值')
[thr,sorh,keepapp,crit] = ddencmp('cmp','wp',x)

＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃

thselect函数

% 当前的延拓模式是补零
% 产生高斯白噪声
init = 2055415866; randn('seed',init);
x = randn(1,1000);
% 寻找各个阈值选择规则的阈值
% 使用SURE的自适应阈值
thr1 = thselect(x,'rigrsure')
% 固定形式阈值
thr2 = thselect(x,'sqtwolog')
% 启发式阈值
thr3 = thselect(x,'heursure')
% 极小极大阈值
thr4 = thselect(x,'minimaxi')
＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃

wbmpen函数

% 装载bumps的噪声信号
load noisbump; x = noisbump;
% 使用sym6执行5层小波分解
wname = 'sym6'; lev = 5;
[c,l] = wavedec(x,lev,wname);
% 由第1层的高频系数估计噪声标准差
sigma = wnoisest(c,l,1);
% 信号降噪时，采用调整参数选择全局阈值
alpha = 2;
thr = wbmpen(c,l,sigma,alpha)
% 使用软阈值和保存的低频信号，进行信号降噪
keepapp = 1;
xd = wdencmp('gbl',c,l,wname,lev,thr,'s',keepapp);
% 画出原始信号和降噪后的信号
figure(1)
subplot(211), plot(x), title('原始信号')
subplot(212), plot(xd), title('降噪后的信号')
＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃

wdcbm2函数

% 装载原始图像
load belmont2;
nbc = size(map,1);
% 使用sym4小波对图像进行3层分解
wname = 'sym4'; lev = 3;
[c,s] = wavedec2(X,lev,wname);
% 图像压缩，使用推荐参数选择各层不同的阈值
alpha = 1.5; m = 2.7*prod(s(1,:));
[thr,nkeep] = wdcbm2(c,s,alpha,m)
% 使用上面的阈值和硬阈值处理方式进行图像压缩
[xd,cxd,sxd,perf0,perfl2] =wdencmp('lvd',c,s,wname,lev,thr,'h');
% 画出原始图像和压缩后的图像
colormap(pink(nbc));
subplot(121), image(wcodemat(X,nbc)),
title('原始图像')
subplot(122), image(

wcodemat(xd,nbc)),
title('压缩后的图像')
disp('恢复百分比: ');
num2str(perfl2)
disp('压缩百分比: ');
num2str(perf0)
＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃

wdcbm函数

% 装载信号并选择其中一段
load leleccum; indx = 2600:3100;
x = leleccum(indx);
% 使用db3在第5层执行信号的小波分解
wname = 'db3'; lev = 5;
[c,l] = wavedec(x,lev,wname);
% 信号压缩，使用wdcbm 选择各层的独立阈值
alpha = 1.5; m = l(1);
[thr,nkeep] = wdcbm(c,l,alpha,m)
% 使用wdencmp和硬阈值进行信号压缩
[xd,cxd,lxd,perf0,perfl2] =wdencmp('lvd',c,l,wname,lev,thr,'h');
% 画出原始信号和压缩后信号
subplot(211), plot(indx,x), title('原始信号');
subplot(212), plot(indx,xd), title('压缩后信号');
xlab1 = ['2-norm rec.: ',num2str(perfl2)];
xlab2 = [' % -- zero cfs: ',num2str(perf0), ' %'];
xlabel([xlab1 xlab2]);
＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃

wdencmp函数

% 当前延拓模式是补零
% 装载原始图像
load sinsin
% X 包含原始图像
subplot(221),image(X),colormap(map);
title('原始图像');
% 产生噪声图像
init=2055615866; randn('seed',init);
x = X + 18*randn(size(X));
subplot(222),image(x),colormap(map);
title('噪声图像');
% 使用wdencmp进行图像降噪
% 寻找默认值
[thr,sorh,keepapp] = ddencmp('den','wv',x);
% 使用全局阈值进行图像降噪
xd = wdencmp('gbl',x,'sym4',2,thr,sorh,keepapp);
subplot(223),image(xd),colormap(map);
title('降噪后图像');
＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃

wden函数

%给定一个正弦信号并作图
t=0:1000;
s=sin(0.03*t);
subplot(3,1,1);
plot(s);
axis([0 1000 -2 2]);
title('原始信号');
%给该信号加噪声
load noissin;
ns=noissin;
subplot(3,1,2);
plot(ns);
title('染噪信号');
%进行降噪处理
xd=wden(ns,'minimaxi','s','one',5,'db3');
subplot(3,1,3);
plot(xd);
title('消噪信号');

＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃

wnoisest函数

% 当前延拓模式是补零
% 产生高斯白噪声
init = 2055415866; randn('seed',init);
x = randn(1,1000);
% 使用db3小波对x进行2层分解
[c,l] = wavedec(x,2,'db3');
% 在第1层和第2层估计系数的标准差
% 既然x是方差为1的高斯白噪声，那么估计值应接近1
sd1=wnoisest(c,l,1:2)
% 现假设x包含10个非正常值
ind = 50:50:500;
x(ind) = 100 * ones(size(ind));
% 使用db3小波对x进行1层分解
[ca,cd] = dwt(x,'db3');
% 普通的cd标准差估计方法
sd2=std(cd)
% cd标准差的鲁棒估计方法
sd3=median(abs(cd))/0.6745
＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃

wnoise函数

% 产生包含2^10点的'Heavy sine'信号
x = wnoise(3,10);
figure(1);
ind = linspace(0,1,2^10);
subplot(5,1,1),plot

(ind,x),title('Heavy sine信号');

% 产生包含2^10点的'Doppler' 信号，信噪比为7
[x,noisyx] = wnoise(4,10,7);
subplot(5,1,2),plot(ind,x),title('Doppler信号');
subplot(5,1,3),plot(ind,x),title('含噪声的Doppler信号');

% 产生随机数种子
init = 2055415866;
[x,noisyx] = wnoise(4,10,7,init);
subplot(5,1,2),plot(ind,x),title('Doppler');
subplot(5,1,4),plot(ind,x),title('由种子产生的Doppler信号');
subplot(5,1,5),plot(ind,x),title('由种子产生的含噪声的Doppler信号2');
% 显示所有的波形
figure(2);
x = wnoise(1,10);
subplot(6,1,1),plot(ind,x),title('Blocks信号');
x = wnoise(2,10);
subplot(6,1,2),plot(ind,x),title('Bumps信号');
x = wnoise(3,10);
subplot(6,1,3),plot(ind,x),title('Heavy sine信号');
x = wnoise(4,10);
subplot(6,1,4),plot(ind,x),title('Doppler信号');
x = wnoise(5,10);
subplot(6,1,5),plot(ind,x),title('HQuadchirp信号');
x = wnoise(6,10);
subplot(6,1,6),plot(ind,x),title('Mishmash信号');


＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃

wpbmpen函数

% 装载噪声chirp信号
load noischir; x = noischir;
% 使用sym6对信号进行5层小波包分解
wname = 'sym6'; lev = 5;
tree = wpdec(x,lev,wname);
% 由第1层高频系数估计噪声的标准差，结点索引号为2
det1 = wpcoef(tree,2);
sigma = median(abs(det1))/0.6745;
% 使用wpbmpen选择阈值
alpha = 2;
thr = wpbmpen(tree,sigma,alpha)
% 上面的阈值，软阈值处理，保存低频，使用wpdencmp进行信号降噪
keepapp = 1;
xd = wpdencmp(tree,'s','nobest',thr,keepapp);
% 画出原始信号和降噪后的信号
figure(1)
subplot(211), plot(x),
title('原始信号')
subplot(212), plot(xd)
title('降噪后信号')

＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃

wpdencmp函数

% 当前延拓模式是补零
% 装载原始信号
load sumlichr; x = sumlichr;
subplot(211), plot(x);
title('原始信号')
% 使用wpdencmp函数进行信号分解
% 寻找默认值（参见ddencmp函数）
[thr,sorh,keepapp,crit] = ddencmp('cmp','wp',x)
% 使用全局阈值进行压缩
[xc,treed,perf0,perfl2] =wpdencmp(x,sorh,3,'db2',crit,thr,keepapp);
subplot(212),plot(xc);
title('压缩后信号');
＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃

wthresh函数

% 产生信号，设置阈值
t = linspace(0,1,100);
y = 0.8*sin(2*pi*t);
subplot(1,3,1),plot(y);
title('原始信号'),grid on;
thr = 0.4;
% 进行硬阈值处理
ythard = wthresh(y,'h',thr);
subplot(1,3,2),plot(ythard );
title('硬阈值处理信号'),grid on;
% 进行软阈值处理
ytsoft = wthresh(y,'s',thr);
subplot(1,3,3),plot(ytsoft);
title('软阈值处理信号'),grid on;
＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃

小波包压缩函数

%装载并显示原始图像
load belmont2;
subplot(1,2,1);
image(X);c

olormap(map);title('原始图像');
%采用默认的全局阈值
[thr,sorh,keepapp,crit]=ddencmp('cmp','wp',X);
%图像压缩
[Xc,treed,perf0,perfl2]=wpdencmp(X,sorh,3,'bior3.1',crit,thr,keepapp);
subplot(1,2,2);image(Xc);colormap(map);title('压缩后的图像');
%给出压缩效率
disp('小波分解系数中置0的系数个数百分比：');
perf0
disp('压缩后图像剩余能量百分比：');
perfl2