The Wavelet Tutorial Part IV

转自http://users.rowan.edu/~polikar/WAVELETS/WTpart3.html

（先传上来，等自己学习时在写些心得，感觉是一种学习的好方式，希望大家能够帮我指点误区，谢谢啦，网络是一个好工具，希望能够遇到志同道合者，嘿嘿）

The Wavelet Tutorial

Part IV

by

ROBI POLIKAR

MULTIRESOLUTION ANALYSIS: THE DISCRETEWAVELET TRANSFORM

 

Why is the Discrete Wavelet Transform Needed?

Although the discretized continuous wavelet transform enables the computation of the continuous wavelet transform by computers, it is not a true discrete transform. As a matter of fact, the wavelet series is simply a sampled version of the CWT, and the information it provides is highly redundant as far as the reconstruction of the signal is concerned. This redundancy, on the other hand, requires a significant amount of computation time and resources. The discrete wavelet transform (DWT), on the other hand, provides sufficient information both for analysis and synthesis of the original signal, with a significant reduction in the computation time(出于信息的冗余和计算时间的考虑).

The DWT is considerably easier to implement when compared to the CWT. The basic concepts of the DWT will be introduced in this section along with its properties and the algorithms used to compute it. As in the previous chapters,examples are provided to aid in the interpretation of the DWT.

 

THE DISCRETE WAVELET TRANSFORM (DWT)

 

The foundations of the DWT go back to 1976 when Croiser,Esteban, and Galand devised a technique to decompose discrete time signals. Crochiere, Weber, and Flanagandid a similar work on coding of speech signals in the same year. They named their analysis scheme as subband coding.In 1983, Burt defined a technique very similar to subband coding and named it pyramidal coding which is also known as multiresolution analysis. Later in 1989, Vetterli andLe Gall made some improvements to the subband coding scheme, removing the existing redundancy in the pyramidal coding scheme. Subband coding is explained below. A detailed coverage of  the discrete wavelet transform and theory of multiresolution analysis can befound in a number of articles and books that are available on this topic, andit is beyond the scope of this tutorial.

 

The Subband Coding and The Multiresolution Analysis

 

The main idea is the same as it is in the CWT. A time-scale representation of a digital signal is obtained using digital filtering techniques. Recall that the CWT is a correlation between a wavelet at different scales and the signal with the scale (or the frequency) being used as a measure of similarity. The continuous wavelet transform was computed by changing the scale of the analysis window, shifting the window in time, multiplying by the signal, and integrating over all times. In the discrete case, filters of different cutoff frequencies are used to analyze the signal at different scales. The signal is passed through a series of high pass filters to analyze the high frequencies, and it is passed through a series of low pass filters to analyze the low frequencies.

The resolution of the signal, which is a measure of the amount of detailinformation in the signal, is changed by the filtering operations, and thescale is changed by upsampling and downsampling (subsampling)operations. Subsampling a signal corresponds toreducing the sampling rate, or removing some of the samples of the signal. Forexample, subsampling by two refers to dropping everyother sample of the signal. Subsampling by a factornreduces the number of samples in the signaln times.

Upsampling a signal corresponds to increasing the sampling rate of a signalby adding new samples to the signal. For example, upsampling by two refers toadding a new sample, usually a zero or an interpolatedvalue, between every two samples of the signal. Upsampling a signal by a factorof n increases the number of samples in the signal by a factor ofn.

Although it is not the only possible choice, DWT coefficients are usuallysampled from the CWT on a dyadic grid, i.e., s₀ = 2 and t₀ = 1, yielding s=2^jand t =k*2^j, as described inPart 3. Since the signal is a discrete time function, the terms function andsequence will be used interchangeably in the following discussion. Thissequence will be denoted by x[n], where n is an integer.

The procedure starts with passing this signal (sequence) through a half banddigital lowpass filter with impulse response h[n]. Filtering a signal correspondsto the mathematical operation of convolution of the signal with the impulseresponse of the filter. The convolution operation in discrete time is definedas follows:

A half band lowpass filter removes all frequencies that are above half ofthe highest frequency in the signal. For example, if a signal has a maximum of1000 Hz component, then half band lowpass filtering removes all the frequenciesabove 500 Hz.

