Wavelets+and+Filter+Banks2015spring_1

1、Wavelets and Filter Banks小波与滤波器组小波与滤波器组（一）（一）彭思龙S中国科学院自动化研究所讲在课程的前面课程的由来；课程的特点；课程的要求；听课的方法；课程参考书；课程的内容。课程的由来小波（wavelet）是信号处理、图像处理乃至数学领域20多年来最重要的学科分支之一；课程为典型的交叉学科：同时产生于诸多工程领域和数学，又同时被工程和数学所推动；地球物理领域的子波分析；数学领域的调和分析；图像处理中的金字塔压缩；最终导致效果多尺度分析理论的出现，形成了独立的学科历来课程偏重于工程或者数学；缺乏交叉2001年来在研究生院讲授该课程。课程的特点理论和应用相结合：一半

2、讲理论，一半讲应用；理论部分：一半讲工程角度，一半讲数学角度。有形和无形相结合：尽量描述思想，而不独侧重于具体理论推导。详细和简略相结合：理论部分详细；应用部分简略。古典和现代相结合：既有基本的小波理论，也有最新进展。课程的特点（续）静态和动态相结合：每年都会更新部分讲课内容，应用内容全部自选，理论节选。理科和文科相结合：我的学习和科研的思考，有理科的，也有文科的。“没有情趣的人，不会有真学问没有情趣的人，不会有真学问”梁漱溟梁漱溟“不识庐山真面目，只缘身在此山中不识庐山真面目，只缘身在此山中”苏轼题西林壁课程的要求线性代数（高等代数）；多项式理论，矩阵理论；数字信号处理；数字图像处理；泛函分

3、析初步；Matlab。注：缺少数字信号处理的同学一定要自己学习信号处理。缺少数学的更要认真听课。听课的方法工程和数学两种思维：对于工程背景的同学：认真理解其应用思想，努力掌握其数学理论；对于数学基础较好的同学：努力掌握其物理背景，认真理解其数学理论；不能掉队，认真听讲，前后有逻辑性；课后认真复习，动手编程序、推导公式。“我们相信，任何一个人的学问成就，都是出我们相信，任何一个人的学问成就，都是出于自学。学校教育不过给学生开一个端，使他于自学。学校教育不过给学生开一个端，使他更容易自学而已更容易自学而已。”梁漱溟：我的自学小史课程参考书课程参考书Wavelet and filter banks,

4、 G. Strang, T. Nguyen, Wellesley-Cambridge Press, 1997 （有MIT的ppt 中文有翻译，瑞士联邦工学院M. Vetterlli的ppt）多抽样率信号处理多抽样率信号处理，宗孔德，清华大学出版社，1996。Multirate systems and filter banks, Vaidyanathan, PP., Englewood Cliffs, New Jersey, Prentice Hall Inc. 1993.A wavelet tour of signal processing, S. Mallat, Academic Press

5、. NY, 1998（信号处理的小波导引，中文有翻译）Ten Lectures on Wavelets, Ingrid Daubechies, 1992（小波十讲，有中文翻译）S课程名称解释Filter banks（滤波器组）=a set of filters, filter is widely used in many fields of engineering and science for a long time.Wavelet （小波）, an old and new tool to produce filter banks, have been thoroughly studied

6、in past 25 years. Here we use wavelets to indicate many kinds of wavelets with different properties.Contents （课程目录）第一章第一章 引言 对小波和滤波器设计的历史做简要回顾和相关概念第二章第二章 滤波器组 抽取和插值。二通道滤波器组重建条件，调制矩阵和多相矩阵方法。第三章第三章 正交滤波器组 滤波器组构造的栅格方法，正交滤波器的几种构造方法。第四章第四章 多尺度分析 泛函分析基础，正交多尺度分析、正交小波和正交滤波器的关系。Contents (cont.)第五章第五章 双正交小波与滤

7、波器 双正交多尺度分析、双正交小波与双正交滤波器第六章第六章 小波滤波器的提升算法 小波滤波器的提升算法介绍第七章第七章 小波应用基础 数字图像处理的线性逆问题、小波域图像特征和统计特性第八章第八章 数字图像处理应用一：图像去噪 小波域收缩算法、基于学习的小波域快速去噪算法Contents (cont.)第九章第九章 数字图像处理应用二：图像编码 EZW、SPHIT和JPEG2000简介第十章第十章 数字图像处理应用三：复原、超分辨率 图像复原和超分辨率算法第十一章第十一章 小波滤波器自适应选取方法 两种小波滤波器选取理论第十二章第十二章 几何小波介绍一：Curvelet Curvelet 小

