load wbarb matlab,matlab小波变换程序

MATLAB2维小波变换经典程序

% FWT_DB.M;

% 此示意程序用DWT实现二维小波变换

% 编程时间2004-4-10，编程人沙威

%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;clc;

T=256; % 图像维数

SUB_T=T/2; % 子图维数

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 1.调原始图像矩阵

load wbarb; % 下载图像

f=X; % 原始图像

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 2.进行二维小波分解

l=wfilters('db10','l'); % db10(消失矩为10)低通分解滤波器冲击响应(长度为20)

L=T-length(l);

l_zeros=[l,zeros(1,L)]; % 矩阵行数与输入图像一致，为2的整数幂

h=wfilters('db10','h'); % db10(消失矩为10)高通分解滤波器冲击响应(长度为20)

h_zeros=[h,zeros(1,L)]; % 矩阵行数与输入图像一致，为2的整数幂

for i=1:T; % 列变换

row(1:SUB_T,i)=dyaddown( ifft( fft(l_zeros).*fft(f(:,i)') ) ).'; % 圆周卷积FFT

row(SUB_T+1:T,i)=dyaddown( ifft( fft(h_zeros).*fft(f(:,i)') ) ).'; % 圆周卷积FFT

end;

for j=1:T; % 行变换

line(j,1:SUB_T)=dyaddown( ifft( fft(l_zeros).*fft(row(j,:)) ) ); % 圆周卷积FFT

line(j,SUB_T+1:T)=dyaddown( ifft( fft(h_zeros).*fft(row(j,:)) ) ); % 圆周卷积FFT

end;

decompose_pic=line; % 分解矩阵

% 图像分为四块

lt_pic=decompose_pic(1:SUB_T,1:SUB_T); % 在矩阵左上方为低频分量--fi(x)*fi(y)

rt_pic=decompose_pic(1:SUB_T,SUB_T+1:T); % 矩阵右上为--fi(x)*psi(y)

lb_pic=decompose_pic(SUB_T+1:T,1:SUB_T); % 矩阵左下为--psi(x)*fi(y)

rb_pic=decompose_pic(SUB_T+1:T,SUB_T+1:T); % 右下方为高频分量--psi(x)*psi(y)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 3.分解结果显示

figure(1);

colormap(map);

subplot(2,1,1);

image(f); % 原始图像 title('original pic');

subplot(2,1,2);

image(abs(decompose_pic)); % 分解后图像

title('decomposed pic');

figure(2);

colormap(map);

subplot(2,2,1);

image(abs(lt_pic)); % 左上方为低频分量--fi(x)*fi(y)

title('\Phi(x)*\Phi(y)');

subplot(2,2,2);

image(abs(rt_pic)); % 矩阵右上为--fi(x)*psi(y)

title('\Phi(x)*\Psi(y)');

subplot(2,2,3);

image(abs(lb_pic)); % 矩阵左下为--psi(x)*fi(y)

title('\Psi(x)*\Phi(y)');

subplot(2,2,4);

image(abs(rb_pic)); % 右下方为高频分量--psi(x)*psi(y)

title('\Psi(x)*\Psi(y)');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 5.重构源图像及结果显示

% construct_pic=decompose_matrix'*decompose_pic*decompose_matrix;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

l_re=l_zeros(end:-1:1); % 重构低通滤波

l_r=circshift(l_re',1)'; % 位置调整

h_re=h_zeros(end:-1:1); % 重构高通滤波

h_r=circshift(h_re',1)'; % 位置调整

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

top_pic=[lt_pic,rt_pic]; % 图像上半部分

t=0;

for i=1:T; % 行插值低频

if (mod(i,2)==0)

topll(i,:)=top_pic(t,:); % 偶数行保持

else

t=t+1;

topll(i,:)=zeros(1,T); % 奇数行为零

end

end;

for i=1:T; % 列变换

topcl_re(:,i)=ifft( fft(l_r).*fft(topll(:,i)') )'; % 圆周卷积FFT

end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

bottom_pic=[lb_pic,rb_pic]; % 图像下半部分

t=0;

for i=1:T; % 行插值高频

if (mod(i,2)==0)

bottomlh(i,:)=bottom_pic(t,:); % 偶数行保持

else

bottomlh(i,:)=zeros(1,T); % 奇数行为零

t=t+1;

end

end;

for i=1:T; % 列变换

bottomch_re(:,i)=ifft( fft(h_r).*fft(bottomlh(:,i)') )'; % 圆周卷积FFT

end;

construct1=bottomch_re+topcl_re; % 列变换重构完毕

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

left_pic=construct1(:,1:SUB_T); % 图像左半部分

t=0;

for i=1:T; % 列插值低频

if (mod(i,2)==0)

leftll(:,i)=left_pic(:,t); % 偶数列保持

else

t=t+1;

leftll(:,i)=zeros(T,1); % 奇数列为零

end

end;

for i=1:T; % 行变换

leftcl_re(i,:)=ifft( fft(l_r).*fft(leftll(i,:)) ); % 圆周卷积FFT

end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

right_pic=construct1(:,SUB_T+1:T); % 图像右半部分

t=0;

for i=1:T; % 列插值高频

if (mod(i,2)==0)

rightlh(:,i)=right_pic(:,t); % 偶数列保持

else

rightlh(:,i)=zeros(T,1); % 奇数列为零

t=t+1;

end

end;

for i=1:T; % 行变换

rightch_re(i,:)=ifft( fft(h_r).*fft(rightlh(i,:)) ); % 圆周卷积FFT

end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

construct_pic=rightch_re+leftcl_re; % 重建全部图像

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 结果显示

figure(3);

colormap(map);

subplot(2,1,1);

image(f); % 源图像显示

title('original pic');

subplot(2,1,2);

image(abs(construct_pic)); % 重构源图像显示

title('reconstructed pic');

error=abs(construct_pic-f); % 重构图形与原始图像误值

figure(4);

mesh(error); % 误差三维图像

title('absolute error display');

clear

clc

%在噪声环境下语音信号的增强

%语音信号为读入的声音文件

%噪声为正态随机噪声

sound=wavread('c12345.wav');

count1=length(sound);

noise=0.05*randn(1,count1);

for i=1:count1

signal(i)=sound(i);

end

for i=1:count1

y(i)=signal(i)+noise(i);

end

%在小波基'db3'下进行一维离散小波变换

[coefs1,coefs2]=dwt(y,'db3'); %[低频 高频]

count2=length(coefs1);

count3=length(coefs2);

energy1=sum((abs(coefs1)).^2);

energy2=sum((abs(coefs2)).^2);

energy3=energy1+energy2;

for i=1:count2

recoefs1(i)=coefs1(i)/energy3;

end

for i=1:count3

recoefs2(i)=coefs2(i)/energy3;

end

%低频系数进行语音信号清浊音的判别

zhen=160;

count4=fix(count2/zhen);

for i=1:count4

n=160*(i-1)+1:160+160*(i-1);

s=sound(n);

w=hamming(160);

sw=s.*w;

a=aryule(sw,10);

sw=filter(a,1,sw);

sw=sw/sum(sw);

r=xcorr(sw,'biased');

corr=max(r);

%为清音(unvoice)时，输出为1；为浊音(voice)时，输出为0

if corr>=0.8

output1(i)=0;

elseif corr<=0.1

output1(i)=1;

end

end

for i=1:count4

n=160*(i-1)+1:160+160*(i-1);

if output1(i)==1

switch abs(recoefs1(i))

case abs(recoefs1(i))<=0.002

recoefs1(i)=0;

case abs(recoefs1(i))>0.002 &

abs(recoefs1(i))<=0.003

recoefs1(i)=sgn(recoefs1(i))*(0.003*abs(recoefs1(i))-0.000003)/0.002;

otherwise recoefs1(i)=recoefs1(i);

end

elseif output1(i)==0

recoefs1(i)=recoefs1(i);

end

end

%对高频系数进行语音信号清浊音的判别

count5=fix(count3/zhen);

for i=1:count5

n=160*(i-1)+1:160+160*(i-1);

s=sound(n);

w=hamming(160);

sw=s.*w;

a=aryule(sw,10);

sw=filter(a,1,sw);

sw=sw/sum(sw);

r=xcorr(sw,'biased');

corr=max(r);

%为清音(unvoice)时，输出为1；为浊音(voice)时，输出为0

if corr>=0.8

output2(i)=0;

elseif corr<=0.1

output2(i)=1;

end

end

for i=1:count5

n=160*(i-1)+1:160+160*(i-1);

if output2(i)==1

switch abs(recoefs2(i))

case abs(recoefs2(i))<=0.002

recoefs2(i)=0;

case abs(recoefs2(i))>0.002 &

abs(recoefs2(i))<=0.003

recoefs2(i)=sgn(recoefs2(i))*(0.003*abs(recoefs2(i))-0.000003)/0.002;

otherwise recoefs2(i)=recoefs2(i);

end

elseif output2(i)==0

recoefs2(i)=recoefs2(i);

end

end

%在小波基'db3'下进行一维离散小波反变换 output3=idwt(recoefs1, recoefs2,'db3');

%对输出信号抽样点值进行归一化处理

maxdata=max(output3);

output4=output3/maxdata;

%读出带噪语音信号，存为'101.wav'

wavwrite(y,5500,16,'c101'); %读出处理后语音信号，存为'102.wav'

wavwrite(output4,5500,16,'c102');

function [I_W , S] = func_DWT(I, level, Lo_D, Hi_D);

%通过这个函数将I进行小波分解，并将分解后的一维向量转换为矩阵形式

% Matlab implementation of SPIHT (without Arithmatic coding

stage)

% Wavelet decomposition

% input: I :

input image

% level : wavelet decomposition level

% Lo_D : low-pass decomposition filter

% Hi_D : high-pass decomposition filter

% output: I_W : decomposed

image vector

% S : corresponding bookkeeping matrix

% please refer wavedec2 function to see

more [C,S] = func_Mywavedec2(I,level,Lo_D,Hi_D); S(:,3) =

S(:,1).*S(:,2); % dim of detail coef nmatrices 求低频和每个尺度中高频的元素个数

%st=S(1,3)+S(2,3)*3+S(3,3)*3;%%%%对前两层加密

%C(1:st)=0;

L = length(S); %a求S的列数

I_W = zeros(S(L,1),S(L,2));%设一个与原图像大小相同的全零矩阵

% approx part

I_W( 1:S(1,1) , 1:S(1,2) ) = reshape(C(1:S(1,3)),S(1,1:2));

%将LL层从C中还原为S(1,1)*S(1,2)的矩阵

for k = 2 : L-1 %将C向量中还原出HL,HH,LH 矩阵

rows =

[sum(S(1:k-1,1))+1:sum(S(1:k,1))];

columns =

[sum(S(1:k-1,2))+1:sum(S(1:k,2))];

% horizontal

part

c_start =

S(1,3) + 3*sum(S(2:k-1,3)) + 1;

c_stop =

S(1,3) + 3*sum(S(2:k-1,3)) + S(k,3);

I_W(

1:S(k,1) , columns ) = reshape( C(c_start:c_stop) , S(k,1:2)

);

% vertical

part

c_start =

S(1,3) + 3*sum(S(2:k-1,3)) + S(k,3) + 1;

c_stop =

S(1,3) + 3*sum(S(2:k-1,3)) + 2*S(k,3);

I_W( rows ,

1:S(k,2) ) = reshape( C(c_start:c_stop) , S(k,1:2) );

% diagonal

part

c_start =

S(1,3) + 3*sum(S(2:k-1,3)) + 2*S(k,3) + 1;

c_stop =

S(1,3) + 3*sum(S(2:k,3));

I_W( rows ,

columns ) = reshape( C(c_start:c_stop) , S(k,1:2) );

end

%%%%%%%mallat algorithm%%%%%

clc;

clear;tic;

%%%%original signal%%%%

f=100;%%frequence

ts=1/800;%%抽样间隔

N=1:100;%%点数

s=sin(2*ts*pi*f.*N);%%源信号

figure(1)

plot(s);%%%源信号s

title('原信号');

grid on;

%%%%小波滤波器%%%%

ld=wfilters('db1','l');%%低通

hd=wfilters('db1','h');%%高通

figure(2)

stem(ld,'r');%%%低通

grid on;

figure(3)

stem(hd,'b')%%%高通

grid on;

%%%%%

tem=conv(s,ld);%%低通和原信号卷积

ca1=dyaddown(tem);%%抽样

figure(4)

plot(ca1);

grid on;

tem=conv(s,hd);%%高通和原信号卷积

cb1=dyaddown(tem);%%抽样

figure(5)

plot(cb1);

grid on;

%%%%%%%%

%[ca3,cb3]=dwt(s,'db1');%%小波变换

%%%%%%%%

[lr,hr]=wfilters('db1','r');%%重构滤波器

figure(6)

stem(lr);

figure(7)

stem(hr);

tem=dyadup(cb1);%%插值

tem=conv(tem,hr);%%卷积

d1=wkeep(tem,100);%%去掉两头的分量

%%%%%%%%%

tem=dyadup(ca1);%%插值

tem=conv(tem,lr);%%卷积

a1=wkeep(tem,100);%%去掉两头的分量

a=a1+d1;%%%重构原信号

%%%%%%%%%

%a3=idwt(ca3,cb3,'db1',100);%%%小波逆变换

%%%%%%%%%

figure(8)

plot(a,'.b');

hold on;

plot(s,'r');

grid on;

title('重构信号和原信号的比较');toc;

%figure(9)

%plot(a3,'.b');

%hold on;

%plot(s,'r');

%grid on;

%title('重构信号和原信号的比较');

通用函数 　 Allnodes 计算树结点 appcoef 提取一维小波变换低频系数 appcoef2 提取二维小波分解低频系数 bestlevt 计算完整最佳小波包树 besttree 计算最佳(优)树 *　 biorfilt 双正交样条小波滤波器组 biorwavf 双正交样条小波滤波器 *　 centfrq 求小波中心频率 cgauwavf Complex Gaussian小波 cmorwavf coiflets小波滤波器 cwt 一维连续小波变换 dbaux Daubechies小波滤波器计算 dbwavf Daubechies小波滤波器 dbwavf(W) W='dbN' N=1,2,3,...,50 ddencmp 获取默认值阈值(软或硬)熵标准 depo2ind 将深度-位置结点形式转化成索引结点形式 detcoef 提取一维小波变换高频系数 detcoef2 提取二维小波分解高频系数 disp 显示文本或矩阵 drawtree 画小波包分解树(GUI) dtree 构造DTREE类 dwt 单尺度一维离散小波变换 dwt2 单尺度二维离散小波变换 dwtmode 离散小波变换拓展模式 *　 dyaddown 二元取样 *　 dyadup 二元插值 entrupd 更新小波包的熵值 fbspwavf B样条小波 gauswavf Gaussian小波 get 获取对象属性值 idwt 单尺度一维离散小波逆变换 idwt2 单尺度二维离散小波逆变换 ind2depo 将索引结点形式转化成深度—位置结点形式 *　 intwave 积分小波数 isnode 判断结点是否存在 istnode 判断结点是否是终结点并返回排列值 iswt 一维逆SWT(Stationary Wavelet Transform)变换 iswt2 二维逆SWT变换 leaves Determine terminal nodes

mexihat 墨西哥帽小波 meyer Meyer小波 meyeraux Meyer小波辅助函数 morlet Morlet小波 nodease 计算上溯结点 nodedesc 计算下溯结点(子结点) nodejoin 重组结点 nodepar 寻找父结点 nodesplt 分割(分解)结点 noleaves Determine nonterminal nodes

ntnode Number of terminal nodes

ntree Constructor for the class NTREE

*　 orthfilt 正交小波滤波器组 plot 绘制向量或矩阵的图形 *　 qmf 镜像二次滤波器 rbiowavf Reverse biorthogonal spline wavelet filters

read 读取二进制数据 readtree 读取小波包分解树 *　 scal2frq Scale to frequency

set shanwavf Shannon wavelets

swt 一维SWT(Stationary Wavelet Transform)变换 swt2 二维SWT变换 symaux Symlet wavelet filter computation.

symwavf Symlets小波滤波器 thselect 信号消噪的阈值选择