The unit of frequency is of particular importance at this time. In discretesignals, frequency is expressed in terms of radians. Accordingly, the samplingfrequency of the signal is equal to 2pradians in terms of radial frequency. Therefore, the highest frequencycomponent that exists in a signal will be pradians, if the signal is sampled at Nyquist’s rate(which is twice the maximum frequency that exists in the signal); that is, the Nyquist’s rate corresponds to p rad/s in the discrete frequencydomain. Therefore using Hz is not appropriate for discrete signals. However, Hzis used whenever it is needed to clarify a discussion, since it is very commonto think of frequency in terms of Hz. It should always be remembered that theunit of frequency for discrete time signals is radians.

After passing the signal through a half band lowpass filter, half of thesamples can be eliminated according to the Nyquist’srule, since the signal now has a highest frequency of p/2 radians instead of pradians. Simply discarding every other sample willsubsamplethe signal by two, and the signal will then have half the number of points. Thescale of the signal is now doubled. Note that the lowpass filtering removes thehigh frequency information, but leaves the scale unchanged. Only the subsampling process changes the scale. Resolution, on theother hand, is related to the amount of information in the signal, andtherefore, it is affected by the filtering operations. Half band lowpassfiltering removes half of the frequencies, which can be interpreted as losinghalf of the information. Therefore, the resolution is halved after thefiltering operation. Note, however, the subsamplingoperation after filtering does not affect the resolution, since removing halfof the spectral components from the signal makes half the number of samplesredundant anyway. Half the samples can be discarded without any loss of information.In summary, the lowpass filtering halves the resolution, but leaves the scaleunchanged. The signal is then subsampled by 2 sincehalf of the number of samples are redundant. Thisdoubles the scale.

This procedure can mathematically be expressed as

Having said that, we now look how the DWT is actually computed: The DWTanalyzes the signal at different frequency bands with different resolutions bydecomposing the signal into a coarse approximation and detail information. DWTemploys two sets of functions, called scaling functions and wavelet functions,which are associated with low pass and highpass filters, respectively. Thedecomposition of the signal into different frequency bands is simply obtainedby successive highpass and lowpass filtering of the time domain signal. Theoriginal signal x[n] is first passed through a halfbandhighpass filter g[n] and a lowpass filter h[n]. After the filtering, half ofthe samples can be eliminated according to the Nyquist’srule, since the signal now has a highest frequency of p /2 radians instead of p . The signal can therefore be subsampled by 2, simply by discarding every other sample.This constitutes one level of decomposition and can mathematically be expressedas follows:

where y_high[k]and y_low[k] are the outputs of thehighpass and lowpass filters, respectively, after subsamplingby 2.

This decomposition halves the time resolution since only half the number ofsamples now characterizes the entire signal. However, this operation doublesthe frequency resolution, since the frequency band of the signal now spans onlyhalf the previous frequency band, effectively reducing the uncertainty in thefrequency by half. The above procedure, which is also known as the subband coding, can be repeated for further decomposition.At every level, the filtering and subsampling willresult in half the number of samples (and hence half the time resolution) andhalf the frequency band spanned (and hence double thefrequency resolution). Figure 4.1 illustrates this procedure, where x[n] is theoriginal signal to be decomposed, and h[n] and g[n] are lowpass and highpassfilters, respectively. The bandwidth of the signal at every level is marked onthe figure as "f".

Figure 4.1. The SubbandCoding AlgorithmAs an example, suppose that the original signalx[n] has 512 sample points, spanning a frequency band of zero to p rad/s. At thefirst decomposition level, the signal is passed through the highpass andlowpass filters, followed by subsampling by 2. Theoutput of the highpass filter has 256 points (hence half the time resolution),but it only spans the frequencies p/2to p rad/s(hence double the frequency resolution). These 256 samples constitute the firstlevel of DWT coefficients. The output of the lowpass filter also has 256samples, but it spans the other half of the frequency band, frequencies from 0to p/2 rad/s.This signal is then passed through the same lowpass and highpass filters forfurther decomposition. The output of the second lowpass filter followed by subsampling has 128 samples spanning a frequency band of 0to p/4 rad/s,and the output of the second highpass filter followed by subsamplinghas 128 samples spanning a frequency band of p/4to p/2 rad/s.The second highpass filtered signal constitutes the second level of DWTcoefficients. This signal has half the time resolution, but twice the frequencyresolution of the first level signal. In other words, time resolution hasdecreased by a factor of 4, and frequency resolution has increased by a factorof 4 compared to the original signal. The lowpass filter output is thenfiltered once again for further decomposition. This process continues until twosamples are left. For this specific example there would be 8 levels of decomposition,each having half the number of samples of the previous level. The DWT of theoriginal signal is then obtained by concatenating all coefficients startingfrom the last level of decomposition (remaining two samples, in this case). TheDWT will then have the same number of coefficients as the original signal.

