Wavelets Applications in Signal and Image Processing Motivation! The Fourier Transform X( f ) = ¥ ò -¥ x(t)e-2ip ft dt Problem The FT of stationary and non stationary signals with the same frequency components are equivalent. X( f ) » Y (F) X( f ) = ¥ ò x(t)e -2ip ft dt -¥ Y( f ) = ¥ ò y(t)e -2ip ft dt -¥ i.e. we are lacking time localization Although FT tells us what frequencies appear in the signal it does not tell us at what time they appear! What has caused this? X( f ) = ¥ ò x(t)e-2ip ft dt -¥ e-2iπf is a function of infinite support / infinite window function Short Term Fourier Transform: STFT STFTx(w) (t ', f ) = ò [x(t)w* (t - t ')]e- j 2 p ft dt t w(t) = e (-a( t2 )) 2 Multiple FT over smaller windows translated in time Compactly supported We now have a timefrequency representation YOU CAN ALL GO HOME NOT! Recall: In the time domain we know exactly the value of the signal at any time (time resolution) In the frequency domain we know exactly the frequencies in the signal (frequency resolution) In STFT the kernel is compact … thus we can only see a band of frequencies based on the size of the kernel Consequence Window size is application specific Narrow window -> good time resolution, poor frequency resolution Wide window -> good frequency resolution, poor time resolution Increasing window width Wavelets to the rescue We would like to develop a method independent of the windowing function that gives us a) Good time resolution and poor frequency resolution at high frequencies b) Good frequency resolution but poor time resolution at low frequencies Low frequency => Signal information High frequency => Excess detail or noise in the signal Continuous Wavelet Transform: CWT yY = ò ¥ -¥ x(t)Y* (t)dt Ψ is the mother wavelet, the shape or choice of this depend on the properties of the signal we wish to analyze Time localization Inspect the signal at different time steps Introduce a translation parameter, t’, that controls the translation of the function: y = Y t' ò ¥ -¥ x(t)Y* (t - t ')dt Frequency Localization Inspect the signal for different frequencies Introduce a dilation parameter, s, that controls the scale of the function: *æt -t'ö x(t) Y ç ÷ dt ò -¥ è s ø 1 æt -t'ö Y s,t ' (t) = Yç ÷ s è s ø y Y s,t ' 1 = s y = Y s,t ' ò ¥ -¥ ¥ x(t)Y *s,t ' (t)dt Result Changing dilation parameter: Frequency Localization Changing translation parameter: Time Localization y = Y s,t ' ò ¥ -¥ x(t)Y*s,t ' (t)dt Result +ve response Low response -ve response 0 response Orthogonality / Orthonormal Orthogonal: b b V (x),W (x) = ò V (x).W (x)dx = åV(x)W * (x) = 0 * a a i.e. 2 functions are, at no place the same or, are symmetric Orthonormal: fi (x), f j (x) = di, j So dilations and translations of a wavelet must be orthonormal to itself so as not to influence the construction of the coefficients These allow for perfect reconstructions of the form yk = x(t), Y (x) = k ò x(t)Y (t)dt x(t) = å x(t), Y k (x) Y(t) k * k Inverse Wavelet Transform: ICWT x(t) = 1 Cg ò -¥ ò 0 ys,t 'Y s,t ' (t) ¥ ¥ dsdt ' s2 Denoise by zeroing out coefficients Frequency to time resolution Low scales / high frequencies have good time resolution but poor frequency resolution. High scales / low frequencies have good frequency resolution but poor time resolution. STFT has constant time to frequency resolution as window size is fixed Discrete Wavelet Transform: DWT The Discrete Wavelet Transform is a sampled version of the Continuous case with discrete dilation and translation parameters Filters or different cut of frequencies are used to analyze the signal at different scales or resolutions We will be requiring a scaling filter/function and a wavelet filter/function in this case Scaling function – low pass filter - approximation Wavelet – high pass filter - details Discrete wavelet Ψ Recall that the CW is defined as: Y s,t ' (t) = 1 æt -t'ö Yç ÷ s è s ø In a continuous transform we find the inner product over all scales S and translates t’. However now we must sample s and t’. Logarithmic sampling of s means we need to move in discrete steps on t’ proportional to the scale s. æ t - nt0' s0m ö -m/2 -m ' Y m,n (t) = Yç = s Y(s t nt ÷ 0 0 0) m m ø s0 è s0 1 Dyadic scaling Dyadic scaling, choose s0=2 and t0’=1 Y m,n (t) = 2-m/2 Y(2-m t - n) Later this will lead to a nice down sampling routine DWT to obtain detail coefficients becomes: Y m,n (t) = 2-m/2 Y(2-m t - n) ym,n = ò ¥ -¥ x(t)Y m,n (t)dt Dyadic scaling Discrete scaling function Φ In the CWT we calculated the set of coefficients ψ over all scales s and translations t’ on the continuous signal x(t) As we are sampling x(t) we cannot have these infinite coefficients. We need some way of keeping track of what the wavelet coefficients don’t express. Therefore we must define how we sample the signal based on the current dilation, m, of the wavelet. This is done via a Scaling function We can convolve the signal with the scaling function to get approximation coefficients fm,n (t) = 2 -m/2 f (2-m t - n) j m,n = ò ¥ -¥ x(t)fm,n (t)dt Discrete scaling function Φ Approximation and detail Approximation coefficients, ϕ, are produced by applying the scaling function to the sampled signal. They express the signal at a lower resolution as if the high frequencies had been removed Detail coefficients, ψ, are produced by applying the wavelet to the sampled signal. They express the higher frequency components in the signal. Thus a signal is represented as the sum of approximation and detail coefficients: x(t) = ¥ åj n=-¥ f m0 ,n m0 ,n (t) + m0 ¥ å åy m=-¥ n=-¥ m,n Y m,n (t) = xm0 (t) + m0 åd m=-¥ m (t) Multi-Resolution Analysis, MRA Haar example f (t) = f (2t) + f (2t -1) ì 1 0 < t <1/ 2 ï Y(t) = í -1 1 / 2 < t < 1 ï 0 else î DWT via Filtering Filter convolution : x[n]* h[n] = ¥ å x[k]h[n - k] k=-¥ H (equivalent to wavelet) is high pass, stripping the signal of its lower band frequencies thus its coefficients represent high frequency components G (equivalent to scaling function) is a low pass, stripping the signal of its higher frequencies thus is passed on to the next scale to remove the next band of high frequency DWT via Lifting Filters can be transformed in the time or frequency domain into distinct in-place processing steps on the signal rather than costly convolutions Expressing a wavelet in terms of lifting steps is know as a Second Generation Wavelet Here rather than low and high pass filters we perform a Prediction step and an Update step Prediction – high pass filter – we predict what the signal is Update – low pass filter – based on the prediction we update the signal Lifting Haar Lifting example Take signal x(t) and split it into odd and even pairs As a prediction step take the odd away from the even: dj-1= oddj-P(evenj) As an update step take the mean value of the odd and even parts Sj-1=evenj+U(dj-1) d j-1[n] = s j [2n +1]- s j [2n] d j-1[n] = 2d j-1[n] 1 s j-1[n] = s j [2n]+ d j-1[n] 2 s j-1[n] = 2s j-1[n] s j-1[n] = d j-1[n] = 1 d j-1[n] 2 1 s j-1[n] 2 1 s j [2n] = s j-1[n]- d j-1[n] 2 s j [2n +1] = s j [2n]+ d j-1[n] 2D DWT Wavelets and scaling functions are orthogonal … hence separable. We can apply the transform in one direction then the other Z-transform Fourier Series: X(e jw ) = å x[n]e- jnw nÎZ Z-Transform: X(z) = å x[n]z-n nÎZ Convolution W(z) = X(z)Y (z) = å x[n]y[n]z-2n Shift Left Xleft (z) = zX(z) Shift Right Xright (z) = z-1 X(z) Down sample Xdown2 (z) = Up sample Xup2 (z) = X(z 2 ) nÎZ 1 X(z1/2 ) + X(z -1/2 )) ( 2 Lifting to Polyphase Split: X(z) = X0 (z 2 )+ z-1 X1 (z 2 ) Prediction: Update: P : X1 (z) -T(z)X0 (z) U : X0 (z)+ S(z)X1 (z) ùé 1 S (z) N úê úûêë 0 1 é H (z) H (z) ù 00 01 ú H (z) = ê êë H10 (z) H11 (z) úû é Y (z) ù é X (z) ù ê 0 ú = H (z)ê 0 ú êë Y1 (z) úû êë X1 (z) úû é k 0 H (z) = ê êë 0 k -1 ùé 1 0 ù é 1 S0 (z) úê ú... ê úûêë -TN (z) 1 úû êë 0 1 ùé 1 0 ù úê ú úûêë -T0 (z) 1 úû Filters to Z-transform 1 H 0 (z1/2 )X(z1/2 ) + H 0 (-z1/2 )X(-z1/2 )) ( 2 1 Y1 (z) = ( H1 (z1/2 )X(z1/2 ) + H1 (-z1/2 )X(-z1/2 )) 2 Y0 (z) = X ' (z) = G0 (z)Y0 (z 2 )+G1 (z)Y1 (z 2 ) Lifting to Filters é 2 Y (z ) 0 ê ê Y1 (z 2 ) ë ù é 2 X (z ) 0 2 ú = H (z )ê ú ê X1 (z 2 ) û ë ù ú ú û Filter Results = Polyphase Lifting H 0 (z) = H 00 (z 2 ) + zH 01 (z 2 ) H1 (z) = H10 (z 2 ) + zH11 (z 2 ) G0 (z) = G00 (z 2 ) + z -1G01 (z 2 ) G1 (z) = G10 (z 2 ) + z -1G11 (z 2 ) Further reading Boundary problems! Vanishing Moments! Wavelet packets Second generation wavelets Multiwavelets Curvelets, ridgelets … Any Questions?