matlab小波变换程序

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)') ) ).'; % 圆周卷积<->FFTend;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,:)) ) ); % 圆周卷积<->FFTend;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); % 奇数行为零 endend;for i=1:T; % 列变换 topcl_re(:,i)=ifft( fft(l_r).*fft(topll(:,i)') )'; % 圆周卷积<->FFTend; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%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; endend;for i=1:T; % 列变换 bottomch_re(:,i)=ifft( fft(h_r).*fft(bottomlh(:,i)') )'; % 圆周卷积<->FFTend; 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); % 奇数列为零 endend;for i=1:T; % 行变换 leftcl_re(i,:)=ifft( fft(l_r).*fft(leftll(i,:)) ); % 圆周卷积<->FFTend; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%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; endend;for i=1:T; % 行变换 rightch_re(i,:)=ifft( fft(h_r).*fft(rightlh(i,:)) ); % 圆周卷积<->FFTend;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%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'); clearclc%在噪声环境下语音信号的增强%语音信号为读入的声音文件%噪声为正态随机噪声sound=wavread('c12345.wav');count1=length(sound);noise=0.05*randn(1,count1);for i=1:count1signal(i)=sound(i);endfor i=1:count1y(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:count2recoefs1(i)=coefs1(i)/energy3;endfor i=1:count3recoefs2(i)=coefs2(i)/energy3;end%低频系数进行语音信号清浊音的判别zhen=160;count4=fix(count2/zhen);for i=1:count4n=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）时，输出为0if corr>=0.8output1(i)=0;elseif corr<=0.1output1(i)=1;endendfor i=1:count4n=160*(i-1)+1:160+160*(i-1);if output1(i)==1switch abs(recoefs1(i))case abs(recoefs1(i))<=0.002recoefs1(i)=0;case abs(recoefs1(i))>0.002 & abs(recoefs1(i))<=0.003recoefs1(i)=sgn(recoefs1(i))*(0.003*abs(recoefs1(i))-0.000003)/0.002;otherwise recoefs1(i)=recoefs1(i);endelseif output1(i)==0recoefs1(i)=recoefs1(i);endend%对高频系数进行语音信号清浊音的判别count5=fix(count3/zhen);for i=1:count5n=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）时，输出为0if corr>=0.8output2(i)=0;elseif corr<=0.1output2(i)=1;endendfor i=1:count5n=160*(i-1)+1:160+160*(i-1);if output2(i)==1switch abs(recoefs2(i))case abs(recoefs2(i))<=0.002recoefs2(i)=0;case abs(recoefs2(i))>0.002 & abs(recoefs2(i))<=0.003recoefs2(i)=sgn(recoefs2(i))*(0.003*abs(recoefs2(i))-0.000003)/0.002;otherwise recoefs2(i)=recoefs2(i);endelseif output2(i)==0recoefs2(i)=recoefs2(i);endend%在小波基'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 partI_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 nodesmexihat   墨西哥帽小波  meyer   Meyer小波  meyeraux   Meyer小波辅助函数  morlet   Morlet小波  nodease   计算上溯结点  nodedesc   计算下溯结点(子结点)  nodejoin   重组结点  nodepar   寻找父结点  nodesplt   分割(分解)结点  noleaves   　 Determine nonterminal nodesntnode   　 Number of terminal nodesntree   　 Constructor for the class NTREE *　 orthfilt   正交小波滤波器组  plot   绘制向量或矩阵的图形  *　 qmf   镜像二次滤波器  rbiowavf   　 Reverse biorthogonal spline wavelet filtersread   读取二进制数据  readtree   读取小波包分解树  *　 scal2frq   　 Scale to frequencyset   　 shanwavf   　 Shannon waveletsswt   一维SWT(Stationary Wavelet Transform)变换  swt2   二维SWT变换  symaux   　 Symlet wavelet filter computation.symwavf   Symlets小波滤波器  thselect   信号消噪的阈值选择  thodes   　 Referencestreedpth   求树的深度  treeord   求树结构的叉数   upcoef   一维小波分解系数的直接重构  upcoef2   二维小波分解系数的直接重构  upwlev   单尺度一维小波分解的重构  upwlev2   单尺度二维小波分解的重构  wavedec   单尺度一维小波分解  wavedec2   多尺度二维小波分解  wavedemo   小波工具箱函数demo  *　wavefun   小波函数和尺度函数  *　 wavefun2   二维小波函数和尺度函数  wavemenu   小波工具箱函数menu图形界面调用函数  *　 wavemngr   小波管理函数  waverec   多尺度一维小波重构  waverec2   多尺度二维小波重构  wbmpen   　 Penalized threshold for wavelet 1-D or 2-D de-noisingwcodemat   对矩阵进行量化编码  wdcbm   　 Thresholds for wavelet 1-D using Birge-Massart strategywdcbm2   　Thresholds for wavelet 2-D using Birge-Massart strategy wden   用小波进行一维信号的消噪或压缩  wdencmp   　De-noising or compression using wavelets wentropy   计算小波包的熵  wextend   　Extend a vector or a matrix *　 wfilters   小波滤波器  wkeep   提取向量或矩阵中的一部分  *　 wmaxlev   计算小波分解的最大尺度  wnoise   产生含噪声的测试函数数据  wnoisest   估计一维小波的系数的标准偏差  wp2wtree   从小波包树中提取小波树    　 wpcoef   计算小波包系数  wpcutree   剪切小波包分解树  wpdec   一维小波包的分解  wpdec2   二维小波包的分解  wpdencmp   用小波包进行信号的消噪或压缩  wpfun   小波包函数  wpjoin   　重组小波包 wprcoef   小波包分解系数的重构  wprec   一维小波包分解的重构  wprec2   二维小波包分解的重构  wpsplt   分割（分解）小波包  wpthcoef   进行小波包分解系数的阈值处理  wptree   　 显示小波包树结构wpviewcf   　 Plot the colored wavelet packet coefficients. wrcoef   对一维小波系数进行单支重构  wrcoef2   对二维小波系数进行单支重构  wrev   向量逆序  write   向缓冲区内存写进数据  wtbo   　 Constructor for the class WTBO wthcoef   一维信号的小波系数阈值处理  wthcoef2   二维信号的小波系数阈值处理  wthresh   进行软阈值或硬阈值处理  wthrmngr   阈值设置管理  wtreemgr   管理树结构