The frequencies that are most prominent in the original signal will appearas high amplitudes in that region of the DWT signal that includes thoseparticular frequencies. The difference of this transform from the Fouriertransform is that the time localization of these frequencies will not be lost.However, the time localization will have a resolution that depends on whichlevel they appear. If the main information of the signal lies in the high frequencies,as happens most often, the time localization of these frequencies will be moreprecise, since they are characterized by more number of samples. If the maininformation lies only at very low frequencies, the time localization will notbe very precise, since few samples are used to express signal at thesefrequencies. This procedure in effect offers a good time resolution at highfrequencies, and good frequency resolution at low frequencies. Most practicalsignals encountered are of this type.

The frequency bands that are not very prominent in the original signal willhave very low amplitudes, and that part of the DWT signal can be discardedwithout any major loss of information, allowing data reduction. Figure 4.2illustrates an example of how DWT signals look like and how data reduction isprovided. Figure 4.2a shows a typical 512-sample signal that is normalized tounit amplitude. The horizontal axis is the number of samples, whereas thevertical axis is the normalized amplitude. Figure 4.2b shows the 8 level DWT ofthe signal in Figure 4.2a. The last 256 samples in this signal correspond tothe highest frequency band in the signal, the previous128 samples correspond to the second highest frequency band and so on. Itshould be noted that only the first 64 samples, which correspond to lowerfrequencies of the analysis, carry relevant information and the rest of thissignal has virtually no information. Therefore, all but the first 64 samplescan be discarded without any loss of information. This is how DWT provides avery effective data reduction scheme.

 

Figure 4.2 Example of a DWT

We will revisit this example, since it provides important insight to how DWTshould be interpreted. Before that, however, we need to conclude ourmathematical analysis of the DWT.

One important property of the discrete wavelet transform is the relationshipbetween the impulse responses of the highpass and lowpass filters. The highpassand lowpass filters are not independent of each other, and they are related by

where g[n] is the highpass, h[n] is the lowpassfilter, and L is the filter length (in number of points). Note that the twofilters are odd index alternated reversed versions of each other. Lowpass tohighpass conversion is provided by the (-1)ⁿterm. Filters satisfying this condition are commonly used in signal processing,and they are known as the Quadrature Mirror Filters(QMF). The two filtering and subsampling operationscan be expressed by

The reconstruction in this case is very easy since halfbandfilters form orthonormal bases. The above procedureis followed in reverse order for the reconstruction. The signals at every levelare upsampled by two, passed through the synthesisfilters g’[n], and h’[n] (highpass and lowpass,respectively), and then added. The interesting point here is that the analysisand synthesis filters are identical to each other, except for a time reversal.Therefore, the reconstruction formula becomes (for each layer)

However, if the filters are not ideal halfband,then perfect reconstruction cannot be achieved. Although it is not possible torealize ideal filters, under certain conditions it is possible to find filtersthat provide perfect reconstruction. The most famous ones are the onesdeveloped by Ingrid Daubechies, and they are known as Daubechies’ wavelets.

Note that due to successive subsampling by 2, thesignal length must be a power of 2, or at least a multiple of power of 2, inorder this scheme to be efficient. The length of the signal determines thenumber of levels that the signal can be decomposed to. For example, if thesignal length is 1024, ten levels of decomposition are possible.

Interpreting the DWT coefficients can sometimes be rather difficult becausethe way DWT coefficients are presented is rather peculiar. To make a real longstory real short, DWT coefficients of each level are concatenated, startingwith the last level. An example is in order to make this concept clear:

