小波实验 (1)

哈尔滨工业大学小波理论及应用上机报告院（系）电气工程及自动化学院学科仪器仪表工程姓名陈鹏宇学号15S101121上机实验内容Butterworth滤波器实验一：滤波器一、实验内容Butterworth 滤波器冲击响应函数为：,0()0,0t Ae t h t t α-⎧≥=⎨<⎩若若 ⑴ 求()ˆhω； ⑵ 判断是否因果；是低通、高通、带通还是带阻？⑶ 对于信号3()(sin 22sin 40.4sin 2sin 40)t f t e t t t t -=++，0t π≤≤，画出()f t 图形；⑷ 画出滤波后图形()f h t *，比较滤波前后图形，你会发现什么，这里取10A α==⑸ 取()(sin5sin3sin sin 40),t f t e t t t t -=+++采用不同的变量值A α=()10A α==初始设定，画出原信号图形与滤波后图形，比较滤波效果 二、实验过程（1） 由连续傅里叶变换公式可求得()ˆhω： ()()0ˆi t t i t A h h t e dt Ae e dt i ωαωωαω+∞+∞----∞===+⎰⎰ (1) （2） 由冲击响应函数 h(t)可知，t<0时h(t)=0;可以得出滤波器为因果滤波器；由()ˆhω可知该信号的幅频特性为： ()H ω= (2) 当α1A ==时，其幅频特性曲线如图1所示，因此为低通滤波器。

图1 滤波器幅频特性（3）()f t，(b)为幅频特性f t的图形如图二所示，(a)为()图2 时域和频域图形（4）画出*()f h t的结果如图2所示，A= α = 10：图3 f*h(t)卷积结果其中图3（a）、（b）分别为原始信号f(t)的时域和频域，（c）、（d）为h t(A =α = 10)卷积后的信号的时域和频域图。

可以看出信号中高频分量被抑制，信号的信噪比明显改善了。

（5）实验中A和的取值，通过实验得到的结果分别如下：1)A= α = 102)A= α = 20，滤波时域和频域图形3)A= α= 100，滤波时域和频域图形4)A= α= 5，滤波时域和频域图形5)A= α= 1，滤波时域和频域图形图8 A= α= 1，滤波时域和频域图形通过对比各种参数滤波器的滤波效果，可以发现A值的大小会影响信号的幅值，A值越大滤波器对信号的放大作用越大，但噪声也被放大，而则影响滤波器的截止频率，越小h(t)的截止频率越小，对高频信号的滤除效果越好。

一、题目：一维Haar 小波2次分解二、目的：编程实现信号的分解与重构三、算法及其实现：离散小波变换离散小波变换是对信号的时－频局部化分析，其定义为：/2200()(,)()(),()()j j Wf j k a f t a t k dt f t L R φ+∞---∞=-∈⎰ 本实验实现对信号的分解与重构：（１）信号分解：用小波工具箱中的dwt 函数来实现离散小波变换，函数dwt 将信号分解为两部分，分别称为逼近系数和细节系数（也称为低频系数和高频系数）,实验中分别记为cA1,cD1，它们的长度均为原始信号的一半，但dwt 只能实现原始信号的单级分解。

在本实验中使用小波函数db1来实现单尺度小波分解，即：[cA1,cD1]＝dwt(s,’db1’)，其中s 是原信号;再通过[cA2,cD2]＝dwt(cA1,’db1’)进行第二次分解，长度又为cA2的一半。

（２）信号重构：用小波工具箱中的upcoef 来实现，upcoef 是进行一维小波分解系数的直接重构,即：A1 = upcoef('a',cA1,'db1'); D1 = upcoef('a',cD1,'db1')。

