Wavelets+and+Filter+Banks2015spring_1

1、Wavelets and Filter Banks小波與濾波器組小波與濾波器組（一）（一）彭思龍S中國(guó)科學(xué)院自動(dòng)化研究所講在課程的前面課程的由來(lái)；課程的特點(diǎn)；課程的要求；聽(tīng)課的方法；課程參考書(shū)；課程的內(nèi)容。課程的由來(lái)小波（wavelet）是信號(hào)處理、圖像處理乃至數(shù)學(xué)領(lǐng)域20多年來(lái)最重要的學(xué)科分支之一；課程為典型的交叉學(xué)科：同時(shí)產(chǎn)生于諸多工程領(lǐng)域和數(shù)學(xué)，又同時(shí)被工程和數(shù)學(xué)所推動(dòng)；地球物理領(lǐng)域的子波分析；數(shù)學(xué)領(lǐng)域的調(diào)和分析；圖像處理中的金字塔壓縮；最終導(dǎo)致效果多尺度分析理論的出現(xiàn)，形成了獨(dú)立的學(xué)科歷來(lái)課程偏重于工程或者數(shù)學(xué)；缺乏交叉2001年來(lái)在研究生院講授該課程。課程的特點(diǎn)理論和應(yīng)用相結(jié)合：一半

2、講理論，一半講應(yīng)用；理論部分：一半講工程角度，一半講數(shù)學(xué)角度。有形和無(wú)形相結(jié)合：盡量描述思想，而不獨(dú)側(cè)重于具體理論推導(dǎo)。詳細(xì)和簡(jiǎn)略相結(jié)合：理論部分詳細(xì)；應(yīng)用部分簡(jiǎn)略。古典和現(xiàn)代相結(jié)合：既有基本的小波理論，也有最新進(jìn)展。課程的特點(diǎn)（續(xù)）靜態(tài)和動(dòng)態(tài)相結(jié)合：每年都會(huì)更新部分講課內(nèi)容，應(yīng)用內(nèi)容全部自選，理論節(jié)選。理科和文科相結(jié)合：我的學(xué)習(xí)和科研的思考，有理科的，也有文科的?！皼](méi)有情趣的人，不會(huì)有真學(xué)問(wèn)沒(méi)有情趣的人，不會(huì)有真學(xué)問(wèn)”梁漱溟梁漱溟“不識(shí)廬山真面目，只緣身在此山中不識(shí)廬山真面目，只緣身在此山中”蘇軾題西林壁課程的要求線性代數(shù)（高等代數(shù)）；多項(xiàng)式理論，矩陣?yán)碚?；?shù)字信號(hào)處理；數(shù)字圖像處理；泛函分

3、析初步；Matlab。注：缺少數(shù)字信號(hào)處理的同學(xué)一定要自己學(xué)習(xí)信號(hào)處理。缺少數(shù)學(xué)的更要認(rèn)真聽(tīng)課。聽(tīng)課的方法工程和數(shù)學(xué)兩種思維：對(duì)于工程背景的同學(xué)：認(rèn)真理解其應(yīng)用思想，努力掌握其數(shù)學(xué)理論；對(duì)于數(shù)學(xué)基礎(chǔ)較好的同學(xué)：努力掌握其物理背景，認(rèn)真理解其數(shù)學(xué)理論；不能掉隊(duì)，認(rèn)真聽(tīng)講，前后有邏輯性；課后認(rèn)真復(fù)習(xí)，動(dòng)手編程序、推導(dǎo)公式?！拔覀兿嘈?，任何一個(gè)人的學(xué)問(wèn)成就，都是出我們相信，任何一個(gè)人的學(xué)問(wèn)成就，都是出于自學(xué)。學(xué)校教育不過(guò)給學(xué)生開(kāi)一個(gè)端，使他于自學(xué)。學(xué)校教育不過(guò)給學(xué)生開(kāi)一個(gè)端，使他更容易自學(xué)而已更容易自學(xué)而已?！绷菏椋何业淖詫W(xué)小史課程參考書(shū)課程參考書(shū)Wavelet and filter banks,

4、 G. Strang, T. Nguyen, Wellesley-Cambridge Press, 1997 （有MIT的ppt 中文有翻譯，瑞士聯(lián)邦工學(xué)院M. Vetterlli的ppt）多抽樣率信號(hào)處理多抽樣率信號(hào)處理，宗孔德，清華大學(xué)出版社，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（信號(hào)處理的小波導(dǎo)引，中文有翻譯）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 （課程目錄）第一章第一章 引言 對(duì)小波和濾波器設(shè)計(jì)的歷史做簡(jiǎn)要回顧和相關(guān)概念第二章第二章 濾波器組 抽取和插值。二通道濾波器組重建條件，調(diào)制矩陣和多相矩陣方法。第三章第三章 正交濾波器組 濾波器組構(gòu)造的柵格方法，正交濾波器的幾種構(gòu)造方法。第四章第四章 多尺度分析 泛函分析基礎(chǔ)，正交多尺度分析、正交小波和正交濾波器的關(guān)系。Contents (cont.)第五章第五章 雙正交小波與濾