Suppose we have a 256-sample long signal sampled at 10 MHZ and we wish toobtain its DWT coefficients. Since the signal is sampled at 10 MHz, the highestfrequency component that exists in the signal is 5 MHz. At the first level, thesignal is passed through the lowpass filter h[n], and the highpass filter g[n],the outputs of which are subsampled by two. Thehighpass filter output is the first level DWT coefficients. There are 128 ofthem, and they represent the signal in the [2.5 5] MHz range. These 128 samplesare the last 128 samples plotted. The lowpass filter output, which also has 128samples, but spanning the frequency band of [0 2.5] MHz, are further decomposedby passing them through the same h[n] and g[n]. The output of the secondhighpass filter is the level 2 DWT coefficients and these 64 samples precedethe 128 level 1 coefficients in the plot. The output of the second lowpassfilter is further decomposed, once again by passing it through the filters h[n]and g[n]. The output of the third highpass filter is the level 3 DWT coefficiets. These 32 samples precede the level 2 DWTcoefficients in the plot.

The procedure continues until only 1 DWT coefficient can be computed atlevel 9. This one coefficient is the first to be plotted in the DWT plot. Thisis followed by 2 level 8 coefficients, 4 level 7 coefficients, 8 level 6coefficients, 16 level 5 coefficients, 32 level 4 coefficients, 64 level 3coefficients, 128 level 2 coefficients and finally 256 level 1 coefficients.Note that less and less number of samples is used at lower frequencies,therefore, the time resolution decreases as frequency decreases, but since thefrequency interval also decreases at low frequencies, the frequency resolutionincreases. Obviously, the first few coefficients would not carry a whole lot ofinformation, simply due to greatly reduced time resolution. To illustrate thisrichly bizarre DWT representation let us take a look at a real world signal.Our original signal is a 256-sample long ultrasonic signal, which was sampledat 25 MHz. This signal was originally generated by using a 2.25 MHz transducer, therefore the main spectral component of thesignal is at 2.25 MHz. The last 128 samples correspond to [6.25 12.5] MHzrange. As seen from the plot, no information is available here,hence these samples can be discarded without any loss of information. Thepreceding 64 samples represent the signal in the [3.12 6.25] MHz range, whichalso does not carry any significant information. The little glitches probablycorrespond to the high frequency noise in the signal. The preceding 32 samplesrepresent the signal in the [1.5 3.1] MHz range. As you can see, the majorityof the signal’s energy is focused in these 32 samples, as we expected to see.The previous 16 samples correspond to [0.75 1.5] MHz and the peaks that areseen at this level probably represent the lower frequency envelope of thesignal. The previous samples probably do not carry any other significantinformation. It is safe to say that we can get by with the 3^rd and 4^thlevel coefficients, that is we can represent this 256sample long signal with 16+32=48 samples, a significant data reduction whichwould make your computer quite happy.

 One area that has benefited the most from this particular property ofthe wavelet transforms is image processing. As you may well know, images,particularly high-resolution images, claim a lot of disk space. As a matter offact, if this tutorial is taking a long time to download, that is mostlybecause of the images. DWT can be used to reduce the image size without losingmuch of the resolution. Here is how:

For a given image, you can compute the DWT of, say each row, and discard allvalues in the DWT that are less then a certain threshold. We then save onlythose DWT coefficients that are above the threshold for each row, and when weneed to reconstruct the original image, we simply pad each row with as manyzeros as the number of discarded coefficients, and use the inverse DWT toreconstruct each row of the original image. We can also analyze the image atdifferent frequency bands, and reconstruct the original image by using only thecoefficients that are of a particular band. I will try to put sample imageshopefully soon, to illustrate this point.

Another issue that is receiving more and more attention is carrying out thedecomposition (subband coding) not only on thelowpass side but on both sides. In other words, zooming intoboth low and high frequency bands of the signal separately. This can bevisualized as having both sides of the tree structure of Figure 4.1. Whatresult is what is known as thewavelet packages. We will not discuss waveletpackages in this here, since it is beyond the scope of this tutorial. Anyonewho is interested in wavelet packages, or moreinformation on DWT can find this information in any of the numerous textsavailable in the market.

And this concludes our mini series of wavelet tutorial. If I could be of anyassistance to anyone struggling to understand the wavelets, I would considerthe time and the effort that went into this tutorial well spent. I would liketo remind that this tutorial is neither a complete nor a through coverage ofthe wavelet transforms. It is merely an overview of the concept of wavelets andit was intended to serve as a first reference for those who find the availabletexts on wavelets rather complicated. There might be many structural and/ortechnical mistakes, and I would appreciate if you could point those out to me.Your feedback is of utmost importance for the success of this tutorial.

Thank you very much for your interest in The Wavelet Tutorial.