四、实现工具：Matlab五、程序代码：%装载leleccum 信号load leleccum;s = leleccum(1:3920);%用小波函数db1对信号进行单尺度小波分解[cA1,cD1]=dwt(s,'db1');subplot(3,2,1);plot(s);title('leleccum 原始信号');%单尺度低频系数cA1向上一步的重构信号A1 = upcoef('a',cA1,'db1');%单尺度高频系数cD1向上一步的重构信号D1 = upcoef('a',cD1,'db1');subplot(3,2,3);plot(A1);title('单尺度低频系数cA1向上一步的重构信号');subplot(3,2,5);plot(D1);title('单尺度高频系数cD1向上一步的重构信号');[cA1,cD1]=dwt(cA1,’db1');subplot(3,2,2);plot(s);title('leleccum 第一次分解后的cA1信号');%第二次分解单尺度低频系数cA2向上一步的重构信号A2= upcoef('a',cA2,'db1',2);%第二次分解单尺度高频系数cD2向上一步的重构信号D2 = upcoef('a',cD2,'db1',2);subplot(3,2,4);plot(A2);title('第二次分解单尺度低频系数cA2向上一步的重构信号');subplot(3,2,6);plot(D2);title('的二次分解单尺度高频系数cD2向上一步的重构信号');六、运行结果：七、结果分析：。

小波分析实验报告
目的：比较传统的AR模型直接估计和经过小波去噪后的AR模型估计哪个更好数据：上海证券交易所A股股票从20031112 到2003112131 (共有240 个交易日) 每交易日的收盘价格。

过程：先选取数据，然后用AR模型直接估计得到一组估计值。

再用选取小波对数据进行分解重构，得到去噪后的数据，在对此数据用AR模型估计，
结论：
通过表3 我们可以看到本文提出的小波预测方法，在预测处理金融数据这类非平稳的时间序列时，同传统的预测方法相比较具有一定的可靠性，具有很好的应用前景。

通过表4 我们可以看到,当分解层数在3～4 层时，预测效果比较好。

实际过程中，如果待预测的时间序列数据量不是很大，我们分解层数一般采用3～4 层。

分工：
孟红波查找相关文献和数据
王骏建模，分析数据。

图像的小波降噪实验报告孙玉祥314113002432一．背景在图像处理过程中，图像的采集、转换和传输常常受到成像设备与外部环境噪声干扰等影响，产生降质。

图像噪声对数字图像的后续处理影响较大，因此对图像噪声的去除有很重要的显示意义。

传统的降噪方法多采用平均或线性方法进行，常用的是维纳滤波，但是降噪效果不够好。

随着小波理论的日益完善，它以自身良好的时频特性在图像降噪领域受到越来越多的关注，开辟了用非线性降噪的先河。

二.原理2.1 小波在图像处理方面的优点小波降噪主要是利用噪声与图像信号在频率上分布的不同，图像信号主要分布在低频区域，而噪声主要分布在高频区域。

小波去噪使得原始图像的结构信息和细节信息很容易被提取是因为小波具有以下特点：（1）低熵性。

小波洗漱的稀疏分布，使得图像变换后的熵降低；多分辨率性。

优于采用了多分辨率分析，因此可以非常好地刻画信号的非平稳特征，如边缘、尖峰、断点等；（2）去相关性。

因为小波变换可以对信号进行去相关，且噪声在变换后有白化趋势，所以小波域比空域更利于去噪；（3）选基灵活性。

优于小波变换可以灵活选取变换基，从而对不同的应用场合，不同的研究对象，可以选用不同的小波母函数，以获得最佳的效果。

2.2 小波去噪方法到目前为止，小波去噪的方法大概分为三大类：第一类方法是基于小波变换模极大值原理，根据信号和噪声在小波变换各尺度上的不同传播特性，剔除由噪声产生的模极大值点，保留信号所对应的模极大值点，然后利用所余模极大值点重构小波系数，进而恢复信号；第二类方法是对含噪信号作小波变换之后，计算相邻尺度间小波系数的相关性，根据相关性的大小区别小波系数的模型，从而进行取舍，然后直接重构信号；第三类方法是阈值方法，该方法就是对小波分解后的各层系数模大于和小于某阈值的系数分别进行处理，然后利用处理后的小波系数重构出降噪后的图像。

2.3 小波阈值去噪小波阈值去噪法有着很好的数学理论支持，实现简单而又非常有效，因此取得了非常大的成功，并吸引了众多学者对其作进一步的研究与改进。

小波分析实验报告课程：小波分析姓名：学院：学号：一、实验目的：1、运用傅里叶变换知识对常用的基本函数做基本变换。

2、通过观察小波变换系数建立对小波变换及其有关性质的感性认识3、加深对因果滤波器的理解，并会判断因果滤波器的类型。

4、运用卷积公式对基本信号做滤波处理并作出分析，以加深理解。

5、熟悉Matlab中相关函数的用法。

二、实验原理：1、“小波”就是小区域、长度有限、均值为0的波形。

所谓“小”是指它具有衰减性；而称之为“波”则是指它的波动性，其振幅正负相间的震荡形式。

与Fourier变换相比，小波变换是时间(空间)频率的局部化分析，它通过伸缩平移运算对信号（函数)逐步进行多尺度细化，最终达到高频处时间细分，低频处频率细分，能自动适应时频信号分析的要求，从而可聚焦到信号的任意细节，解决了Fourier变换的困难问题，成为继Fourier变换以来在科学方法上的重大突破。

小波转换分成两个大类：离散小波变换(DWT)和连续小波转换(CWT)。

两者的主要区别在于，连续转换在所有可能的缩放和平移上操作，而离散转换采用所有缩放和平移值的特定子集。

小波变换的公式有内积形式和卷积形式，两种形式的实质都是一样的。

它要求的就是一个个小波分量的系数也就是“权”。

其直观意义就是首先用一个时窗最窄，频窗最宽的小波作为尺子去一步步地“量”信号，也就是去比较信号与小波的相似程度。

信号局部与小波越相似，则小波变换的值越大，否则越小。

当一步比较完成后，再将尺子拉长一倍，又去一步步地比较，从而得出一组组数据。

如此这般循环，最后得出的就是信号的小波分解（小波级数）。

当尺度及位移均作连续变化时，可以理解必将产生大量数据，作实际应用时并不需要这么多的数据，因此就产生了离散的思想。

将尺度作二进离散就得到二进小波变换，同时也将信号的频带作了二进离散。

当觉得二进离散数据量仍显大时，同时将位移也作离散就得到了离散小波变换。

2、二维离散小波变换常用函数三、实验内容：1. 对信号noissin 分别采用图形接口和命令行两种方式进行单尺度小波分解重构和多尺度小波分解重构层数为4，并显示各层低频高频图形，加以比较。

实验四一、实验目的理解小波阈值去噪法原理。

对所得的去噪效果进行分析。

二、实验要求在载入原始图片后，对图片进行含噪和消噪处理，再对所得的图片效果进行分析。

三、主要内容载入原始图片，对原始图片添加一个随机噪声，得出含噪图片。

用sym6小波对图像进行1层分解，设置一个全局阈值，对图像分解系数，将低频系数进行重构，得出消噪后的图像。

再与原图像，含噪图像一起进行分析比较。

运行代码如下clear all;load woman;subplot(2,2,1);image(X);colormap(map);xlabel('(a)原始图像');axis square;init=2055615866;randn('seed',init);x=X+48*randn(size(X));subplot(2,2,2);image(x);colormap(map);xlabel('(b)含噪图像');axis square;%用sym6小波对图像进行1层分解t1=wpdec2(x,1,'sym6');%设置一个全局阈值thr=10.358;%对图像分解系数t2=wpthcoef(t1,0,'s',thr);%对低频系数进行重构x1=wprcoef(t1,1);subplot(2,2,3);image(x1);运行结果四、思考体会小波去噪的根本任务是在小波域将信号的小波变换与噪声的小波变换有效的分离。

噪声的能量分布于整个小波域内，小波分解后，信号的小波系数幅值要大于噪声的系数幅值，也可以认为，幅值比较大的小波系数一般以信号为主，而比较小的系数在很大程度上是噪声。

于是，采用阈值的方法可把信号系数保留，而把大部分噪声系数减少至零。

将含噪信号在各尺度上进行小波分解，保留大尺度（低分辨率）下的全部系数，对于小尺度（高分辨率）下的小波系数，设定一个阈值，幅值不超过阈值的小波系数设置为零，幅值高于该阈值的小波系数或者完整保留，或者做相应的收缩处理，最后将处理后的小波系数利用逆小波变换进行重构，恢复出有效信号。

