小波分析实验：维离散小波变换

实验目的:在理解离散小波变换原理和 Mallat快速算法的基础上，通 过编程对图像进行二维离散小波变换，从而加深对二维小波分解 和重构的理性和感性认识，并能提高编程能力，为今后的学习和 工作奠定基础。实验工具：计算机，matlab6.5重构算法二讨二工歡肖一十三畧3 - 2k)dk *分解算法写成矩阵的形式 Cli, g的长度为4).汕叩4-弹0粗为疑习0-00 乙0 _00闊 1A2A3-0加卩1闫2i殆0190 0B0机0|i殆L-勺 _1山一1.诃 砒【1*呂【00目1 S2 g3 0-0T.r加卫_01曲&【1吕【2 g3I-ia駅T_g2由00 0 0g。glj鼻【0血1吃h3 000 莎记010s1吐=S1凰1】疑2 * 上13 110即】1虫2咎00 00应0址1勺堆-1.g0gl审g3 000 心0】s2=0-90BI煎01-或1烈习 川B4109D产心1圍2护00 00或0g【】L如-1于是Malht分解公式为矩阵变换！丄重构算法写成矩阵变换：三三出片y 峰值信噪比计算公式：P碗二 Q产f ME豆i MSE（4）.（4-）分别表示原始圏像和重建圈像.本实验采取的_些外技巧；卍（1）分蹄法：0凤1243 0 i| 0 0 T0!也A23 0r： q?1砒曙0 00-0 A0 |1【0注刃2fi3 00叮10和规1g2g301:：-盘1g2即0 00-0 g0J|編程时用如下思想，（h, g的长度为4）4如0】勺t 1将数据Up和低通（高通）滤波盎进行添零到数据长度出，再2抽样，相卽得到小波后一半系数“內插尋咼適 重构 滤波 器重构结果編程时用如下思想0心1*IV川（）o*4-0000000*4-4000G0“3*iioTgIT*t00 0 00 曲0 00* 4 * * 4s* 0 0 00004 *或4*或1肩_ 如。1 汕1=:由【151+一 0 沖1i-1H/J-lH-l实验结果：屮圉tl）三纽小谀分解匿10015020C25010015010015G 200 25JE1C2 7原嫦图爍和重构图慷峰值信噪比三psnr-252.8923lb附录:（1）二维小波分解函数%二维小波分解函数fun cti on Y=mallatdec2(X,w name,level)%输入:X载入的二维图像像数值level小波分解次（级）数设定值（如果设定值超过最高可分解次数,按最高分解次数分解）wname 小波名字 wavelet name%输岀：丫多极小波分解后的小波系数矩阵h,g=wfilters(w name,d);%h,g分别为低通和高通滤波器X=double(X);t=1;hh=size(X,2);while t=level%先进行行小波变换for row=1:hhY(row,1:hh)=mdec1(X(row,1:hh),h,g);end%再进行列小波变换for col=1:hhtemp=mdec1( Y(1:hh,col),h,g);Y(1:hh,col)=temp;endt=t+1;hh=hh/2;X=Y;end实现%内部子函数，对一行(row)矢量进行一次小波变换，利用fftfun cti on y=mdec1(x,h,g)%输入：x行数组%h为低通滤波器%g为高通滤波器%输岀：y进行一级小波分解后的系数len x=size(x,2);len h=size(h,2);rh=h(e nd:-1:1);rrh=zeros(1,(le nx-le nh),rh;rrh=circshift(rrh,1);rg=g(e nd:-1:1);rrg=zeros(1,(le nx-le nh),rg;rrg=circshift(rrg,1);r1=dyaddow n( ifft(fft(x).*fft(rrh,le nx),1);%use para 1r2=dyaddow n(ifft(fft(x).*fft(rrg,le nx),1);y=r1,r2;(2) 二维小波重构函数%二维小波重构函数fun cti on Y=mallatrec2(X,w name,level)%输入:X载入的小波系数矩阵；,按最高分解次%level 小波分解次(级)数设定值(如果设定值超过最高可分解次数数分解)%wname小波名字 wavelet name%输岀：Y重构图像矩阵h,g=wfilters(wname,d);%h,g分别为重构低通滤波器和重构高通滤波器hz=size(X,2);h1=hz/(2A(level-1);while h1(/*(一CD1 +. _ j2 f广;W（ k力当小波函娄编程时先离散化,在通过求和代替积分。了7比（巳）加k怎G2站匸/（血）X0严） 根据彩亍a1%巧金闻（字）本实验采取的一些小技巧：卩1 L J因为实验中的数据通过lgd（ MattLinaf ）先将数据拷到work文箱 下）命令后，数据存在了矩阵中所以兀知就是血;由于数扮样间隔主为003常量人所以可以把这个系数忽略，实际上本程序彳 -1走_屈的是 WAa.b）=a 3 V fk）产）、丿i I 茁 ? I 拓 V Rd）本实验中墨西哥帽函数为：= C和0, otherwiseC3）如果按照（巧式计算小波系数矩阵，则算法时间复杂度为口涉）在 適的电脑上要运行5分钟左右才能出正确结果.通过硏可以技 对于同一尺度与枫H）有如下关系：呦匕口）从 aaa二项开始与财二相同，只是第一项不同，（在程序中有说明人所D a把算法时间复杂度降为巩沪），程序执行时间在$秒钟内口U）程序的正确性可由xvt来验证口亠（?）本实验还编了个朴前程序来说明存在哪些分量】即对正弦波函裳 样后进行连续小波变换，分析图象结果口和宾验结果：“(1)实骑数据图像和小波系数图像(2)怕哉数据图像和小波系数图像心正弦波的小波系数1015202530354020406080 IX 120140163180200程序附录:(1) 墨西哥帽小波函数,按照(*)式编程 % mexh.m fun cti on Y=mexh(x)if abs(x)=M+N-1%按照圆周卷积的定义，相关和线性卷积的关系(原始算法和线性卷积的关系) 冰意画图理解clc;clear;t仁cputime;load(data.mat);len=len gth(dat);lna=70;% a (尺度)的长度a=zeros(1. In a);% a表示尺度wfab=zeros(l na,le n);%小波系数矩阵mexhab=zeros(1,2*le n-1);data=zeros(1,le n-1),dat;Ydata=fft( data ,4*le n);for s=1:l nafor k=1:2*le n-1mexhab(k)=mexh(k-le n) /s);endtemp=ifft( Ydata.*fft( mexhab,4*le n ) ,4*le n);为什么要取实部而wfab(s,:)=real(temp(2*le n-1:3*le n-2)/sqrt(s); %不是取模，我也不是很清楚，可是有种感觉endfigure(1);plot(dat);title(原始数据图);figure(2);%小波系数谱image(wfab);colorm ap(pi nk(128);title(小波系数谱);cputime-tl4)fft快速计算cwt%按照圆周卷积的定义,冰意画图理解clc;clear;t仁cputime;load(data.mat);len=len gth(dat);Ina=70;% a (尺度)的长度a=5;data=dat,zeros(1,le n);Ydata=fft(dat,2*le n);for s=1:l namexhab=zeros(1,2*le n);k=-a*s:1:a*s;mexhab(k+le n)=mexh2(k./s);temp=ifft( Ydata.*fft( mexhab,2*len ) ,2*le n);要取实部而不是取模,wfab(s,:)=real(temp(len+1:2*len)/sqrt(s); % 呵呵endfigure(1);plot(dat);title(原始数据图);figure(2);%小波系数谱image(wfab);coIorm ap(pi nk(128);title( 小波系数谱);cputime-tl5)保存为 mexh2.m fun cti on Y=mexh2(x)Y=exp(-x.*x/2).*(1-x.A2);