8、波理论Contents (cont.)第十三章第十三章 几何小波介绍二：Bandelet和Countourlet Bandelet和countourlet理论第十四章第十四章 图像的非线性表示 Matching pursuit算法、稀疏表示和Basis pursuit算法介绍第十五章第十五章 后小波时代：EMD和NSP EMD和NSP等非线性信号分析方法第十六章第十六章 前沿选讲 有关非线性信号分析的前沿进展介绍Contents (cont.)Some ideas in life and research （杂谈）博客：晴朗的天空博客：晴朗的天空-彭思龙彭思龙http:/ bit of his

9、tory:from Fourier to Haar to wavelets(from the material of M. Vetterlli)傅立叶早在1807年就写成关于热传导的基本论文热的传播，向巴黎科学院呈交，但经拉格朗日、拉普拉斯和勒让德审阅后被科学院拒绝，1811年又提交了经修改的论文，该文获科学院大奖，却未正式发表。傅立叶在论文中推导出著名的热传导方程 ，并在求解该方程时发现解函数可以由三角函数构成的级数形式表示，从而提出任一函数都可以展成三角函数的无穷级数。傅立叶级数（即三角级数）、傅立叶分析等理论均由此创始。来自花花公子彩页瑞典人基本概念Signal （信号）（信号）: x(

10、t) or x(n)注：信号处理虽然久远，但是仍然没有成熟，基本概念不明确。Filter（滤波器）（滤波器）: a sequence, h=h(n)Filtering （滤波）（滤波） y=h*x, where * is the convolution operator: FIR=Finite Impulse Response=finite length（有限脉冲响应滤波器）IIR=Infinite Impulse Response=infinite length （无限脉冲响应滤波器）kknxkhny)()()(020406080100120-1.5-1-0.500.511.50204060

11、80100120-1.5-1-0.500.511.5Example of filtering:x=sin(t)+0.1*randn(1,101);h=1 1 1 1/4; y=x*hContinuous Fourier Transform（连续傅里叶变换连续傅里叶变换） Some basic properties:Linearity （线性性）Parseval Identity: （Parseval 等式）22( )( ) |( )|Fourier Transform:1( )( ) or ( ) 2Inverse Fourier Transform:11( )( ) or ( ) 22Ri

12、ti tRRititRRf tL Rf tdtff t edtf t edtf tfedfed 1,; 1 or 2where , Rf gcf gcf gf gZ transform （Z-变换）and DTFTGiven a signal or filter ( ), Z transform is defined as: ( )( )Discrete Time Fourier Transform (DTFT)( )( ); -1( ) is also called frequency response of ( ).( ) | ( )|njnjs nS zs n zSs n ejSs nS

13、Se(), , | ( )|: , ( )*( )( )( )( )( )( )|( )| |( )|*|( )|( )( )( )yxhSyx hY zX z H zYXHYXH 幅频响应：相频响应, =0低频， = 最高频Poisson Summation Formula（Poisson公式）We use Dirac function as a sampling function:Poisson summation formula:Equivalent form:0( )0other( ) ()( ), ( )1ttx tta dtx at dt特别的ZkZnitnekt21)2(knnG

14、kG)(21)2(注：凡是涉及到用Poisson公式证明的理论都需要借助于广义函数理论才能得到严格的证明。对于数学专业的同学大多困惑于此。是广义函数的典型代表Problem:Signal class：x(t)How to choose Classical sampling theory:For band-limited and energy limited Limited: time finite signal is not band-limited.Wavelet sampling theory: Energy limited New theory: Finite rate of innov

15、ation注：永无止境注：永无止境Signal sampling 信号采样012( )( ), ( ), ( ),x tx tx tx t ( )( )( )iix tx tt and ( )iittSampling theory historyBY H. NYQUIST, Certain Topics in Telegraph Transmission Theory, American Telephone and Telegraph, February 13-17, 1928.CLAUDE E. SHANNON, Communication in the Presence of Noise,

16、 1949Theorem 1: If a function contains no frequencies higher than w cps, it is completely determined by giving its ordinates at a series of points spaced 1/2w seconds apart.Other materialsJ. M. Whittaker, Interpolatory Function Theory, Cambridge Tracts in Mathematics and Mathematical Physics, no. 33

17、. Cambridge, U.K.: Cambridge Univ. Press, ch. IV, 1935.W. R. Bennett, “Time division multiplex systems,” Bell Syst. Tech. J., vol. 20, p. 199, Apr. 1941, D. Gabor, “Theory of communication,” J. Inst. Elect. Eng. (London),vol. 93, pt. 3, no. 26, p. 429, 1946.多查阅原始文献教科书往往割裂了结论与问题教科书往往割裂了结论与问题经典理论并不是一蹴

18、而就的经典理论并不是一蹴而就的教科书往往是后人不断整理的结果教科书往往是后人不断整理的结果看看原始文献能够看到大家思考问题的看看原始文献能够看到大家思考问题的逻辑逻辑直接查阅原始文献，了解其发展的历史，再直接查阅原始文献，了解其发展的历史，再现理论的原貌，了解原创者的思想现理论的原貌，了解原创者的思想Shannon sampling theorem(香农采样定理)If signal f(t) satisfies: supp( ) is included in the interval -T, T, and sampling rate is r, thenIf r , we can not re

19、construct signal f(t);If r= , we haveWhere is called Nyquist sampling rate.TT)()()(krxkrfrtf. 1| /r, /r,-) supp( satisfieshich function wany is T -T,Tx x)sin(T :lyPurticularf注：“好读书，不求好读书，不求甚解；每有会意，便甚解；每有会意，便欣然忘食。欣然忘食。”-陶渊明陶渊明itikrikriksupp , , Tlet r=/,1therefore (f(t)f()e d)21f(kr)f()e d21 f()e d 2

20、1f(kr) f()ed (*)2ranother point of vRRfTTr ,iew of f()() :r f(), ,rjkkkc e 周 期 化Simple Proof:ik-,-1f()ed =() 2r()(), for , that is()(), for ,For any function with |1 and |0()()() (kikZikrZT TRikrZcrfkrfrfkr erfrfkr efrfkr eRT , ,)Inverse transform will be( )()()ZTT TRf trfkrtkr 0( )0 or 1,( )0|( )|

21、|( )|*|( )|HcHYXHh=h(n), Examples:Simplest: h=1 1/2, differenceDual spline: -1 2 1/4;General:h(n) sum of h is 0.bandpass Filter（带通滤波器）：（带通滤波器）：passing band between All-pass filter（全通滤波器）（全通滤波器）：passing all band.11( )0,( ) 0 or c=1,|( )| |( )|*|( )|zzH zH zcYXH0 and 020406080100120-1.5-1-0.500.511.50

22、20406080100120-1.5-1-0.500.511.5Example of low pass filter:x=sin(t)+0.1*randn(1,101);h=1 1 1 1/4; y=x*hxy020406080100120-1.5-1-0.500.511.5020406080100120-0.4-Example of high pass filter:x=sin(t)+0.1*randn(1,101);h=-1 2 1/4; y=x*hxy-3-2-1012301Magnitude response（幅

23、频响应） of 1 1/2 |( )|HMagnitude response of 1 2 1/4-3-2-1012301|( )|HMagnitude response of -1 2 -1/4-3-2-1012301|( )|HPhase（相位），（相位），Magnitude（幅度）（幅度） is called the magnitude response of H. is called the phase of HIf , we say H has linear phase （线性

24、相位）（线性相位）H has linear phase is equivalent to say H is symmetric or antisymmetric （反对称）小练习，自己证明；()( ) |( )|jHHe)(ba)(|() |H/2/2/222filters with linear phase(1)/ 2cos(/ 2) (1)/ 2sin(/ 2) (12)/ 4(1cos( )/ 2(12)/ 4(1cos( )/ 2jjjjjjjjjjeeeeeeeeeePhase estimation for real signalFor real signal, how to est

25、imate its phase?( ) is real, we must find a analytical signal ( ) such that( )( )( )theoretically, ( )( ) where ( ) is the Hilbert transform of s.Example.( )(2 cos( )*sin(10* )the amplitude is (2 coss tc tc ts tid td tH sH ss ttt( ), and the phase is 10* , but for descrete data,s(n), we use the hilb

26、ert transform. (see testphase.m)tt%example of phase and amplitude of real signalt=0:0.01:4*pi;s=(2+cos(t).*sin(10*t);figure;plot(s); c = hilbert(s);ac = abs(c);figure;plot(ac);pc =unwrap( angle(c);figure;plot(pc) Invertibility （可逆性）（可逆性）Y(z)=H(z)X(z), If H does not equal to 0 at any |z|=1, we say H

27、is invertible, that is to say, we can reconstruct X by: X(z)=Y(z)/H(z), which is a inverse filtering, the filter is 1/H(z). But in most cases, H equals to 0 at some points, we can not reconstruct X exactly.Example: H(z)=1+0.5z is invertible, but H(z)=(1+z)/2 is not invertible. Classical useful inver

28、tible filter includes: all-pass filter, discrete Hilbert transform.Inverse filteringY(z)=H(z)X(z)Problem: if H(z) is not invertible, how to get Y?Ill-posed problem: most modern signal or image processing problems are ill-posed.Mathematically, can not be solved. Physically, we must and can solve them

29、.Another path: filter bank.注：物理可解，就一物理可解，就一定有数学解定有数学解Filter bank （滤波器组）（滤波器组）=Lowpass + Highpass (互补), Simplest idea: H0 and H1, where H0 is a lowpass filter, and H1 is a highpass filter, the lost information in the process of lowpass filtering can be fund in the output of the highpass filter.Some p

30、roblems: How to reconstruct the signal?How to find such filter bank?How to reduce the computation?How to reduce the cost (storage and/or hardware)? Any more properties beside reconstruction?Delay （延时、延迟）Sx(n)=x(n-1)Advance（提前）S-1x(n)=x(n+1)SS-1=S-1S=I, where I represent the unit operator.Time-invari

31、ant filters: H is a linear filter, if H(Sx)=S(Hx): a shift of the input produces a shift of the output.Ideal filters, （理想滤波器）Ideal lowpass:Ideal Highpass:For ideal lowpass fitler: But in practice, we must use finite filter for computation, the first idea is to use finite part of this filter, it lead

32、s to Gibbs phenomenon。ZkikekhH|/2 , 0/2|0 1, )()(ZkikekhH|/2 , 1/2|0 0, )()(sin2( )kh kk-11k11-3-2-10123-0.60.81-21k21-3-2-10123-0.60.81-41k41-3-2-10123-0.60.81-81k81-3-2-10123-0.60.81-201k201-3-2-10123-0.60.81-1.62-1.6-1.58-1.56-1.54-1.52-1.5-1.48-1

33、.460.920.940.960.9811.021.041.061.081.1Overshot: 0.089490.-Gibbs, J. Willard, Fourier Series. Nature 59, 200 (1898) and 606 (1899). Traditional filter design methods:Firstly, we note that we only need to construct lowpass filter, and shift in phase by we can get corresponding highpass filters.lowpas

34、s: (0)1,()0( )()highpass: (0)( )0;( )(2 )(0)1HHHHHHHHHWindow method（窗口法）: h(n)=hI(n)w(n)Hamming window: Hanning window:Kaiser window:1)/2-(N|n| )2cos()1 ()(Nnnw1)/2-(N|n| )2cos(2/12/1)(Nnnw1220020!)5 . 0(1)(IwhereN/2|n| )(/2121)(kkxxINnInwEquiripple method（等纹波方法）:The filter with the smallest maximum

35、 error in passband and stopband is an equiripple filter. Means the ripples in passband and stopband is equal height.Remez exchange algorithmWeighted least squares（加权最小二乘） (eigenfilters（特征滤波器）)This method is to minimize the function:response.frequency desired is )D( whereweight)d(| )H(e - )D(| 2jE纹波减纹波减弱，或弱，或者控制者控制在特定在特定形式下，形式下，但永在。但永在。Heisenbergs Uncertainty Principle:Define two window width:time and frequency:Then: if |f|=1, we haveIf f is the Gaussian function, the minim