使用MATLAB设计小波变换程序中的若干问题(转)

在使用MATLAB完成小波变换程序和通过阈值来压缩图像的过程中,我和许多同学都是边学边用,是从一个接一个的问题中逐步理解小波和MATLAB编写程序的。因此我愿意就个人遇到和解决问题的经验与大家讨论,希望能够对遇到同样问题的人有所帮助。在清华大学林福宗老师倡导的网上互动的学习方式中,老师同学的开诚布公的讨论,尤其是林老师启发大家对出现‘问题’采取的态度和做法,对我今后成长为一名合格的清华的研究生意义重大,谨以此文表示对他指导关心的敬意!  
  
内容简介 
  
本文分为三部分: 如何使用MATLAB设计小波标准与标准分解;如何完成使用小波变换压缩图像;仍需探讨的问题。每一部分主要以问题和例子的形式讲述,为了便于参照,附录部分给出部分源代码供大家参考指正。 
  
我个人在完成作业的时候,走了许多弯路,最后,反复读老师的第3章讲义的第3--5节,明白了小波变换的一些非常基础性的知识,正是如此,我主要是按照第三章的例子,做出使用MATLAB进行Haar变换的实例,这至少对完成非标准Haar小波的3级分解与合成没有任何问题,而且计算速度远远比使用一维变换快,甚至超过dwt2和wavedec2。再者,是通过这些简短的例子的学习,也能从实践的角度理解小波变换的基础知识和道理,我觉得掌握小波的知识要比简单使用dwt2直接分析出结果更重要。为了让使用一维变换或打算采用卷积完成分解重构程序的同学有所参考,也给出了我以前写的代码和设想供有兴趣者参考。在图像压缩的任务中,我个人建议采用wdencmp,原因是简单可靠,能够生成老师要求小波压缩与重构的演示及PNG文件。在最后部分着重对使用小波压缩之后重构图像的PNG文件会变大问题进行了简单的猜测性的解释和讨论,希望感兴趣的同学能深入研究。 
  
如何使用MATLAB设计小波标准与标准分解  
  
使用Haar和Db9小波编写图像文件的标准和非标准分解重构程序  
实质上,任务书中规定首先要编写4种程序之一:哈尔(haar)小波的标准分解重构程序;Daubechies 9小波的标准分解重构程序;哈尔(haar)小波的非标准分解重构程序;Daubechies 9小波的非标准分解重构程序。如果采用标准方法,需要生成图3-25的系列图像。如果采用非标准方法则需要生成图3-26的系列图像。 
  
我个人编写了非标准的分解重构程序,并总结网上同学们公开的方法,认为可以有3种方法实现。但前提是必须仔细阅读林老师推荐的补充教材第三章“小波与小波变换”的3.3、3.4和3.5节(14-27页)。如果能使用MATLAB简单地实践一下教材的例子,至少完成哈尔小波的标准与非标准程序相当简单。为此我主要介绍如何实现老师教材给我们的例子是如何在MATLAB中实现的。然后,再引申到使用MATLAB的函数实现方法。为了描述问题简便,使用黑体表示需在MATLAB的命令窗口(command window)的输入部分。如果不熟悉下面出现的函数功能和使用方法,请再命令窗口中使用 help 函数名,如help dwt2,能够得到英文的功能和用法说明。 
  
哈尔小波变换的MATLAB实例  
先介绍MATLAB的矩阵(Matrix)和向量(Vector)的赋值方法:在command window的>>符号下输入: 
  
一维向量: I = [ 9 7 3 5] 
  
      F = [2, 5, 8, 9, 7, 4, -1, 1] 
  