7、波器 雙正交多尺度分析、雙正交小波與雙正交濾波器第六章第六章 小波濾波器的提升算法 小波濾波器的提升算法介紹第七章第七章 小波應(yīng)用基礎(chǔ) 數(shù)字圖像處理的線性逆問(wèn)題、小波域圖像特征和統(tǒng)計(jì)特性第八章第八章 數(shù)字圖像處理應(yīng)用一：圖像去噪 小波域收縮算法、基于學(xué)習(xí)的小波域快速去噪算法Contents (cont.)第九章第九章 數(shù)字圖像處理應(yīng)用二：圖像編碼 EZW、SPHIT和JPEG2000簡(jiǎn)介第十章第十章 數(shù)字圖像處理應(yīng)用三：復(fù)原、超分辨率 圖像復(fù)原和超分辨率算法第十一章第十一章 小波濾波器自適應(yīng)選取方法 兩種小波濾波器選取理論第十二章第十二章 幾何小波介紹一：Curvelet Curvelet 小

8、波理論Contents (cont.)第十三章第十三章 幾何小波介紹二：Bandelet和Countourlet Bandelet和countourlet理論第十四章第十四章 圖像的非線性表示 Matching pursuit算法、稀疏表示和Basis pursuit算法介紹第十五章第十五章 后小波時(shí)代：EMD和NSP EMD和NSP等非線性信號(hào)分析方法第十六章第十六章 前沿選講 有關(guān)非線性信號(hào)分析的前沿進(jìn)展介紹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年就寫(xiě)成關(guān)于熱傳導(dǎo)的基本論文熱的傳播，向巴黎科學(xué)院呈交，但經(jīng)拉格朗日、拉普拉斯和勒讓德審閱后被科學(xué)院拒絕，1811年又提交了經(jīng)修改的論文，該文獲科學(xué)院大獎(jiǎng)，卻未正式發(fā)表。傅立葉在論文中推導(dǎo)出著名的熱傳導(dǎo)方程 ，并在求解該方程時(shí)發(fā)現(xiàn)解函數(shù)可以由三角函數(shù)構(gòu)成的級(jí)數(shù)形式表示，從而提出任一函數(shù)都可以展成三角函數(shù)的無(wú)窮級(jí)數(shù)。傅立葉級(jí)數(shù)（即三角級(jí)數(shù)）、傅立葉分析等理論均由此創(chuàng)始。來(lái)自花花公子彩頁(yè)瑞典人基本概念Signal （信號(hào)）（信號(hào)）: x(

10、t) or x(n)注：信號(hào)處理雖然久遠(yuǎn)，但是仍然沒(méi)有成熟，基本概念不明確。Filter（濾波器）（濾波器）: a sequence, h=h(n)Filtering （濾波）（濾波） y=h*x, where * is the convolution operator: FIR=Finite Impulse Response=finite length（有限脈沖響應(yīng)濾波器）IIR=Infinite Impulse Response=infinite length （無(wú)限脈沖響應(yīng)濾波器）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（連續(xù)傅里葉變換連續(xù)傅里葉變換） 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 幅頻響應(yīng)：相頻響應(yīng), =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公式證明的理論都需要借助于廣義函數(shù)理論才能得到嚴(yán)格的證明。對(duì)于數(shù)學(xué)專業(yè)的同學(xué)大多困惑于此。是廣義函數(shù)的典型代表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注：永無(wú)止境注：永無(wú)止境Signal sampling 信號(hào)采樣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.多查閱原始文獻(xiàn)教科書(shū)往往割裂了結(jié)論與問(wèn)題教科書(shū)往往割裂了結(jié)論與問(wèn)題經(jīng)典理論并不是一蹴

18、而就的經(jīng)典理論并不是一蹴而就的教科書(shū)往往是后人不斷整理的結(jié)果教科書(shū)往往是后人不斷整理的結(jié)果看看原始文獻(xiàn)能夠看到大家思考問(wèn)題的看看原始文獻(xiàn)能夠看到大家思考問(wèn)題的邏輯邏輯直接查閱原始文獻(xiàn)，了解其發(fā)展的歷史，再直接查閱原始文獻(xiàn)，了解其發(fā)展的歷史，再現(xiàn)理論的原貌，了解原創(chuàng)者的思想現(xiàn)理論的原貌，了解原創(chuàng)者的思想Shannon sampling theorem(香農(nóng)采樣定理)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注：“好讀書(shū)，不求好讀書(shū)，不求甚解；每有會(huì)意，便甚解；每有會(huì)意，便欣然忘食。欣然忘食。”-陶淵明陶淵明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、頻響應(yīng)） 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 （反對(duì)稱）小練習(xí)，自己證明；()( ) |( )|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.注：物理可解，就一物理可解，就一定有數(shù)學(xué)解定有數(shù)學(xué)解Filter bank （濾波器組）（濾波器組）=Lowpass + Highpass (互補(bǔ)), 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 （延時(shí)、延遲）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（加權(quán)最小二乘） (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