小波分析实验：实验2 二维离散小波变换(Mallat快速算法)实验目的：在理解离散小波变换原理和Mallat快速算法的基础上，通过编程对图像进行二维离散小波变换，从而加深对二维小波分解和重构的理性和感性认识，并能提高编程能力，为今后的学习和工作奠定基础。

实验工具：计算机，matlab6.5附录：（1）二维小波分解函数%二维小波分解函数function Y=mallatdec2(X,wname,level)%输入:X 载入的二维图像像数值;% level 小波分解次(级)数设定值(如果设定值超过最高可分解次数,按最高分解次数分解)% wname 小波名字wavelet name%输出：Y 多极小波分解后的小波系数矩阵[h,g]=wfilters(wname,'d'); %h,g分别为低通和高通滤波器X=double(X);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)矢量进行一次小波变换，利用fft实现function y=mdec1(x,h,g)%输入：x 行数组% h为低通滤波器% g为高通滤波器%输出： y 进行一级小波分解后的系数lenx=size(x,2);lenh=size(h,2);rh=h(end:-1:1);rrh=[zeros(1,(lenx-lenh)),rh];rrh=circshift(rrh',1)';rg=g(end:-1:1);rrg=[zeros(1,(lenx-lenh)),rg];rrg=circshift(rrg',1)';r1=dyaddown(ifft(fft(x).*fft(rrh,lenx)),1); %use para 1r2=dyaddown(ifft(fft(x).*fft(rrg,lenx)),1);y=[r1,r2];（2）二维小波重构函数%二维小波重构函数function Y=mallatrec2(X,wname,level)%输入:X 载入的小波系数矩阵;% level 小波分解次(级)数设定值(如果设定值超过最高可分解次数,按最高分解次数分解)% wname 小波名字wavelet name%输出：Y 重构图像矩阵[h,g]=wfilters(wname,'d'); %h,g分别为重构低通滤波器和重构高通滤波器hz=size(X,2);h1=hz/(2^(level-1));while h1<=hz% 对列变换for col=1:h1temp=mrec1(X(1:h1,col)',h,g)';X(1:h1,col)=temp;end%再对行变换for row=1:h1temp=mrec1(X(row,1:h1),h,g);X(row,1:h1)=temp;endh1=h1*2;endY=X;%内部子函数，对一行小波系数进行重构function y=mrec1(x,h,g)%输入：x 行数组% h为低通滤波器% g为高通滤波器%输出： y 进行一级小波重构后值lenx=size(x,2);r3=dyadup(x(1,1:lenx*0.5),0); %内插零use para 0r4=dyadup(x(1,(lenx*0.5+1):lenx),0); %use para 0y=ifft(fft(r3,lenx).*fft(h,lenx))+ ifft(fft(r4,lenx).*fft(g,lenx));（3）测试函数（主函数）%测试函数（主函数）clc;clear;X=imread('E:\Libin的文档\Course\Course_wavelet\实验2要求\exp2\LENA.bmp');%路径X=double(X);A = mallatdec2(X,'sym2',3);image(abs(A));colormap(gray(255));title('多尺度分解图像');Y= mallatrec2(A,'sym2',3);Y=real(Y);figure(2);subplot(1,2,1);image(X);colormap(gray(255));title('原始图像');subplot(1,2,2);image(Y);colormap(gray(255));title('重构图像');csize=size(X);sr=csize(1);sc=csize(2);mse=sum(sum( (Y-X).^2,1))/(sr*sc);psnr=10*log(255*255/mse)/log(10)小波分析实验：实验1 连续小波变换实验目的：在理解连续小波变换原理的基础上，通过编程实现对一维信号进行连续小波变换，（实验中采用的是墨西哥帽小波），从而对连续小波变换增加了理性和感性的认识，并能提高编程能力，为今后的学习和工作奠定基础。

一、题目：双树复小波变换二、目的：双树复小波和实小波变换的比较三、算法及其实现：提取阶梯型边界点1．算法。

幅值:相位：2．代码实现。

我采用Matlab函数编程实现。

具体程序见shift_test_2D.m，drawcirc.m，setfig.m，dtwavexfm2.m，dtwaveifm2.m,waveifm2.m,wavexfm2.mSkelMap.m。

设和分别是双正交对偶尺度函数与对偶小波，，，和实部：虚部：双树复小波变换可以通过离散小波变换DWT实现：一个DWT产生实部，另一个产生虚部。

四、实现工具：Matlab五、程序代码：（1）shift_test_2D.m：% shift_test_2D.m%% M-file to perform a 4-level wavelet transform on a circle using Q-shift% dual wavelet tree and DWT, and to compare shift invariance properties.%% Nick Kingsbury, Cambridge University, May 2002.clear allclose all% Draw a circular disc.x = round((drawcirc(64,1,0,0,256) - 0.5) * 200);setfig(1);colormap(gray(256))image(min(max(x+128,1),256));set(gca,'position',[0.1 0.25 .25 .5]);()ctψ=,,j,c h gf fψψψ=+|(,)|cd j n=(,)arctan(,)gchd j nd j n∠= ⎪⎝⎭,h hφφ%,h hψψ%()h n()h n%()1h n()1h n%11()()(2)()()(2)()()(2)()()(2)h hnh hnh hnh hnt h n t nt h n t nt h n t nt h n t nφφψφφφψφ=-=-=-=-%%%%%%11()()(2)()()(2)()()(2)()()(2)g gng gng gng gnt g n t nt g n t nt g n t nt g n t nφφψφφφψφ=-=-=-=-%%%%%%axis('off');axis('image');% draw(xx);title('Input (256 x 256)','FontSize',14); drawnow% Do 4 levels of CWT.[Yl,Yh] = dtwavexfm2(x,4,'near_sym_b','qshift_b');% Loop to reconstruct output from coefs at each level in turn.% Starts with the finest level.titl = ['1st';'2nd';'3rd';'4th';'Low'];yy = zeros(size(x) .* [2 3]);yt1 = 1:size(x,1); yt2 = 1:size(x,2);for mlev = 1:5,mask = zeros(6,5);mask(:,mlev) = 1;z = dtwaveifm2(Yl*mask(1,5),Yh,'near_sym_b','qshift_b',mask);figure;draw(z);drawnowyy(yt1,yt2) = z;yt2 = yt2 + size(x,2)/2;end% disp('Press a key ...')% pause% Now do same with DWT.% Do 4 levels of Real DWT using 'antonini' (9,7)-tap filters. [Yl,Yh] = wavexfm2(x,4,'antonini');yt1 = [1:size(x,1)] + size(x,1); yt2 = 1:size(x,2);for mlev = 1:5,mask = zeros(3,5);mask(:,mlev) = 1;z = waveifm2(Yl*mask(1,5),Yh,'antonini',mask);figure;draw(z);drawnowyy(yt1,yt2) = z;yt2 = yt2 + size(x,2)/2;endfigure;setfig(gcf);colormap(gray(256))image(min(max(yy+128,1),256));set(gca,'position',[0.1 0.1 .8 .8]);axis('off');axis('image');hold onplot(128*[[1;1]*[1:4] [0;6]]+1,128*[[0;4]*[1 1 1 1] [2;2]]+1,'-k');hold offtitle('Components of reconstructed ''disc'' images','FontSize',14);text(-0.01*size(yy,2),0.25*size(yy,1),'DT CWT','horiz','r');text(0.02*size(yy,2),1.02*size(yy,1),'wavelets:','horiz','r','vert','t');text(-0.01*size(yy,2),0.75*size(yy,1),'DWT','horiz','r');for k=1:4, text(k*128-63,size(yy,1)*1.02,sprintf('level %d',k),'FontSize',14,'horiz','c','vert','t'); endtext(5*128+1,size(yy,1)*1.02,'level 4 scaling fn.','FontSize',14,'horiz','c','vert','t');drawnow% print -deps circrecq.epsdisp('Press a key to see perfect reconstruction property ...')pause% Accumulate the images from lowband upwards to show perfect reconstruction.sy = size(x,2)/2;for mlev = 4:-1:1,yt2 = [1:sy] + (mlev-1)*sy;yy(:,yt2) = yy(:,yt2) + yy(:,yt2+sy);endfigure;setfig(gcf);colormap(gray(256))image(min(max(yy+128,1),256));set(gca,'position',[0.1 0.1 .8 .8]);axis('off');axis('image');title('Accumulated reconstructions from each level of DT CWT ','FontSize',14);text(size(yy,2)*0.5,size(yy,1)*1.02,'Accumulated reconstructions from each level of DWT ','FontSize',14,'hor','c','vert','t');drawnowreturn（2）function p = drawcirc(r,w,dx,dy,N)% function p = drawcirc(r,w,dx,dy,N)% Generate an image of size N*N pels, containing a circle% radius r pels and centred at dx,dy relative% to the centre of the image. The edge of the circle is a cosine shaped% edge of width w (from 10 to 90% points).x = ones(N,1) * (([1:N] - (N+1)/2 - dx)/r);y = (([1:N]' - (N+1)/2 - dy)/r) * ones(1,N);p = 0.5 + 0.5 * sin(min(max((exp(-0.5 * (x + y )) - exp(-0.5))*(r*3/w), -pi/2), pi/2));return（3）function Z = dtwaveifm2(Yl,Yh,biort,qshift,gain_mask);% Function to perform an n-level dual-tree complex wavelet (DTCWT)% 2-D reconstruction.%% Z = dtwaveifm2(Yl,Yh,biort,qshift,gain_mask);%% Yl -> The real lowpass image from the final level% Yh -> A cell array containing the 6 complex highpass subimages for each level.%% biort -> 'antonini' => Antonini 9,7 tap filters.% 'legall' => LeGall 5,3 tap filters.% 'near_sym_a' => Near-Symmetric 5,7 tap filters.% 'near_sym_b' => Near-Symmetric 13,19 tap filters.%% qshift -> 'qshift_06' => Quarter Sample Shift Orthogonal (Q-Shift) 10,10 tap filters,% (only 6,6 non-zero taps).% 'qshift_a' => Q-shift 10,10 tap filters,% (with 10,10 non-zero taps, unlike qshift_06).% 'qshift_b' => Q-Shift 14,14 tap filters.% 'qshift_c' => Q-Shift 16,16 tap filters.% 'qshift_d' => Q-Shift 18,18 tap filters.%% gain_mask -> Gain to be applied to each subband.% gain_mask(d,l) is gain for subband with direction d at level l.% If gain_mask(d,l) == 0, no computation is performed for band (d,l).% Default gain_mask = ones(6,length(Yh)).%% Z -> Reconstructed real image matrix%%% For example: Z = dtwaveifm2(Yl,Yh,'near_sym_b','qshift_b');% performs a 3-level reconstruction from Yl,Yh using the 13,19-tap filters% for level 1 and the Q-shift 14-tap filters for levels >= 2.%% Nick Kingsbury and Cian Shaffrey% Cambridge University, May 2002a = length(Yh); % No of levels.if nargin < 5, gain_mask = ones(6,a); end % Default gain_mask.if isstr(biort) & isstr(qshift) %Check if the inputs are stringsbiort_exist = exist([biort '.mat']);qshift_exist = exist([qshift '.mat']);if biort_exist == 2 & qshift_exist == 2; %Check to see if the inputs exist as .mat files load (biort);load (qshift);elseerror('Please enter the correct names of the Biorthogonal or Q-Shift Filters, see help DTW A VEIFM2 for details.');endelseerror('Please enter the names of the Biorthogonal or Q-Shift Filters as shown in help DTW A VEIFM2.');endcurrent_level = a;Z = Yl;while current_level >= 2; ; %this ensures that for level -1 we never do the following lh = c2q(Yh{current_level}(:,:,[1 6]),gain_mask([1 6],current_level));hl = c2q(Yh{current_level}(:,:,[3 4]),gain_mask([3 4],current_level));hh = c2q(Yh{current_level}(:,:,[2 5]),gain_mask([2 5],current_level));% Do even Qshift filters on columns.y1 = colifilt(Z,g0b,g0a) + colifilt(lh,g1b,g1a);y2 = colifilt(hl,g0b,g0a) + colifilt(hh,g1b,g1a);% Do even Qshift filters on rows.Z = (colifilt(y1.',g0b,g0a) + colifilt(y2.',g1b,g1a)).';% Check size of Z and crop as required[row_size col_size] = size(Z);S = 2*size(Yh{current_level-1});if row_size ~= S(1) %check to see if this result needs to be cropped for the rows Z = Z(2:row_size-1,:);endif col_size ~= S(2) %check to see if this result needs to be cropped for the cols Z = Z(:,2:col_size-1);endif any(size(Z) ~= S(1:2)),error('Sizes of subbands are not valid for DTW A VEIFM2');endcurrent_level = current_level - 1;endif current_level == 1;lh = c2q(Yh{current_level}(:,:,[1 6]),gain_mask([1 6],current_level));hl = c2q(Yh{current_level}(:,:,[3 4]),gain_mask([3 4],current_level));hh = c2q(Yh{current_level}(:,:,[2 5]),gain_mask([2 5],current_level));% Do odd top-level filters on columns.y1 = colfilter(Z,g0o) + colfilter(lh,g1o);y2 = colfilter(hl,g0o) + colfilter(hh,g1o);% Do odd top-level filters on rows.Z = (colfilter(y1.',g0o) + colfilter(y2.',g1o)).';endreturn%==================================================================== ======================% ********** INTERNAL FUNCTION **********%==================================================================== ======================function x = c2q(w,gain)% function z = c2q(w,gain)% Scale by gain and convert from complex w(:,:,1:2) to real quad-numbers in z. %% Arrange pixels from the real and imag parts of the 2 subbands% into 4 separate subimages .% A----B Re Im of w(:,:,1)% | |% | |% C----D Re Im of w(:,:,2)sw = size(w);x = zeros(2*sw(1:2));if any(w(:)) & any(gain)sc = sqrt(0.5) * gain;P = w(:,:,1)*sc(1) + w(:,:,2)*sc(2);Q = w(:,:,1)*sc(1) - w(:,:,2)*sc(2);t1 = 1:2:size(x,1);t2 = 1:2:size(x,2);% Recover each of the 4 corners of the quads.x(t1,t2) = real(P); % a = (A+C)*sc;x(t1,t2+1) = imag(P); % b = (B+D)*sc;x(t1+1,t2) = imag(Q); % c = (B-D)*sc;x(t1+1,t2+1) = -real(Q); % d = (C-A)*sc;endreturn（3）function [Yl,Yh,Yscale] = dtwavexfm2(X,nlevels,biort,qshift);% Function to perform a n-level DTCWT-2D decompostion on a 2D matrix X %% [Yl,Yh,Yscale] = dtwavexfm2(X,nlevels,biort,qshift);%% X -> 2D real matrix/Image%% nlevels -> No. of levels of wavelet decomposition%% biort -> 'antonini' => Antonini 9,7 tap filters.% 'legall' => LeGall 5,3 tap filters.% 'near_sym_a' => Near-Symmetric 5,7 tap filters.% 'near_sym_b' => Near-Symmetric 13,19 tap filters.%% qshift -> 'qshift_06' => Quarter Sample Shift Orthogonal (Q-Shift) 10,10 tap filters,% (only 6,6 non-zero taps).% 'qshift_a' => Q-shift 10,10 tap filters,% (with 10,10 non-zero taps, unlike qshift_06).% 'qshift_b' => Q-Shift 14,14 tap filters.% 'qshift_c' => Q-Shift 16,16 tap filters.% 'qshift_d' => Q-Shift 18,18 tap filters.%%% Yl -> The real lowpass image from the final level% Yh -> A cell array containing the 6 complex highpass subimages for each level.% Yscale -> This is an OPTIONAL output argument, that is a cell array containing% real lowpass coefficients for every scale.%%% Example: [Yl,Yh] = dtwavexfm2(X,3,'near_sym_b','qshift_b');% performs a 3-level transform on the real image X using the 13,19-tap filters% for level 1 and the Q-shift 14-tap filters for levels >= 2.%% Nick Kingsbury and Cian Shaffrey% Cambridge University, Sept 2001if isstr(biort) & isstr(qshift) %Check if the inputs are stringsbiort_exist = exist([biort '.mat']);qshift_exist = exist([qshift '.mat']);if biort_exist == 2 & qshift_exist == 2; %Check to see if the inputs exist as .mat filesload (biort);load (qshift);elseerror('Please enter the correct names of the Biorthogonal or Q-Shift Filters, see help DTW A VEXFM2 for details.');endelseerror('Please enter the names of the Biorthogonal or Q-Shift Filters as shown in help DTW A VEXFM2.');endorginal_size = size(X);if ndims(X) >= 3;error(sprintf('The entered image is %dx%dx%d, please enter each image sliceseparately.',orginal_size(1),orginal_size(2),orginal_size(3)));end% The next few lines of code check to see if the image is odd in size, if so an extra ...% row/column will be added to the bottom/right of the imageinitial_row_extend = 0; %initialiseinitial_col_extend = 0;if any(rem(orginal_size(1),2)), %if sx(1) is not divisable by 2 then we need to extend X by adding a row at the bottomX = [X; X(end,:)]; %Any further extension will be done in due course.initial_row_extend = 1;endif any(rem(orginal_size(2),2)), %if sx(2) is not divisable by 2 then we need to extend X by adding a col to the leftX = [X X(:,end)]; %Any further extension will be done in due course.initial_col_extend = 1;endextended_size = size(X);if nlevels == 0, return; end%initialiseYh=cell(nlevels,1);if nargout == 3Yscale=cell(nlevels,1); %this is only required if the user specifies a third output component.endS = [];sx = size(X);if nlevels >= 1,% Do odd top-level filters on cols.Lo = colfilter(X,h0o).';Hi = colfilter(X,h1o).';% Do odd top-level filters on rows.LoLo = colfilter(Lo,h0o).'; % LoLoYh{1} = zeros([size(LoLo)/2 6]);Yh{1}(:,:,[1 6]) = q2c(colfilter(Hi,h0o).'); % Horizontal pairYh{1}(:,:,[3 4]) = q2c(colfilter(Lo,h1o).'); % Vertical pairYh{1}(:,:,[2 5]) = q2c(colfilter(Hi,h1o).'); % Diagonal pairS = [ size(LoLo) ;S];if nargout == 3Yscale{1} = LoLo;endendif nlevels >= 2;for level = 2:nlevels;[row_size col_size] = size(LoLo);if any(rem(row_size,4)), % Extend by 2 rows if no. of rows of LoLo are divisable by 4;LoLo = [LoLo(1,:); LoLo; LoLo(end,:)];endif any(rem(col_size,4)), % Extend by 2 cols if no. of cols of LoLo are divisable by 4;LoLo = [LoLo(:,1) LoLo LoLo(:,end)];end% Do even Qshift filters on rows.Lo = coldfilt(LoLo,h0b,h0a).';Hi = coldfilt(LoLo,h1b,h1a).';% Do even Qshift filters on columns.LoLo = coldfilt(Lo,h0b,h0a).'; %LoLoYh{level} = zeros([size(LoLo)/2 6]);Yh{level}(:,:,[1 6]) = q2c(coldfilt(Hi,h0b,h0a).'); % HorizontalYh{level}(:,:,[3 4]) = q2c(coldfilt(Lo,h1b,h1a).'); % VerticalYh{level}(:,:,[2 5]) = q2c(coldfilt(Hi,h1b,h1a).'); % DiagonalS = [ size(LoLo) ;S];if nargout == 3Yscale{level} = LoLo;endendendYl = LoLo;if initial_row_extend == 1 & initial_col_extend == 1;warning(sprintf(' \r\r The image entered is now a %dx%d NOT a %dx%d \r The bottom row and rightmost column have been duplicated, prior to decomposition. \r\r ',...extended_size(1),extended_size(2),orginal_size(1),orginal_size(2)));endif initial_row_extend == 1 ;warning(sprintf(' \r\r The image entered is now a %dx%d NOT a %dx%d \r Row number %d has been duplicated, and added to the bottom of the image, prior to decomposition. \r\r',...extended_size(1),extended_size(2),orginal_size(1),orginal_size(2),orginal_size(1)));endif initial_col_extend == 1;warning(sprintf(' \r\r The image entered is now a %dx%d NOT a %dx%d \r Col number %d has been duplicated, and added to the right of the image, prior to decomposition. \r\r',...extended_size(1),extended_size(2),orginal_size(1),orginal_size(2),orginal_size(2)));endreturn%==================================================================== ======================% ********** INTERNAL FUNCTION **********%==================================================================== ======================function z = q2c(y)% function z = q2c(y)% Convert from quads in y to complex numbers in z.sy = size(y);t1 = 1:2:sy(1); t2 = 1:2:sy(2);j2 = sqrt([0.5 -0.5]);% Arrange pixels from the corners of the quads into% 2 subimages of alternate real and imag pixels.% a----b% | |% | |% c----d% Combine (a,b) and (d,c) to form two complex subimages.p = y(t1,t2)*j2(1) + y(t1,t2+1)*j2(2); % p = (a + jb) / sqrt(2)q = y(t1+1,t2+1)*j2(1) - y(t1+1,t2)*j2(2); % q = (d - jc) / sqrt(2)% Form the 2 subbands in z.z = cat(3,p-q,p+q);return（4）function Z = dtwaveifm2(Yl,Yh,biort,qshift,gain_mask);% Function to perform an n-level dual-tree complex wavelet (DTCWT)% 2-D reconstruction.%% Z = dtwaveifm2(Yl,Yh,biort,qshift,gain_mask);%% Yl -> The real lowpass image from the final level% Yh -> A cell array containing the 6 complex highpass subimages for each level.%% biort -> 'antonini' => Antonini 9,7 tap filters.% 'legall' => LeGall 5,3 tap filters.% 'near_sym_a' => Near-Symmetric 5,7 tap filters.% 'near_sym_b' => Near-Symmetric 13,19 tap filters.%% qshift -> 'qshift_06' => Quarter Sample Shift Orthogonal (Q-Shift) 10,10 tap filters, % (only 6,6 non-zero taps).% 'qshift_a' => Q-shift 10,10 tap filters,% (with 10,10 non-zero taps, unlike qshift_06).% 'qshift_b' => Q-Shift 14,14 tap filters.% 'qshift_c' => Q-Shift 16,16 tap filters.% 'qshift_d' => Q-Shift 18,18 tap filters.%% gain_mask -> Gain to be applied to each subband.% gain_mask(d,l) is gain for subband with direction d at level l.% If gain_mask(d,l) == 0, no computation is performed for band (d,l). % Default gain_mask = ones(6,length(Yh)).%% Z -> Reconstructed real image matrix%%% For example: Z = dtwaveifm2(Yl,Yh,'near_sym_b','qshift_b');% performs a 3-level reconstruction from Yl,Yh using the 13,19-tap filters% for level 1 and the Q-shift 14-tap filters for levels >= 2.%% Nick Kingsbury and Cian Shaffrey% Cambridge University, May 2002a = length(Yh); % No of levels.if nargin < 5, gain_mask = ones(6,a); end % Default gain_mask.if isstr(biort) & isstr(qshift) %Check if the inputs are stringsbiort_exist = exist([biort '.mat']);qshift_exist = exist([qshift '.mat']);if biort_exist == 2 & qshift_exist == 2; %Check to see if the inputs exist as .mat files load (biort);load (qshift);elseerror('Please enter the correct names of the Biorthogonal or Q-Shift Filters, see help DTW A VEIFM2 for details.');endelseerror('Please enter the names of the Biorthogonal or Q-Shift Filters as shown in help DTW A VEIFM2.');endcurrent_level = a;Z = Yl;while current_level >= 2; ; %this ensures that for level -1 we never do the following lh = c2q(Yh{current_level}(:,:,[1 6]),gain_mask([1 6],current_level));hl = c2q(Yh{current_level}(:,:,[3 4]),gain_mask([3 4],current_level));hh = c2q(Yh{current_level}(:,:,[2 5]),gain_mask([2 5],current_level));% Do even Qshift filters on columns.y1 = colifilt(Z,g0b,g0a) + colifilt(lh,g1b,g1a);y2 = colifilt(hl,g0b,g0a) + colifilt(hh,g1b,g1a);% Do even Qshift filters on rows.Z = (colifilt(y1.',g0b,g0a) + colifilt(y2.',g1b,g1a)).';% Check size of Z and crop as required[row_size col_size] = size(Z);S = 2*size(Yh{current_level-1});if row_size ~= S(1) %check to see if this result needs to be cropped for the rows Z = Z(2:row_size-1,:);endif col_size ~= S(2) %check to see if this result needs to be cropped for the cols Z = Z(:,2:col_size-1);endif any(size(Z) ~= S(1:2)),error('Sizes of subbands are not valid for DTW A VEIFM2');endcurrent_level = current_level - 1;endif current_level == 1;lh = c2q(Yh{current_level}(:,:,[1 6]),gain_mask([1 6],current_level));hl = c2q(Yh{current_level}(:,:,[3 4]),gain_mask([3 4],current_level));hh = c2q(Yh{current_level}(:,:,[2 5]),gain_mask([2 5],current_level));% Do odd top-level filters on columns.y1 = colfilter(Z,g0o) + colfilter(lh,g1o);y2 = colfilter(hl,g0o) + colfilter(hh,g1o);% Do odd top-level filters on rows.Z = (colfilter(y1.',g0o) + colfilter(y2.',g1o)).';endreturn%==================================================================== ======================% ********** INTERNAL FUNCTION **********%==================================================================== ======================function x = c2q(w,gain)% function z = c2q(w,gain)% Scale by gain and convert from complex w(:,:,1:2) to real quad-numbers in z.%% Arrange pixels from the real and imag parts of the 2 subbands% into 4 separate subimages .% A----B Re Im of w(:,:,1)% | |% | |% C----D Re Im of w(:,:,2)sw = size(w);x = zeros(2*sw(1:2));if any(w(:)) & any(gain)sc = sqrt(0.5) * gain;P = w(:,:,1)*sc(1) + w(:,:,2)*sc(2);Q = w(:,:,1)*sc(1) - w(:,:,2)*sc(2);t1 = 1:2:size(x,1);t2 = 1:2:size(x,2);% Recover each of the 4 corners of the quads.x(t1,t2) = real(P); % a = (A+C)*sc;x(t1,t2+1) = imag(P); % b = (B+D)*sc;x(t1+1,t2) = imag(Q); % c = (B-D)*sc;x(t1+1,t2+1) = -real(Q); % d = (C-A)*sc;endreturn(5)function [Yl,Yh,Yscale] = wavexfm2(X,nlevels,biort);% Function to perform a n-level DWT-2D decompostion on a 2-D matrix X.%% [Yl,Yh,Yscale] = dtwavexfm2(X,nlevels,biort);%% X -> real 1-D signal column vector (or matrix of vectors)%% nlevels -> No. of levels of wavelet decomposition%% biort -> 'antonini' => Antonini 9,7 tap filters.% 'legall' => LeGall 5,3 tap filters.% 'near_sym_a' => Near-Symmetric 5,7 tap filters.% 'near_sym_b' => Near-Symmetric 13,19 tap filters.%% Yl -> The lowpass subband from the final level.% Yh -> A cell array containing the highpass subband for each level.% Yscale -> This is an OPTIONAL output argument, that is a cell array containing% the lowpass coefficients at every scale.%%% Example: [Yl,Yh] = wavexfm2(X,4,'near_sym_b');% performs a 4-level 2-D DWT on the real image X using the 13,19-tap filters.%% Nick Kingsbury, Cambridge University, May 2002if isstr(biort) % Check if the biort input is a stringbiort_exist = exist([biort '.mat']);if biort_exist == 2, % Check to see if the filter exists as a .mat fileload (biort);elseerror('Please enter the correct name of the Biorthogonal Filter, see help W A VEXFM2 for details.');endelseerror('Please enter the name of the Biorthogonal Filter as shown in help W A VEXFM2.');endL = size(X);if any(rem(L,2)), % ensure that X is an even length, thus enabling it to be extended if needs be.error('Size of X must be a multiple of 2');end%initialiseYh=cell(nlevels,1);if nargout == 3Yscale=cell(nlevels,1); % This is only required if the user specifies a third output component.endLoLo = X;for level = 1:nlevels;if rem(size(LoLo,1),4), % Check to see if height of LoLo is divisable by 4, if not extend.LoLo = [LoLo(1,:); LoLo; LoLo(end,:)];endif rem(size(LoLo,2),4), % Check to see if height of LoLo is divisable by 4, if not extend.LoLo = [LoLo(:,1) LoLo LoLo(:,end)];end% Do filters on rows.Lo = coldwtfilt(LoLo,h0o,0).';Hi = coldwtfilt(LoLo,h1o,1).';% Do filters on columns.LoLo = coldwtfilt(Lo,h0o,0).'; %LoLoYh{level} = zeros([size(LoLo) 3]);Yh{level}(:,:,1) = coldwtfilt(Hi,h0o,0).'; % HorizontalYh{level}(:,:,3) = coldwtfilt(Lo,h1o,1).'; % VerticalYh{level}(:,:,2) = coldwtfilt(Hi,h1o,1).'; % Diagonalif nargout == 3Yscale{level} = LoLo;endendYl = LoLo;return(6)function Z = waveifm2(Yl,Yh,biort,gain_mask);% Function to perform an n-level DWT 2-D reconstruction.%% Z = waveifm2(Yl,Yh,biort,gain_mask);%% Yl -> The real lowpass subband from the final level% Yh -> A cell array containing the complex highpass subband for each level.%% biort -> 'antonini' => Antonini 9,7 tap filters.% 'legall' => LeGall 5,3 tap filters.% 'near_sym_a' => Near-Symmetric 5,7 tap filters.% 'near_sym_b' => Near-Symmetric 13,19 tap filters.%% gain_mask -> Gain to be applied to each subband.% gain_mask(l) is gain for wavelet subband at level l.% If gain_mask(l) == 0, no computation is performed for band (l).% Default gain_mask = ones(1,length(Yh)).%% Z -> Reconstructed real signal vector (or matrix).%%% For example: Z = waveifm2(Yl,Yh,'near_sym_b');% performs a reconstruction from Yl,Yh using the 13,19-tap filters.%% Nick Kingsbury, Cambridge University, May 2002nlevels = length(Yh); % No of levels.if nargin < 4, gain_mask = ones(3,nlevels); end % Default gain_mask.if isstr(biort) % Check if the biort input is a stringbiort_exist = exist([biort '.mat']);if biort_exist == 2, % Check to see if the filter exists as a .mat fileload (biort);elseerror('Please enter the correct name of the Biorthogonal Filter, see help W A VEIFM2 for details.');endelseerror('Please enter the name of the Biorthogonal Filter as shown in help W A VEIFM2.');endLoLo = Yl;for level = nlevels:-1:1, % Reconstruct levels in reverse order.if size(LoLo,1) ~= size(Yh{level},1) % If LoLo is not the same height as the next Yh => t1 was extended.LoLo = LoLo(2:size(LoLo,1)-1,:); % Therefore we have to clip LoLo so it is the same height as the next Yh.endif size(LoLo,2) ~= size(Yh{level},2) % If LoLo is not the same width as the next Yh => t1was extended.LoLo = LoLo(:,2:size(LoLo,2)-1); % Therefore we have to clip LoLo so it is the same width as the next Yh.endif any([size(LoLo) 3] ~= size(Yh{level})),error('Yh sizes are not valid for WA VEIFM2');endlh = Yh{level}(:,:,1) * gain_mask(1,level);hl = Yh{level}(:,:,3) * gain_mask(3,level);hh = Yh{level}(:,:,2) * gain_mask(2,level);% Do even Qshift filters on columns.y1 = coliwtfilt(LoLo,2*g0o,0) + coliwtfilt(lh,2*g1o,1);y2 = coliwtfilt(hl,2*g0o,0) + coliwtfilt(hh,2*g1o,1);% Do even Qshift filters on rows.LoLo = (coliwtfilt(y1.',2*g0o,0) + coliwtfilt(y2.',2*g1o,1)).';endZ = LoLo;return(7)function setfig(f)% function setfig(f)% Set figure dimensions for correct display of text to agree with% print -deps.figure(f)set(gcf,'position',[220 40 790 526],'DefaultTextFontSize',14);return六、运行结果：七、结果分析：。

图像拼接实验报告一、实验目的选用适当的拼接算法实现两幅图像的拼接。

二、实验原理图像拼接技术就是将数张有重叠部分的图像（可能是不同时间、不同视角或者不同传感器获得的）拼成一幅大型的无缝高分辨率图像的技术。

图像配准和图像融合是图像拼接的两个关键技术。

图像配准是图像融合的基础,而且图像配准算法的计算量一般非常大,因此图像拼接技术的发展很大程度上取决于图像配准技术的创新。

早期的图像配准技术主要采用点匹配法,这类方法速度慢、精度低,而且常常需要人工选取初始匹配点,无法适应大数据量图像的融合。

图像拼接的方法很多,不同的算法步骤会有一定差异,但大致的过程是相同的。

一般来说,图像拼接主要包括以下五步:（1）图像预处理：包括数字图像处理的基本操作(如去噪、边缘提取、直方图处理等)、建立图像的匹配模板以及对图像进行某种变换(如傅里叶变换、小波变换等)等操作。

（2）图像配准：就是采用一定的匹配策略,找出待拼接图像中的模板或特征点在参考图像中对应的位置,进而确定两幅图像之间的变换关系。

（3）建立变换模型：根据模板或者图像特征之间的对应关系,计算出数学模型中的各参数值,从而建立两幅图像的数学变换模型。

（4）统一坐标变换：根据建立的数学转换模型,将待拼接图像转换到参考图像的坐标系中,完成统一坐标变换。

（5）融合重构：将带拼接图像的重合区域进行融合得到拼接重构的平滑无缝全景图像。

图像拼接技术主要包括两个关键环节即图像配准和图像融合对于图像融合部分，由于其耗时不太大，且现有的几种主要方法效果差别也不多，所以总体来说算法上比较成熟。

而图像配准部分是整个图像拼接技术的核心部分，它直接关系到图像拼接算法的成功率和运行速度，因此配准算法的研究是多年来研究的重点。

目前的图像配准算法基本上可以分为两类:基于频域的方法(相位相关方法)和基于时域的方法。

相位相关法对拼接的图像进行快速傅立叶变换，将两幅待配准图像变换到频域，然后通过它们的互功率谱直接计算出两幅图像间的平移矢量，从而实现图像的配准。

一、题目：一维Haar 小波１次分解二、目的：编程实现信号的分解与重构三、算法及其实现：离散小波变换离散小波变换是对信号的时－频局部化分析，其定义为：/2200()(,)()(),()()j j Wf j k a f t a t k dt f t L R φ+∞---∞=-∈⎰ 本实验实现对信号的分解与重构：（１）信号分解：用小波工具箱中的dwt 函数来实现离散小波变换，函数dwt 将信号分解为两部分，分别称为逼近系数和细节系数（也称为低频系数和高频系数）,实验中分别记为cA1,cD1，它们的长度均为原始信号的一半，但dwt 只能实现原始信号的单级分解。

在本实验中使用小波函数db1来实现单尺度小波分解，即：[cA1,cD1]＝dwt(s,’db1’)，其中s 是原信号;（２）信号重构：用小波工具箱中的upcoef 来实现，upcoef 是进行一维小波分解系数的直接重构,即： A1 = upcoef('a',cA1,'db1'); D1 = upcoef('a',cD1,'db1')。

四、实现工具：Matlab五、程序代码：%装载leleccum 信号load leleccum;s = leleccum(1000:3920);%用小波函数db1对信号进行单尺度小波分解[cA1,cD1]=dwt(s,'db1');subplot(3,2,1);plot(s);title('leleccum 原始信号');subplot(3,2,2);plot(s);title('leleccum 原始信号');subplot(3,2,4);plot(cA1);title('单尺度低频系数cA1');subplot(3,2,6);plot(cD1);title('单尺度高频系数cD1');%单尺度低频系数cA1向上一步的重构信号cA1 = upcoef('a',cA1,'db1');%单尺度高频系数cD1向上一步的重构信号cD1 = upcoef('a',cD1,'db1');subplot(3,2,3);plot(cA1);title('单尺度低频系数cA1向上一步的重构信号');subplot(3,2,5);plot(cD1);title('单尺度高频系数cD1向上一步的重构信号');六、运行结果：七、结果分析：。