二维矩阵: A = [  
  
64 2 3 61 60 6 7 57  
  
9 55 54 12 13 51 50 16 
  
17 47 46 20 21 43 42 24 
  
40 26 27 37 36 30 31 33 
  
32 34 35 29 28 38 39 25 
  
41 23 22 44 45 19 18 48 
  
49 15 14 52 53 11 10 56 
  
8 58 59 5 4 62 63 1 
  

  
  
  
I和F现在就是[例3.1][例3.2]一维向量, A就是3.5.1中的图像矩阵。在进行Haar小波的一维与二维变换(分解)之前,首先介绍一下MATLAB的矩阵运算。 
  
MATLAB是Matrix Laboratory 的合写,它对矩阵运算之功能堪称一流。由于使用矩阵描述问题更象数学表达式,所以编写的程序不仅高效,更易读。所以MATLAB程序应该尽量使用矩阵直接描述。 如 M是一4X4的矩阵,则B=I*M就完成了使用两个for循环编写乘法的程序。 
  
分析老师的[例3.1],那么可以得到如下结果: 
  
[(9+7)*1/2 (3+5)* 1/2 (9-7)*1/2 (3-5)* 1/2] 
  
如果把其看作矩阵方式的乘法,那么令 M为其表述求取平均和差值的系数矩阵,则输入: 
  
M = [ 
  
1/2 0 1/2 0 
  
1/2 0 -1/2 0  
  
0 1/2 0 1/2 
  
0 1/2 0 -1/2 
  

  
把上述的M输入到MATLAB中,然后使用: 
  
C=I*M 
  
得到的结果就是 [ 8 4 1 -1]。这就是非规范化Haar小波的第一级分解系数。注意C的前两项[ 8 4 ] 就是近似系数(Approximation Coefficients), 后两项[1 _1]就是细节系数(Detail Coefficients). M就是Haar小波的非规范化(non-normalization)系数矩阵。如果我们执行 : 
  
M*M_ 
  
(M_代表的是M的转置矩阵),你会发现得到对角线为0.5其它为0的4X4矩阵,所以如果让它变成单位矩阵(对角线为1),必须把原来的1/2增大,变为1/sqrt(2), 所以执行: 
  
N=M*sqrt(2) 
  
后得到新的系数矩阵 
  
[ 0.7071 0 0.7071 0 
  
0.7071 0 -0.7071 0 
  
0 0.7071 0 0.7071 
  
0 0.7071 0 -0.7071 
  

  
再执行 
  
N * N_ 
  
就会得到单位矩阵。N就是Haar小波的规范化(normalization)系数矩阵。现在执行: 
  
C*M_ 
  
将得到 [4.5 3.5 1.5 2.5],显然是因为M*M_的对角矩阵的值为0.5。现在使用  
  
Cn=I*N 
  
得到 [11.3137 5.6569 1.4142 -1.4142 ] 规范化的一维Haar小波一级分解系数。 
  
I1=Cn*N_ 
  
得到 [9 3 7 7] 规范化的一维Haar小波逆变换(重构)结果。可见得到这样规范化矩阵的好处是为了逆变换时不需要额外的运算,只是需要分解系数矩阵乘以规范化系数矩阵的转置。 
  
。 
  
如果需要对分解系数矩阵进行第2级分解,那么此时只对C或Cn的前2项的低频系数进行,此时的Haar的非规范化与规范化系数矩阵M1和N1分别为: 
  
M1=[ 
  
1/2 1/2 
  
1/2
  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
f1=50; % 频率1 f2=100; % 频率2 fs=2*(f1+f2); % 采样频率 Ts=1/fs; % 采样间隔 N=120; % 采样点数 n=1:N; y=sin(2*pi*f1*n*Ts)+sin(2*pi*f2*n*Ts); % 正弦混合 figure(1) plot(y); title('两个正弦信号') figure(2) stem(abs(fft(y))); title('两信号频谱') %% 2.小器谱分析 h=wfilters('db30','l'); % 低通 g=wfilters('db30','h'); % 高通 h=[h,zeros(1,N-length(h))]; % 补零(圆周卷积,且增大分辨率变于观察) g=[g,zeros(1,N-length(g))]; % 补零(圆周卷积,且增大分辨率变于观察) figure(3); stem(abs(fft(h))); title('低通滤器图'); figure(4); stem(abs(fft(g))); title('高通滤器图') %% 3.MALLET分解算法(圆周卷积的快速傅里叶变换实现) sig1=ifft(fft(y).*fft(h)); % 低通(低频分量) sig2=ifft(fft(y).*fft(g)); % 高通(高频分量) figure(5); % 信号图 subplot(2,1,1) plot(real(sig1)); title('分解信号1') subplot(2,1,2) plot(real(sig2)); title('分解信号2') figure(6); % 频谱图 subplot(2,1,1) stem(abs(fft(sig1))); title('分解信号1频谱') subplot(2,1,2) stem(abs(fft(sig2))); title('分解信号2频谱') %% 4.MALLET重构算法 sig1=dyaddown(sig1); % 2抽取 sig2=dyaddown(sig2); % 2抽取 sig1=dyadup(sig1); % 2插值 sig2=dyadup(sig2); % 2插值 sig1=sig1(1,[1:N]); % 去掉最后一个零 sig2=sig2(1,[1:N]); % 去掉最后一个零 hr=h(end:-1:1); % 重构低通 gr=g(end:-1:1); % 重构高通 hr=circshift(hr',1)'; % 位置调整圆周右移一位 gr=circshift(gr',1)'; % 位置调整圆周右移一位 sig1=ifft(fft(hr).*fft(sig1)); % 低频 sig2=ifft(fft(gr).*fft(sig2)); % 高频 sig=sig1+sig2; % 源信号 %% 5.比较 figure(7); subplot(2,1,1) plot(real(sig1)); title('重构低频信号'); subplot(2,1,2) plot(real(sig2)); title('重构高频信号'); figure(8); subplot(2,1,1) stem(abs(fft(sig1))); title('重构低频信号频谱'); subplot(2,1,2) stem(abs(fft(sig2))); title('重构高频信号频谱'); figure(9) plot(real(sig),'r','linewidth',2); hold on; plot(y); legend('重构信号','原始信号') title('重构信号与原始信号比较') f1=50; % 频率1 f2=100; % 频率2 fs=2*(f1+f2); % 采样频率 Ts=1/fs; % 采样间隔 N=120; % 采样点数 n=1:N; y=sin(2*pi*f1*n*Ts)+sin(2*pi*f2*n*Ts); % 正弦混合 figure(1) plot(y); title('两个正弦信号') figure(2) stem(abs(fft(y))); title('两信号频谱') %% 2.小器谱分析 h=wfilters('db30','l'); % 低通 g=wfilters('db30','h'); % 高通 h=[h,zeros(1,N-length(h))]; % 补零(圆周卷积,且增大分辨率变于观察) g=[g,zeros(1,N-length(g))]; % 补零(圆周卷积,且增大分辨率变于观察) figure(3); stem(abs(fft(h))); title('低通滤器图'); figure(4); stem(abs(fft(g))); title('高通滤器图') %% 3.MALLET分解算法(圆周卷积的快速傅里叶变换实现) sig1=ifft(fft(y).*fft(h)); % 低通(低频分量) sig2=ifft(fft(y).*fft(g)); % 高通(高频分量) figure(5); % 信号图 subplot(2,1,1) plot(real(sig1)); title('分解信号1') subplot(2,1,2) plot(real(sig2)); title('分解信号2') figure(6); % 频谱图 subplot(2,1,1) stem(abs(fft(sig1))); title('分解信号1频谱') subplot(2,1,2) stem(abs(fft(sig2))); title('分解信号2频谱') %% 4.MALLET重构算法 sig1=dyaddown(sig1); % 2抽取 sig2=dyaddown(sig2); % 2抽取 sig1=dyadup(sig1); % 2插值 sig2=dyadup(sig2); % 2插值 sig1=sig1(1,[1:N]); % 去掉最后一个零 sig2=sig2(1,[1:N]); % 去掉最后一个零 hr=h(end:-1:1); % 重构低通 gr=g(end:-1:1); % 重构高通 hr=circshift(hr',1)'; % 位置调整圆周右移一位 gr=circshift(gr',1)'; % 位置调整圆周右移一位 sig1=ifft(fft(hr).*fft(sig1)); % 低频 sig2=ifft(fft(gr).*fft(sig2)); % 高频 sig=sig1+sig2; % 源信号 %% 5.比较 figure(7); subplot(2,1,1) plot(real(sig1)); title('重构低频信号'); subplot(2,1,2) plot(real(sig2)); title('重构高频信号'); figure(8); subplot(2,1,1) stem(abs(fft(sig1))); title('重构低频信号频谱'); subplot(2,1,2) stem(abs(fft(sig2))); title('重构高频信号频谱'); figure(9) plot(real(sig),'r','linewidth',2); hold on; plot(y); legend('重构信号','原始信号') title('重构信号与原始信号比较')
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值