# Full text of "Acoustic noise removal by combining wiener and wavelet filtering techniques"

## See other formats

Institutional Archive of the Naval Postgraduate School Calhoun: The NPS Institutional Archive DSpace Repository Theses and Dissertations 1998-06 1. Thesis and Dissertation Collection, all items Acoustic noise removal by combining wiener and wavelet filtering techniques Forney, Fredric D. Monterey, California. Naval Postg http://ndl.handle.net/10945/7845 raduate School Downloaded from NPS Archive: Calhoun atthe | DUDLEY WN] | ciseaRy http://www.nps.edu/library Calhoun is the Naval Postgraduate School's public access digital repository for research materials and institutional publications created by the NPS community. Calhoun is named for Professor of Mathematics Guy K. Calhoun, NPS's first appointed — and published — scholarly author. Dudley Knox Library / Naval Postgraduate School 411 Dyer Road / 1 University Circle Monterey, California USA 93943 oe ate Pe ea ar “aby eh o ig t ; esi oh O se uneteaa sail SE ag ea tlt 1a Gud ta “= eer ee hd By ” CPteLY, ae adm " A as 1d ch te ade sf IAT A eos Prd is es . yh lag ays ni at 6 Ak 1 iat SP ne ls? R oer eee mile tis = arareenas teat, whanes ) cs can wane 5! Res nid teeth sctistiaa mes - 8. 06 pa ata atag | are pei 4 . ral seuss e , athe 7H ae tf tg Bot ty pre v8 ite Aah stein! re) ” 5 2 ous Piotseecaanden heb ‘pate : aves -oten 7 ows of 4s oarft Sinaia ie ois athe 42.8 € 0 6% al oF sigh tidal whe iat ad Fad HF ait a Me By parse See m3 iis ree! . Ae * FORNEY, F. 7 i! ins ae eens i an id ‘a coe ells wae siggy et is roth git dou bats " tet ve Ted tones gestae act has chee mgt J Re gr 0% re waege alec al a@natae aS Heh at Pa hae eee saton m4 ‘a Ten pte ote tp sta. palite wn a {itiige Aevad ue Relor, i Pere goyeuee ae Te Vedas ae He eer a eh ial, eer Beats ut wens Nae ahah d, hat oy ege sey dete tas sy Satan 3 a wae Fab Dec eMal ees See eaia Se ea ee te Tih yt deka oS wha) Gilg2 ete cote Ms ors ae bend 8 a0 wl a ‘ at ae 1 Ya det A Ooty gdeg platen eye tet uns ‘ wo i re tg Nas ran nt gde pita te ytateday! 1pyy ane Bae fa we C4 * > ° it agboeinisginie a Tass bac ches Ch care "3 ais ob LOS hare ad yenee Satiaragniteeey Hoe . eee aes : he athe 3 a a Bane he ange 7g Reet acs ie ry seth ra eg P Se Nekahgahe Pet Tog ahd pa pen et pata obetal sae te ta . v vatehne ats % auledader vase eet a iat ata! Nita dae ae i fabio af ny Sphhhics H as ae + mee ens 4 at 4 + fg lntstated padet® CS re * ata batedateds a acter td ie inte af ae wie tgieat sete» - eva bad eral ef 4 peas 2 igs aa - nary A een i is Be arte + a E , Ce ek Har ah Muted Made Qe biFe B ames 8 et LaSade tencla te otehee 1 So® eee em teh e sFe8s eeatedet it” A a ee MST ser Bie ey Ny Sakae wf Men eh 3a" 1s eae weg srylole'e aSeaedes ony 0 Babe te obst sete ut of Ny tat of Gane te DeMoivie iets vont Od gg pp evade! Seat aoe” watets etattede pee dote i td 9% ee teh G04 wake 4 As wt as sate $s aanieatt.: aie Cae ere Oe ot erate bene® rr ee ie Le "tk rer Ag rs) “pear ~eed Te tetes tea oe yt 8 ee Seueee onade ant a 3 nf . .° = awals & Seted Bie * ‘aoe ares : taSeaeate © oa eae 8 sea haty td A, i pene «= . ns Pie rae ™ hf s entg nN bases = 2, e * hate of Tate ka a sor any fie aa os a ee tergtee a8 gtyaee eee a” te 4 3° ¢ NM Ca Behe 4 il ee ae 2 oe ate . oor. 4 ‘ ig etree re | . pene ‘ . 4h yi ae Per a roe oem 148 rn ,." “uh 2 4A Rar ea ee Gonea gna te ere’ os a se pare ie ee Oye a ge etant” Serene 2 UN aS ane eat tis POON had af bie ry aa poe: Ser anes *, 1%e* - Ay a Ra oats aes sd CaN 1 SP ar ne qa am ue agen ee oF <4 ie ee id a Det paneme cy wa sl Pie fa ry soncaoip ees eta 7105 842g Be + ‘4 ea" “| oe ase PS ath Lokal de tet seers : a : a = zl Pree yt he dh at “86 pachecty aA gre eg Oe. 8 be fe urate ee bad ae a — Re « pe ST) Re eer : . ya ghg ath i voatety ata bk ae t) seryigterei Beer: pee a “7% ee LP etd aren meee Nt ae yay Cae nity Ait ite pakistan fa sass steet om git es at ait incur! ree seit Se fights rae ¥, pepe yong > a tie abe ae cee sto tate wer a hee i Ne PNT Atay an cabetarn secherg tary meet store B08! rs tte at aetiye Eder bho gurees PP ak dee he Sa pa iL Sy Tie Baars Tie it acd bt ae a ihn Le es Pitan ear ahty rate coats fa “| ral mt al urhety tata we Me eap escnepeaetes iw eet ne cl etatensteadant at te? Be st ut a6 taegaree fet ea rages te pick Ter el ey Seat) stab ety" sghyegine hee Nitaeey r Ae seh by ted binest " ak Sieee akong i; .¥ geuaze ba ps4 oe an sone 4 atewutts gents HY anette i ee et ts ; ‘ eh isa cee ie ae es SS eat ananassae on ete wate Rosod ey ace rep Piey ig ern get te) Se Bi be PRE be ah bot ts" A ney preter 74 y tt daptunaute ecerud wens areey” ee % . : fate al Seat i a astitey crea ee a H edtoee" erie Le a et rc te i ory Meee eee a & a) aie ess A che ‘ ahs aa Ae aaa atet Dy ay he aytets "4 eyee em, bade he ord ware fT. ; a ost ‘ ie ake A ‘Sescethe ee tye A amr Pine Sh Reet tt % seat wot erg tar Aa “ene wih a patel a eh4e ie = Hifi en ae eae a es Dita ipls fe seit! ot rT Detetatece tenure ae hearin ie *esetg Re ite = apys Pits PRL Ted ie, pe Ee ty tgeesbrene WZ teteteta wegen Seo ' OL aif ies ’ tater HA whoa ayer pepe Mes me shee ta'n a lela dA i eoraresece cnr tgng-'* ar rereryet SW Serates eT ar, bp Sees tarig's meee ~ es ae etry oe i ries ence po A AR A ae a eR at fahy es eens i be tp siel ge are are efile ett ote? raga lay baal abe tes dite! aa vag ea) Lee tate pee 28 egw TIME STE Se = st wi 48 es ser ppacmaceae jo"s. io” “ tPA é FH < u iat meee . 7 Crit. tid = pray eee cigetley iets et inlet ames ead sites inetses eee whe "ase weeps Tar Ft Py ndagenig tees 1 Ait yiate, eee ster atad ! -y‘e" 4 se "2 eet se Dee att ‘8 Ges oes op ne 2 ed ink ve y re ed be stk Oe wid igi tees De ee Me me Pirie. herds atere Paine > cecets pe ce) baum Spacey f sb aetsee 8 ebatae st itt, ie eet Syagstewtty Aceee” : Uk aa Pity ope pera sere Beare oe i a ao ense NAVAL POSTGRADUATE SCHOOL Monterey, California THESIS ACOUSTIC NOISE REMOVAL BY COMBINING WIENER AND WAVELET FILTERING TECHNIQUES by Fredric D. Forney, Jr. June 1998 | Thesis Advisor: Monique P. Fargues Co-Advisor: Ralph Hippenstiel Approved for public release; distribution is unlimited. : REPORT DOCUMENTATION PAGE | Public reporung burden for this collection of oeeaen is estimated to average | hour per response, mivalreiny: the ume for reviewing instruction, searching existing }} data sources, gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for | Information Operations and eae ae Jefferson Davis ae Suite 1204, Arlington, VA 22202-4302, and to the Office of Management and Budget, AGENCY USE ONLY i blank) REPORT DATE REPORT TYPE AND DATES COVERED June 1998 Master’s Thesis TITLE AND SUBTITLE: ACOUSTIC NOISE REMOVAL BY FUNDING NUMBERS COMBINING WIENER AND WAVELET FILTERING TECHNIQUES 6. AUTHOR(S) Fredric D. Forney, Jr. 7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) Naval Postgraduate School ORGANIZATION Monterey CA 93943-5000 REPORT NUMBER 9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSORING/MONITORING AGENCY REPORT NUMBER 11. SUPPLEMENTARY NOTES The views expressed in this thesis are those of the author and do not reflect the official policy or position of the Department of Defense or the U.S. Government. 12a. DISTRIBUTION/AVAILABILITY STATEMENT 12b. DISTRIBUTION CODE Approved for public release; distribution is unlimited. 13. ABSTRACT (maximum 200 words) This thesis investigates the application of Wiener filtering and wavelet techniques for the removal of noise from underwater acoustic signals. Both FIR and IIR Wiener filters are applied in separate methods which involve the filtering of wavelet coefficients which have been produced through a discrete wavelet decomposition of the acoustic signal. The effectiveness of the noise removal methods is evaluated by applying them to simulated data. The combined Wiener wavelet filtering methods are compared to traditional denoising techniques which include Wiener filtering and wavelet thresholding methods. PERFORMING 14. SUBJECT TERMS Wavelet Analysis, Wiener Filtering, Denoising, Acoustic 15. NUMBER OF Eiris. PAGES 148 16. PRICE CODE 17. SECURITY CLASSIFI- 18. SECURITY CLASSIFI- 19. SECURITY CLASSIFI- 20. LIMITATION OF CATION OF REPORT CATION OF THIS PAGE | CATION OF ABSTRACT | ABSTRACT Unclassified Unclassified Unclassified UL NSN 7540-01-280-5500 Standard Form 298 (Rev. 2-89) i Prescribed by ANSI Std. 239-18 298-102 Approved for public release; distribution is unlimited. ACOUSTIC NOISE REMOVAL BY COMBINING WIENER AND WAVELET FILTERING TECHNIQUES Fredric D. Forney, Jr. Lieutenant Commander, United States N avy B.S., United States Naval Academy, 1983 M.S., The George Washington University, 1991 Submitted in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE IN ELECTRICAL ENGINEERING from the NAVAL POSTGRADUATE SCHOOL June 1998 ABSTRACT This thesis investigates the application of Wiener filtering and wavelet techniques for the removal of noise from underwater acoustic signals. Both FIR and IIR Wiener filters are applied in separate methods which involve the filtering of Wavelet coefficients which have been produced through a discrete wavelet decomposition of the acoustic signal. The effectiveness of the noise removal methods is evaluated by applying them to simulated data. The combined Wiener wavelet filtering methods are compared to more traditional denoising techniques which include Wiener filtering and wavelet thresholding methods. - —— | 2 TABLE OF CONTENTS _ LT INTRODUCTION ... 2.6.6. eee ence eee eee een eens l 1. OCEAN ACOUSTIC NOMS Bo vicar ee nin eect S A. NOISE SOURGES oo: vs. 55 5 cre ere ciel + l. AMDICHEINOISE. .... 2:5. << << 400 cacee eee ee 4 a. Ultra-Low Band (<1 Hz) ....... Sing ep ie PI 4 b. Infrasonic Band (1 to 20 1iz).......... ee ee 4 Cc. Low Sonic Bana(20to Z0UTHz) «ee ss 5 d. High Sonic Band (200 to SU UUIZ) ee ee 5 e. Ultrasonie Band (> 50,000:HzZ)... . 7 eee. . ee 6 Ze SelONoise ec... . Sh es ee eee di B. TRANSMISSION LOSS AND WATER MASS CHARACTERISTICS ee errr re we oh . fd eas Ske 8 c. WOU Sle) (C)} 15] 5 EE sj on coe eee 8 mime ei NEN FUMING gees 8 es ee ee oes no see a eo eee cs en Hi A. MODEIMDES@RIDTIOIN........:.:+-- +... ee. Os es 1] B. FIR WHENER FIDTER DEV ELOPVIE Nig. ee ee ees 2 12 IV. WAVELET ANALYSIS AND MULTIRATE SYSTEMS .................... 19 Vil A. Pelt 1 ecatibes Ee CIN Meena eee OE re ee ee eee ie B. THE CONTINUOUS WAVELET TRANSFORM ...............-. 22 C. THE DISCRETE WAVELET TRANSFORM .................... 4) D. IVETE RATES SwY.S EVES AUN SS ames sas ss ws See 30 Jel DAUBEGHITE SW Ay ECE ro IOS. ss. sds ae S/S) V. DENOISING METHODS USING WAVELET THRESHOLDING ............. 43 A. WAVELET COEFFICIENT THRESHOLDING .................. 43 B. TREES IOI Sees CO aaa ks se es ee 45 l. LOPES EUG CT lore oe ee ee 45 2. SIRE i itesO ld eee ei: ss << s « Se le 46 8. IVD BIG AN TESMNO| dae cepa wi sie... Semone aecy cee ee 46 4. Minimax Primeaple Vntesheld) . . “aiteicc.u as... sss sesso 47 >: Threshold Discussion ... et ee 47 pple ReUNe wey NORIO ER . eee tees ee ee aks eee ay A. M@IPTMIRATE SYSTEMS ANAIEYSIS 2...) . i scugenn 00 ee 60 B. AMPIAS ‘CANCELEA TION 9.) igs ee eee 66 Pe elimliicen te Nite tn WW AVERT TIGDIIR ox 2 eee ss in a 98S ee tie we 2 A. De AIEEE |) 6) 522 ook il 2) See rrr ane Ome 71 B. IIR WIENER FILTER APPLIED TO THE WAVELET DOMAIN _... 73 Vill @F WIENERSHRINK. . peepee 719 VI. COMPARISON OF DENOISING METHODS .......................005. 85 A. WIENER FILTER ........ ... 0.0%... oS tee 86 B WAVELET THRESHOLDIIG & «27ers esc 87 C FIR WIENER WAVELET FILTER .. se ee ere 87 D TTR WIENEReWAVELET FILTER «22a oe ee eee 88 E WIENERSHRINK = oi: esos Se ee. ee 88 iPemveAVELE |) PACKETS: 2... soe Seri ic..... lee ae 101 A WAVELET PACKET THEORY Sue a0 262 een eee 102 B WAVELET PACKETS APPLIED TO TEST SIGNAL ............ 104 l. Short-time Wiener Filter 2... 245... 22 eee. 106 ip Wavelet Ehresholding .........iah.. 2. oe eee ee 106 3. FIR Witeper Wavelct Filter... 0: soe ne 106 4. IR Wiener Wavelet Milter... 0.22... ear see 107 Be Wavelet Packet [hresholdime =.) acer nen es: 107 6. BIR OW ienen Wavelet Packet Rilietc--pa es oe 6 ble a eee 108 7. IIR Wiener Wavelet Packet Filter ..... ics oe 108 8. SUMIRIARY 9)... eens I A Ee i ae 108 PCO CILUSIONS oi cn nt ee eile se ea eee He 1X I. INTRODUCTION Accurate analysis of ocean acoustic signals has proven to be a difficult task due to unwanted noise which is usually present. Significant effort has been directed to the problem of noise removal from underwater acoustic signals. In the 1940s, Norbert Wiener conducted extensive research in the area of linear minimum mean-square error (MMSE) filtering. Provided the spectral characteristics of an additive combination of signal and noise are known, the Wiener filter will operate in a linear manner to yield the best separation of the signal from the noise. Here, “best’”” means minimum mean-square error. The Wiener filter is today a well tested method of noise reduction. Wavelet techniques are, conversely, a more modern approach to noise reduction, although their origin lies in the timeless methods developed by Fourier. Wavelet analysis decomposes the signal into a family of basis functions and provides two significant advantages over the traditional Fourier analysis. First, in wavelet analysis there is a wide variety of basis functions to choose from, and second, a multi-resolution capability is provided in the time-frequency domain which is critical to the identification and elimination of noise in a non-stationary environment. Wavelet techniques have proven to be a viable tool for the denoising of acoustic signals. Recently, a question has been posed concerning the combination of these two invaluable techniques. The approach presented applies a discrete wavelet transform to the noisy signal. Next, rather than thresholding the wavelet coefficients, a Wiener filter is applied to separate the noise from the signal. While this approach may intuitively seem reasonable, there are problems associated with aliasing and perfect reconstruction of the denoised signal which must be considered. Given that the aliasing problem can be dealt with, the results still may not be superior to either Wiener filtering alone or wavelet-based techniques alone but the possibilities are certainly worthy of investigation. In this thesis, an approach is presented which combines the Wiener filter, in both the causal finite impulse response (FIR) and the non-causal infinite impulse response (IR) forms, with wavelet-based techniques. The various noise removal methods are applied to simulated data which offers the advantage of providing ground truth to the analysis. Additionally, the signal-to-noise ratio (SNR) level can be easily modified to evaluate the effectiveness of the denoising algorithms. The combined Wiener wavelet based denoising methods demonstrate promise in the recovery of ocean acoustic signals from the noisy ocean environment in which ships operate. | The problem of denoising underwater acoustic signals is addressed in nine additional chapters summarized as follows. Chapter II presents background information on the characteristics of ocean acoustic noise . In Chapter III the theory and application of Wiener filtering methods is discussed. Chapter IV presents the theoretical development of wavelet analysis and and mutirate systems. Standard wavelet threshold based denoising techniques are presented in Chapter V. FIR Wiener wavelet noise removal methods are developed in Chapter VI followed by IR Wiener wavelet methods in Chapter VI. A comparison of the various denoising methods is presented in Chapter VIII. Wavelet packet denoising techniques are applied to the problem of high frequency signals of short duration in Chapter [X. Finally, conclusions are presented in Chapter X. Il. OCEAN ACOUSTIC NOISE Though significant resources have been devoted to non-acoustic detection methods in the undersea warfare environment, sonar remains the primary and most reliable method of detection. However, as modern technology continues to enable significant reduction in radiated source levels the problem of early detection and classification of underwater acoustic signals gains greater significance. A possible solution lies in the area of improved digital signal processing and filtering methods which would allow detection and classification of signals at lower signal-to-noise ratios. The passive sonar equation states that the source level of the target minus the loss due to propagation through the medium, minus the sum of all interfering noises plus improvement by the spatial processing gain of the receiver, must be greater than the detection threshold for a sonar system to sense the presence of a target with a fifty percent probability of detection [1] SL-TL>NL- DI+DT, (2.1) where: SL_ = Source Level of the target being detected passively TL = Transmission Loss as the signal propagates to the detector NL = Noise Spectrum Level of the ambient noise and self noise in the ocean DI =Directivity of the detecting system DT = Detection threshold for 50% probability of detection and all terms are in dB referenced to 1pPa. A reduction in the noise spectrum level would certainly result in an improved detection probability. The noise spectrum level includes both self noise and ambient noise. aX NOISE SOURCES 1. Ambient Noise Urick identifies ambient noise as that part of the total noise background observed with a nondirectional hydrophone which is not due to that hydrophone and its manner of mounting called “self-noise”, or to some identifiable localized noise source [2]. Ambient noise 1s produced by a variety of sources and may be found in the frequency range from 1 Hz to 100 KHz. The noise frequency range is typically divided into five frequency bands of interest, which are discussed next. a. Ultra-Low Band (<1 Hz) Though little is known about the exact contributors at the lower end of the spectrum, it 1s surmised that these sources include tides and hydrostatic effects of waves, Seismic disturbances, and oceanic turbulence [3]. Tides and waves cause hydrostatic pressure changes resulting in a low frequency, high amplitude contribution to the ambient noise spectrum. The constant seismic activity measured on land extends into the ocean environment causing low frequency, high amplitude contributions which add to those produced by tides and waves. Oceanic turbulence gives rise to varying dynamic pressures which are detected by pressure sensitive hydrophones. b. Infrasonic Band (1 to 20 Hz) This band has gained importance with the emergence of improved low frequency narrowband passive sonar systems. It contains the strong blade-rate fundamental frequency of propeller-driven vessels and its accompanying harmonics. A steep negative spectral slope of 10 dB per octave is common in the region from | to5 Hz. This slope goes positive from 5 to 20 Hz as shipping noise begins to become a more significant factor. In the absence of ship traffic this region continues to fall off and is more affected by wind speed. C. Low Sonic Band (20 to 200 Hz) Studies have shown that distant ship traffic is the dominant source of noise at 100 Hz and has a signineent effect in the low thc band [3]. In areas of low shipping intensity, wind speed continues to be the major factor just as it is in the infrasonic and high sonic bands. Thus, an area of heavy shipping such as the North Atlantic, where on average 1,100 ships are underway, will see a much greater effect than less traveled areas such as the South Pacific. d. High Sonic Band (200 to 50,000 Hz) The well known acoustician V.O. Knudsen conducted extensive studies in this band during World War II. In these studies, he was able to correlate noise with wind speed in the frequency range 500 Hz to 25 KHz. His results, published in 1948, are best summarized by the curves shown in Figure 2.1. These curves show a clear relationship between wind speed (or sea state which isn’t measured as accurately) and spectrum levels. Subsequent studies have shown the spectrum to be flatter in the range 200 to 800 Hz but have generally confirmed Knudsen’s results [4]. Crouch and Burt [5] have developed an expression to model the noise spectrum level in dB which is given as NL(f) = Bf) + 20 n logy V, (2.2) where NL is the noise spectrum level in dB referenced to 1 Pa at frequency f, B(f) is the noise level at a wind speed of 1 knot at a particular frequency, n is an empirical coefficient, 100 Moderate: Light Spectrum Level (dB re 1 micro-Pa} 10" 10 Frequency (kHz) Figure 2.1: Average Ambient Noise Spectra. From Ref. [6]. and V is the wind speed in knots. For n = 1 the noise level increases as 20 log,, V, and the noise intensity will increase as the square of the wind speed. é. Ultrasonic Band (> 50,000 Hz) At frequencies above 50,000 Hz thermal noise becomes the predominant contributor to the noise background. In 1952 Mellen [7] showed theoretically that the thermal noise of the molecules of the sea affects hydrophones and places a limit on their sensitivity at high frequencies. Based on principles of classical statistical mechanics he formulated the following expression for the spectrum level NL in dB referenced to | Pa of the thermal noise at frequency f in kHz given by NL(f) = -15 + 20 log, f- (2.3) Though measurements have been recorded in this band they have not conclusively substantiated the above expression due to excessive shipping noise in nearby ports. No measurements in this band in deep, quiet open ocean water appear to have been made until now. Zs Self Noise Self noise includes all noise created by the receiving platform and usually falls into one of two categories which include equipment noise and platform noise. Equipment noise includes electronic or thermal noise produced within the sonar electronic system. Platform noise is produced from the same sources as radiated noise except that the source of platform noise is the receiving platform vice the source or target platform. These sources include machinery noise, hydrodynamic noise, propeller noise and transients. Platform noise reaches the receiving transducer by a variety of methods including vibration via an all-hull path, all- water direct path, all-water back scattered path from volume scatterers, and all-water bottom reflected path [2]. Machinery noise occurs principally as low frequency tonals which are relatively independent of speed. Hydrodynamic noise which includes all sources of noise resulting from the flow of water past the hydrophone and its support and the outer hull structure of the vessel, becomes more significant as speed increases. At high speeds propeller noise becomes the dominant contributor to self noise. B. TRANSMISSION LOSS AND WATER MASS CHARACTERISTICS Properties of the water mass such as temperature, salinity, and density affect the transmission path of sound in water and therefore alter the signal received at the hydrophone. Additionally, the depth of water and bottom structure influence the path traveled and could result in multiple transmission paths between the source and receiver. The affects of absorption and attenuation cause the ocean to filter out the high frequency spectrum while passing the low frequency spectrum. Thus, ocean acoustic noise appears to occur more predominantly in the lower frequency region. CG NOISE MODEL Modeling ocean acoustic noise as a stationary white Gaussian random process substantially simplifies the analysis and study of denoising. However, the variety of sources which contribute to self and ambient noise and the acoustic transmission characteristics of the ocean result in a noise contribution which is colored vice white in nature. A common method around this obstacle is to pre-whiten the acoustic signal. Stationarity is another assumption which must be applied carefully. Certainly the Ocean environment is one in which sources of acoustic signals change regularly. Moreover, the water mass properties are in a constant state of flux. This hardly meets the criterion for Stationarity at first glance. However, on a time scale of several seconds, the ocean environment can indeed be considered stationary. Finally, evaluation of ocean acoustic noise reveals that the assumption of a Gaussian random process appears to hold true in most cases. An artificial computer generated white Gaussian noise sample and its power spectral density are pictured in Figure 2.2. This assumption can be proved by comparing the histogram of a noise sample to the Gaussian probability density function with appropriate sample mean and variance or by checking a normal probability distribution plot of the noise sample as seen in Figure 2.3. Barsanti’s results have verified this assumption with many actual ocean acoustic signals [8]. Frack discusses several more involved methods of verifying this critical modeling assumption [9]. All noisy signals in this thesis are generated by adding white Gaussian noise with zero mean (as shown in Figure. 2.2) to the noise-free signal. When analyzing actual ocean acoustic signals a noise sample must be available to determine the statistical properties necessary to produce optimal filtering. A noise sample with statistical properties which accurately reflect the noise embedded in the noisy signal produces a more effective filter and amore accurate denoised signal. Frack discusses methods which can be used to model noise for filtering purposes when adequate noise samples are not available [9]. Noise Amplitude Oo 100 200 300 400 500 600 700 800 900 1000 Sample PSD of Noise Amplitude (dB) O O.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Normalized Freq Figure 2.2: Artificial White Gaussian Noise Sample. Histogram of Noise and Superimposed Normal PDF Normal Probability Plot of Noise 8:88 = 885 B 0.50 = es Bee | = —-2 —1 Oo 1 2 3 Data Figure 2.3: Histogram and Normal Probability Plot of Artificial White Gaussian Noise Sample. The superimposed normal PDF has been scaled to match the histogram for 1024 samples. 10 III. WIENER FILTERING Wiener filtering is one of the earliest methods used to separate noise from a desired signal. Developed by Norbert Wiener in the 1940's, this form of optimal filtering is used to denoise ocean acoustic signals with the objective of improving signal detection and classification. The FIR Wiener filter and its application to denoising is developed in this chapter. A. MODEL DESCRIPTION The optimum linear discrete-time estimation model is illustrated below in Figure 3.1. The filter input signal x(n), consists of a signal s(m) and an additive noise w(n). The estimator is constrained to be a linear filter with impulse response h(n), which is designed to remove the noise while preserving the characteristics of the signal. The filter impulse response is obtained by minimizing the estimation error e(n), defined as the difference between the filter output and the desired signal d(n), which is taken as the original signal s(n), for filtering applications. , Desired Signal lige mesposié Discrete-Time ; 4 Filter E 7 a) + 4 7 *@)=s@)+wa) 4) y(n) = (2) + d(n) = s(n) ; Estimation Noise Error w(n) e(n) = s(n) - s(n) Figure 3.1: Linear Discrete-Time Estimation Model. 1] A time-varying filter can be derived for nonstationary signals but the complications are significant [10]. A more tractable approach involves segmenting the input signal into blocks which can be considered stationary for the short duration of a few milliseconds to a few seconds. In developing the filter for this model, the input signal and noise spectral characteristics or equivalently, auto- and cross-correlation functions must be known. Although typically not known a priori, this information can be estimated from the data assuming that the noise can be isolated sometime during the observation interval. B. FIR WIENER FILTER DEVELOPMENT Given the spectral characteristics of an additive combination of signal and noise, Wiener proposed a scheme which best separates signal from the noise [11]. His procedure involves solving for the “best” filter coefficients by minimizing the difference between the desired ideal noise-free signal and the filter output in a mean square error sense. The mean square error is most often used due to its mathematical simplicity, as it leads to a second order cost function with a unique minimum obtained by using optimal filter coefficients [10]. The FIR Wiener filter is constrained to length P with coefficients h, (0 < k < P-1). Its output depends on the finite input data record x(n) = [x(n), x(n-1), ......x(n-P+1)]’ and may be expressed as the convolutional sum: P-1 §(n) = y(n) = s h *(k)x(n-k). (3.1) The mean square value of the error, e(n), between the desired output s(m), and the filter output S(7), can then be expressed as: P-1 o-=E{le(n)?} =E{Is(n)- > h*(k)x(n-k)P}. (3.2) k=0 The principle of orthogonality states that the optimal filter coefficients h(n), for n = 0,]1...., P-1, minimize the mean square error if chosen such that E{x(n-i) e'(n)} =0, fori = 0,1...., P-1, that is if the error 1s orthogonal to the observations [12]. The minimum mean square error can then be given by o*, = E{s(n) e'(n)}. Now applying the orthogonality principle results in the following set of equations: E{x(n-ie *(n)} | xn-i| 5 “(n) y+ h(k)x oh) | ="05" “forr=Onterre (B38) k=0 Equation (3.3) can also be expressed in terms of the auto-correlation function of the observations R,.(k) and the cross-correlation function between the signal and the observation sequence R,,(i) which leads to: x ieee?) (eeu), for i= 0,],,...,P-1. (3.4) k=0 Assuming the signal and noise are independent and have a zero mean, it can be shown that: R (A=R (A) ES) and R(A)=R (A) +R (4), where R,(k) and R,,(k) are the signal and noise correlation functions respectively. The set of discrete Wiener-Hopf equations for the causal stationary filter can be expressed in matrix form as Kee © _., (3.6) where R, = E{x(n)x "T(n)} (3.7) and r.. = E{s(n)x(n)} (3.8) 66 699 ow an represents the reversal operator applied to a vector [12]. Equation (3.6) may be 13 solved for the filter coefficients using matrix methods. The minimum mean square error is found using the second part of the orthogonality theorem: o. = E{s(n)e*(n)} = | sn| s w-¥ h encr-b (3.9) or again in terms of the correlation function: o. = R{(0) >? h*(k)R_(k). (3.10) The Wiener filter is now applied to a synthetically generated noisy sinusoidal signal of varying frequency with a SNR of 0 dB. The original signal and noisy signal, shown in Figure 3.2(a), are compared to the original signal and denoised signal in Figure 3.2(b). The standard measure used to compare denoised signals to the noise free original signal is a modified Mean Squared Error (MSE), defined as: N MSE = > a el (3.11) 1 norm(s) norm(y) where s(7) is the noise free original signal, y(n) is the denoised output, and N is the length of the signal. The signals are energy normalized by dividing by their norms. This represents a denoised signal amplitude gain which is applied to account for losses incurred during the filtering process. This normalized MSE will be used throughout the thesis as a measure of performance. It is not normalized by the signal length. However, a standard signal length of 1024 or 16384 points is used exclusively for all denoising trials. The Wiener filtering results for various filter orders is depicted in Figures 3.3 and 3.4. A single-frequency noisy sinusoidal signal (0 dB SNR) is denoised through ten trials and an average normalized MSE 14 is computed. This is repeated for test signals with normalized frequencies which range from 0.01 to 0.49 in steps of 0.01 and the results are displayed. The noise sample supplied to the Wiener filter is an independent noise sample in Figure 3.3 and the original noise sample used to create the noisy sinusoidal signal in Figure 3.4. Although the original noise sample is not practically available, it is shown for comparison purposes. A table comparing the MSE values averaged across the spectrum is shown following the Figures. The MSE performance improves for higher filter orders. Additionally, the MSE value is rather flat across the spectrum for a filter orders higher than four. These results are as expected for the Wiener filter and provide a benchmark standard of performance against which other denoising methods can be compared. 15 Original Signal and Noisy Signal: SNR = OdB, Normalized Freq = O.1 Amplitude —3 Ss 10 18 20 25 30 35 40 45 50 Time (sample number) Original Signal and Wiener Filtered Signal: Filter Order = 12 3 2 a | <b> = = e ° << ~—1 —2 =—3 oO Ss 10 18 20 25 30 35 40 45 50 Time (sample number) Figure 3.2: Wiener Filter (a) Original (dashed) and Noisy Sinusoidal Signal at normalized frequency=0.1, sampling frequency=1. (b) Original (dashed) and filtered signal. 16 Wienr Filter Results: SNR = 0dB 0.45 Filter Order 0.3 S us © 3S 0.25 Cop) = cS i: 8) > 0.15 0.05 Oo 0.05 O.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 Normalized Freq Figure 3.3: Mean-square error comparison for Wiener filters of various order. Independent noise sample. Ten trial average. ly VViener Filter Results: SNR =» Oodb Flliter Order Mean Square Erros 3 0.35 o.4 0.45 o.5 Oo o0.0S o.1 o.1S eo.2 0.258 oO. Normalized Freq Figure 3.4: Mean-square error comparison for Wiener filters of various order. Actual noise sample. Ten trial average. Independent Noise Average MSE 18 IV. WAVELET ANALYSIS AND MULTIRATE SYSTEMS Wavelet analysis has found a myriad of applications in areas ranging from financial research to engineering analysis. It is particularly well suited for signal processing applications such as image compression, sound enhancement, and statistical analysis [13]. This chapter presents the underlying theory of wavelet analysis including both the mathematical basis, which stems from Fourier analysis, and the multirate implementation in filter banks. A. FOURIER ANALYSIS Fourier and wavelet analysis are similar in that the goal is to express an observed real time-series signal, x(t), as a linear decomposition of the form x(t) = >> a, €,(0), (4.1) k where k is an integer index for the finite or infinite sum, a, are real-valued expansion coefficients, and &,(t) are a set of functions called the expansion set [14]. If an appropriate expansion set is chosen, the resulting expansion coefficients will provide useful information about the signal allowing better analysis and processing. The expansion set is called a basis for the signal space if every signal in that space can be expressed as a unique set of linearly independent functions. Further, if the basis functions are orthogonal so that their inner product is zero: GLOLEIG = i E(te,(t)dt = 0, j+k (4.2) then the coefficients can be calculated by the inner product 1) a, = (x(t),€,(f) = [ x(t) €;(t) dt. (4.3) In addition to simplifying the computation of the expansion coefficients, an expansion utilizing an orthogonal basis function results in a sparse representation where most of the coefficients are zero or near zero. This has great utility in the nonlinear noise reduction methods developed in later chapters. A well known transform which is used to extract spectral information from time- series signals is the Fourier transform, NG) [xe a dg (4.4) which uses orthogonal sinusoidal ee functions contained in the complex exponential function. Since the Fourier transform is integrated over all time, it is only suitable for stationary signals which have no time varying frequencies. However, most real-world signals are non-stationary, and therefore, a time-varying spectral representation is essential in conceptualizing the temporal relationship between the spectral components. The short-time Fourier transform (STFT) is the key to adapting the Fourier transform to time-frequency analysis. Let x(#) be the nonstationary signal of interest. A temporal window, w(t-T), is applied to the signal. The window, which is centered at t, segments the signal into short intervals which can be considered stationary. The traditional Fourier transform is then applied to the windowed signal. Thus, the STFT may be expressed as follows: peer (c.f ) -f x(t) w(t-t) exp(-j2mft)dt. (4.5) As a result, the STFT maps a one-dimensional time series into a two-dimensional 20 time-frequency function STFT(t,f). Equation (4.5) can be considered the inner (scalar) product of the signal x(t) with a two parameter family of basis functions w(t-t) exp(-j27/t) for varying t and f [15]. The squared modulus of the short-term Fourier transform of a signal is known as a spectrogram and is expressed as SPEC(t,f ) = | STFT (tf) 1*. The spectrogram is a measure of the signal energy in the time-frequency plane and provides powerful insight when analyzing nonstationary signals. A common band-pass filter implementation of the STFT is shown in Figure 4.1. In this model, the frequency f might be considered the midband frequency of the bandpass filter. This model results in a filter bank setup with a set of analysis filters that perform a STFT for each specific value of frequency f [15]. Band-pass Filter x(t) w(-t)exp(j2xft) STFT( tf) f=nf, n=(1,2,..) exp( -j2xft) Figure 4.1: Band-pass Filter Implementation of STFT. In terms of time and frequency resolution the STFT must satisfy the time-bandwidth limitations of the uncertainty principle, namely: AtAf > — (4.6) 21 where At is the time resolution and Af is the frequency resolution or bandwidth. This equality can be met by using a Gaussian window [16]. However, once the window is specified a constant bandwidth filter bank results with fixed time-frequency resolution over the entire time-frequency plane. The time-frequency plane has a uniform tiling as shown in Figure 4.2. Although the Gaussian window minimizes the time-bandwidth product, the STFT still doesn’t provide sufficient time resolution without a significant degradation in frequency resolution. A narrow or compactly supported window function provides better time resolution at the expense of poorer frequency resolution. Similarly, a wider window function results in better frequency resolution but poorer time resolution. An alternative approach to the STFT is required, one which analyzes the signal at different frequencies with different resolutions. B. THE CONTINUOUS WAVELET TRANSFORM The wavelet transform is a multi-resolution analysis (MRA) method designed to STFT Time-Frequency Display Frequency (A f) A) > N 0 1 2 3 4 5 Time (A t) Figure 4.2: STFT Time-Frequency Display. 22 provide good time resolution and poor frequency resolution at high frequencies and good frequency resolution and poor time resolution at low frequencies. Fortunately, many practical ocean acoustic signals are composed of high frequency components of short duration plus low frequency components of long duration, and thus are well suited to this form of analysis [17]. In terms of the filter bank model described earlier, the time resolution must increase as the central frequency of the analysis filters increases. The filterbank is then composed of band-pass filters with constant relative bandwidth otherwise known as constant- Q filters where the constant relative bandwidth is defined as: Sf = constant = Q. (4.7) f A clearer understanding may be obtained by comparing the division of the frequency domain, as shown in Figure 4.3 below. The STFT analysis filters in Figure 4.3(a) are regularly spaced over the frequency axis while WT multi-resolution analysis filters in Figure STFT with Constant Bandwidth Filters Amplitude 0 {. Zt, 3f, 4f, 5f, 6f, ce Sf, Frequenc Figure 4.3(a): STFT decomposition with Constant Bandwidth Filters. 2 WTF with Constant-Q Factor Filters Amplitude O f, 2f, 4f, 8f, Frequency Figure 4.3(b): WT decomposition with Constant-Q Factor Filters. 4.3(b) are regularly spread in a logarithmic scale. The improved frequency resolution at low frequencies and better time resolution at higher frequencies 1s clearly illustrated for the WT case. The continuous wavelet transform (CWT) is a powerful multi-resolution analysis method which has become a widely used alternative to the STFT. It may be viewed as a signal decomposition onto a set of basis functions where the basis functions are called wavelets [15]. The CWT has some similarities to the STFT in that a segmented time-domain signal is transformed by multiplying it by a wavelet function which is similar to a windowed function in the STFT. However, the basis functions for the two transforms are quite | UN) ee Figure 4.4: Basis Functions - Sine Wave and Daubechies-8 Wavelet. 24 different. Whereas, the Fourier transform breaks the signal into a series of sine waves of different frequencies, the wavelet transform breaks the signal into its “wavelets”, scaled and shifted versions of the “mother wavelet”. The two basis functions are illustrated in Figure 4.4, where a sine wave and a Daubechies-8 wavelet are compared. The smooth sine wave of finite length is a direct contrast to the wavelet which is irregular in shape. Its irregular shape improves the analysis of signals with discontinuities or sharp changes, while its compact support allows temporal localization of a signal’s features [18]. The two-parameter family of basis functions known as wavelets may be expressed I faut as: (= ( — | (4.8) ? Ja a where a is a Scale factor or dilation parameter and t is a time delay. This leads to the definition of the CWT: wy { t-t CWT (t,a) =— [ x(ONw [2 dt. (4.9) a aa a In contrast to the Short Time Fourier transform, the CWT offers improved frequency resolution at low frequency while the time resolution is improved at higher frequency. _ Figure 4.5 illustrates the time-frequency plane for the CWT. The wavelet coefficients, determined by the translation and dilation operations shown in Equation (4.8), represent the 25 Wavelet Time-Frequency Display Frequency (A f) Time (A t) Figure 4.5: CWT Time-Frequency Display. correlation between the wavelet and a localized section of the signal. The process of translation or shifting and dilation or scaling is illustrated in Figure 4.6. AN fi Nin VAY AA a) ia WA Wy Wry ¥ Wy "V a WC!) WE) Ww (¢-3k) Shi fling —_—_» Shifting Shi fringe : -__oo Scaling | Shifting Repeat Shifting Operation NIN WAY IN rfAAANIA AANA NYS VY NYVVY VV KV VWVVY ‘ail WU(E/ 2) Figure 4.6: Scaling and Shifting Process of the Wavelet Transform (Used with permission from Joshua Altmann, [18]). 26 G THE DISCRETE WAVELET TRANSFORM The advantages of multi-resolution analysis and the CWT have made the wavelet transform an invaluable tool in signal processing. However, real-time signal processing and denoising operations using the CWT would require a substantial expenditure in computational time and resources [19]. A much simpler implementation of the wavelet transform is available which reduces the redundant information experienced in CWT analysis. Based on work by Croisier, Esteban, and Galand [20], the discrete wavelet transform (DWT) is an efficient discrete signal processing technique which lends itself to digital computational methods. The DWT efficiently performs signal decomposition and reconstruction. A complete mathematical interpretation of wavelets is based on functional analysis, which is well defined in reference [21]. Basic concepts from both functional and multi- resolution analysis, as discussed in reference [14], will be used to define the scaling function and the wavelet function and describe their significance in the DWT of a signal. The objective, expressed by (4.1), 1s to represent a signal in a given vector space as a linear combination of basis functions in a transformed vector space. These basis functions include dilated and shifted versions of the orthonormal scaling and wavelet functions. A set of scaling functions which span an arbitrary vector subspace V, may be defined as 62) = Ot-k, keEZ, be L’, (4.10) where Z denotes the set of all integers and L’ is the space of square integrable functions. Equation (4.1) may be rewritten as, ag x(t) = X a,0,(t) for any 5 (an fe (4.11) The subspace spanned by the scaling function may be expanded to V; by forming a set of scaling functions which are both scaled and translated as follows 2) = 2? b(2/t-k). (4.12) Consequently, (4.11) becomes x(t) = yy a,,,(t) for any Gls VA (4.13) For j > 0, @,, is shorter and translates in smaller steps allowing the representation of signals with finer details. A nested set of spanned spaces 1s defined as * CV_CV_ CV CV CVi6 oc ibe (4.14) or Vic Vi forall j € Z. (4.15) The recursive relationship ot) = Yo h(n) ¥2b(2t-n), ne Z (4.16) allows the lower scale (2) to be expressed in terms of a weighted sum of the shifted higher scale (21). The coefficients h(n) are the scaling function coefficients which form a filter to be used in the DWT. Considering the scaling function as a coarse approximation of a signal, the wavelet function, which spans the differences between the spaces spanned by the various scales of the scaling function, provides the finer details. The two functions are orthogonal at a given scale and the wavelet function can be expressed in terms of the weighted sum of shifted scaling function (2f) defined in (4.16) by 28 WO) = DY) a(n) ¥2H(2r-n), 1 € Z, (4.17) where g(n) are the wavelet function coefficients. The scaling function coefficients and wavelet function coefficients are required by orthogonality to be related by g(n) = (-1)" h(N-1-n), (4.18) where JN is the finite even length of the signal. The relationship Wi) = 2!" W(2/t-k), (4.19) which is similar to (4.12), determines the two-dimensional family of wavelet functions by scaling and translating the mother wavelet given in (4.17). Figure 4.7 below depicts the relationship between the scaling function and wavelet function spaces. With the scaling and wavelet functions defined, any signal x(#) € L? (R) can be written W,1W,1W,1V, V2 ),2).2. LQ Figure 4.7: Scaling function and wavelet function space orthogonal relationship [14]. as the discrete series xp =) c, (2 *6(2"t-) + > d (kK) y(2t-b). (4.20) 7 é kK J=Jo The choice of jy sets the coarsest scale spanned by the scaling function @,,(¢) providing a Des coarse approximation of the signal. The rest of the space is spanned by the wavelet function which provides the high resolution details of the signal. The set of coefficients from (4.20) is called the discrete wavelet transform of the signal. D. MULTIRATE SYSTEMS ANALYSIS Multirate systems and filter banks are the workhorse of the DWT. A deeper understanding of their role in signal analysis and reconstruction is essential in exploring the effects of —_—e the Wiener filter with wavelet analysis. Multirate system components used to implement the DWT include decimators, interpolators, analysis filters, and synthesis filters. These operations will be discussed in the context of the DWT and then a simple multirate system will be analyzed. Decimation involves the subsampling or downsampling of a discrete sequence by a factor of two for the DWT. This is equivalent to discarding every other sample and results in reducing the sampling rate by a factor of two. Mathematically this process is represented by BA) 24) (4.21) Taking the Fourier transform of (4.21) yields the frequency domain expression for decimation X,(w) = -[X(S) + X(S+n)). (4.22) The Z-transform of (4.21) provides the Z-domain expression X12) = s1X@")+X(-2'")]. (4.23) The frequency and Z-domain expressions show that in addition to reproducing the input at 30 half the sampling rate, X(>), decimation produces a second term shifted by 7 radians, X(>+T). This second term is responsible for aliasing which causes a distorted frequency response unless the original signal is properly bandlimited. If X(w) is bandlimited to|w| < = , then there is no overlap between the two terms; and the aliasing term (i.e. shifted term) can be removed by filtering. Interpolation, also known as expansion, upsamples a discrete sequence, again by a factor of two for the DWT. The insertion of a zero between each sample in the sequence leads to a doubling of the sampling rate. Mathematically this is expressed as x(n/2) for even n 2) = 0 for odd n eee The Fourier and Z-transforms of (4.24) produce the following results X,,(W) = X(20) os (4.25) X,(z) = X(z°). Doubling of the sampling rate causes the original frequency response to be compressed by a factor of two. Additionally, interpolation causes an effect known as imaging where one input frequency causes two output frequencies, one at = radians and a second image at =n radians. This is different from aliasing where two input frequencies, @ and w + 7 result in the same output [22]. Block diagrams of the decimation and interpolation operators are Decimator Interpolator Figure 4.8: Decimation and Interpolation Operators. shown in Figure 4.8. The analysis/synthesis filter bank employs the decimation and interpolation operators in sequential order in order to recover the original sequence, however, the effects of aliasing and imaging must be considered. Decimation followed by interpolation produce the following frequency and Z-domain results Xiryin(@) = [X(@)+X(o+7)] i (4.26) X12y12)(2) = 5lA@) +X(-z)]; which represent a scaled reproduction of the original sequence plus an additional term caused by aliasing and imaging 122), A lowpass filter H,(z) and highpass filter H,(z), known as analysis filters, are now introduced prior to decimation in order to bandlimit the input signal. These maximally decimated FIR filters form a Quadrature Mirror Filter (QMF) bank in which the filters exhibit frequency responses which are a mirror image of each other as shown in Figure 4.9. Since the analysis filters have a nonzero transition bandwidth and stop-band gain, some Quadrature Mirror Fitters: Daubechies~-20 Wavelet Nomalzed Gain oO 0.0$ 0.1 0.18 0.2 0.2 0.35 0.4 0.45 0.5 5 0.3 Normalizad Frequancy Figure 4.9: QMF bank: Analysis Filters H) and A. a7 aliasing does occur due to the decimation and must be eliminated by other means. The imaging problem is remedied by a second QMF bank following the interpolation operation. These filters, known as synthesis filters, remove the images produced by the interpolation operation. However, they suffer from the same non-ideal characteristics as the analysis filters and thus are not capable of completely removing all images. If the analysis and synthesis filters meet specific conditions, they are able to remove all aliasing and imaging and perfectly reconstruct the signal. The analysis and synthesis filters as well as the decimator and interpolator operators are now assembled to form a filter bank as shown below in Figure 4.10. Perfect reconstruction is defined mathematically as A) NOx (Ham) (4.27) A x(n) y,(n) v,(n) u,(n) x(n) Input Analysis Decimators Interpolators Synthesis Output Figure 4.10: Two-channel Quadrature Mirror Filter bank [23]. where c is a non-zero constant and n, a positive integer. Thus x(m) is merely a scaled and delayed version of x(n) [23]. A Z-transform analysis of the filter bank in Figure 4.10 provides the conditions necessary for perfect reconstruction. The output of the analysis filter 38 for channel k (k = 0,1) 1s expressed as amit a) ix(Z) (4.28) According to (4.23) decimation produces the result AOS AGEO) (4.29) = =A, (2 I NGG ne) +H ( ~7 x ~7 at uh The interpolator output, determined by applying (4.25) and (4.26), becomes Uz) = V,{z*) (4.30) = ~[H(2)X@) +H (-z)X(-2)]. Finally, the reconstructed signal, which is the sum of the synthesis filter bank channels, is X(2) = ~(Fo()Ho(2)+F | @H,@)X@)+{Fo(2)Ho(-z)+F @)A(-2))X(-z). (4.31) In matrix form, (4.31) becomes (4.32) Kee) X(- = AO ACM ay) H(-2) [Fi The second component in (4.31) is due to aliasing and imaging and produces replica of the Ho) H,@) ts first term shifted 7t radians to the right on the unit circle. It is commonly referred to as the alias term [23]. The following must be true to eliminate aliasing and imaging La) td (ey elt AZ) id eee) OU. (4.33) Additionally, the remaining term must conform to the following to remove distortion caused by the presence of the analysis and synthesis filters F (2H (2+F (H(z) = 227. (4.34) Having met the alias cancellation and no distortion conditions, perfect reconstruction results in the following, where / is the delay for this filter bank based on the filter length: X(z) = z'X(z) (4.35) The analysis and synthesis filters can now be selected based on (4.33) and (4.34). 34 Assuming the following lowpass analysis filter H(z) = y h(n)z ” , the remaining filters are obtained as follows to produce a two-channel perfect ee filter bank [23] TE) ere EI aca) PiZe 2 ae) (4.36) F(z) = z-“H,(z"'). E. DAUBECHIES WAVELET The Daubechies wavelet, named after Ingrid Daubechies of Bell Laboratories, is a common choice for the analysis and synthesis filters, because it possesses several nice properties. First, the filters produced by this family of wavelets are orthogonal with compact Byaport providing many advantages previously mentioned. Second, the frequency response has maximum flatness atw =O andw=m. For a filter of length N, placing half the filter zeros at z= -1 (or W = 7), results in a maximum flat filter of regularity K = N/2. This second property leads to excellent results when Daubechies filters are used for the DWT decomposition and reconstruction of a large class of signals. The following theorem from Daubechies 1988 paper on orthonormal wavelets [24] summarizes her findings: Theorem 1. A discrete-time Fourier transform of the filter coefficients h(n) having K zeros at w = 7 of the form Hw) = cs * Go) (4.37) satisfies IH(w)? + IA(w+n)P? = 2 (4.38) if and only if L(w) =I{(w)l? can be written as 35 L(w) = P(sin* w/2), (4.39) with P(y) = Py) + y* R(1/2 - y), (4.40) K-1 where Py) = aad eg (4.41) k=0 \ and R is an odd polynomial, chosen such that P(y) 20 for O<y<1. There is no explicit expression for the scaling or wavelet functions. A brief explanation of Daubechies development from references [14] and [22] follows. Equation (4.37) is expressed in terms of the frequency response IH(w)l* as follows: lH(w)*I = 0), (4.42) where L(w) = I{(w)I?. In trigonometric terms, (4.42) may be expressed as polynomials in cos(W@): lH(sin?(w/2))? = Icos*(w/2) P(sin?(w/2)) (4.43) which, after the change of variables of y = sin?(w/2)=1-cos?(w/2), becomes lH(y)? = (1-y)*PO), (4.44) where P(y) is an (N-K) order polynomial. A similar derivation for [H(w+m)!* enables (4.38) to be expressed as: (1-y)KP(y)+y KP(1-y) = 2. (4.45) Daubechies then applied Bezout’s theorem with R=0 in (4.40) resulting in minimum length N for a given regularity K (N=K/2) to determine the explicit solution given by (4.41). P() is the binomial series for (1-y)*, truncated after K terms, with degree K-1 and K coefficients: Pf aa Py) = | +Ky oy 24 RK 2) ye = (1-y)*+0(y *). (4.46) 36 Hence the response in terms of y becomes: K-1 IHF = (1-yyk a yt (4.47) The filter frequency response is now expressed as: (P= (Hse) (esa (4.48) where the following transformations were used: y= and iy = (4.49) Equation(4.48) has K zeros at w=m and K-1 zeros due to the binomial term P (@). Transforming (4.48) to the Z-domain using z + z’ = 2 cos(w) simplifies the calculations. The minimum phase transfer function H(z) and its transpose H(-z) are defined as follows: H(z) = Shen mandi ea)e= y h(n)z”. (4.50) - Equation(4.48) becomes: F(z) = H@)H() = 2) — Z p> ad el 1) (4.51) which has 2K zeros at z = -1, half of which belong to H(z), and 2K-2 zeros due to P(z), half of which are inside the unit circle and belong to H(z). The maximum flat Daubechies filter coefficients are solved for by applying spectral factorization methods to (4.51) to obtain A(z). References [14] and [25] provide software methods for determining the scaling and wavelet functions from which the Daubechies wavelet filter coefficients can be determined. The results using MATLAB [25] are shown for the Daubechies-2 and Daubechies-20 wavelet (filter lengths of four and forty respectively) in Figures 4.11 and 4.12. Notice the improved regularity and smoothness as the filter order increases. The DWT is implemented through a series of half-band highpass and lowpass filters a of varying cutoff frequencies which analyze the signal at different frequency resolution. The filtering operations, which repeatedly half the signal spectral content, determine the frequency resolution. Following this filtering operation, the number of signal samples remains unchanged and consequently half the samples are redundant based on Nyquist’s sampling criterion. Thus, the filtering operation is followed by a downsampling operation to eliminate the unnecessary samples. Downsampling a signal by a factor of two corresponds to reducing the sampling rate by a factor of two which results in halving the time resolution. Thus each successive filtering and downsampling operation improves the frequency resolution by a factor of two while the time resolution suffers by a factor of two. The DWT of signal x[n] is determined as follows. The signal is passed through a half band digital lowpass filter with impulse response h[n] and a half band digital highpass filter with impulse response g[n]. The highpass filter is associated with the wavelet function while the lowpass filter is determined from the scaling function. This is equivalent to convolving the signal with the filter impulse response and is expressed as: Viol] = x[n]*h[n] = Do x[k]A[n-k] = D0 ALA) x[n-K] k k Yrigal] = x[n]*gln] = D/ x[k]gin-k] = )/ glk]xin-d1. Cm k k The filtering process divides the frequency band of y[”] in half forming two signals, y,,,,[7] and y,ie,17] which now have redundant information based on Nyquist’s sampling criterion. These signals are then decimated or subsampled by a factor of two without any loss of information. The combined filtering and subsampling process reduces the time resolution 38 Scaling function phi 0.5 1 1:5 2 25 Analysis low-pass filter fe) 1 2 3 Synthesis low-pass filter Oo 1 2 3 0.2 0.4 0 1 2 3 Wavelet function psi 0 0.5 1 io 2 2.5 Analysis high-pass filter Oo 1 2 3 Synthesis high-pass filter Figure 4.11: Daubechies-2 Scaling and Wavelet Functions [25]. Scaling function phi 0.5 0 10 20 30 Analysis low-pass filter 0.6 0.4 0.2 =0:2 =0.4 0 3 6 9 12 15 18 21 24 27 30 33 36 39 Synthesis low-pass filter 0.6 0.4 0.2 0.2 ~0.4 0 3 6 9 12 15 18 21 24 27 30 33 36 39 Wavelet function psi 0.5 0 10 20 30 Analysis high-pass filter 0.6 0.4 Or2 -—0.2 -0.4 0 3 6 9 12 15 18 21 24 27 30 33 36 39 Synthesis high-pass filter 0.6 0.4 0.2 -0.2 ~0.4 0 3 6 9 12 15 18 21 24 27 30 33 36 39 Figure 4.12: Daubechies-20 Scaling and Wavelet Functions [25]. 59 by a factor of two and doubles the scale which improves the frequency resolution by a factor of two. Mathematically this is expressed as Viowlk] = > x[n] h[2k-n] n Yrignlk] = 0 x[n] g[2k-n]. (4.43) n The combined operation of filtering and subsampling is repeated to decompose the signal to some desired level which is limited by the length of the signal. At every level, filtering and subsampling operations produce data with half the number of samples (reducing the time resolution by a factor of two) and half the frequency spectrum bandwidth (which doubles the frequency resolution). The filterbank representation of the DWT depicted in Figure 4.13 can be efficiently implemented, and with a reasonable amount of computational effort most acoustic signals can be decomposed for denoising purposes. The DWT coefficients resulting from the process described above provide a spectral analysis of the signal with a time resolution that varies depending on the level of decomposition. Large amplitude DWT coefficients indicate significant spectral content and the position of these coefficients within the DWT vector provides time localization. The time resolution is precise at high frequencies and gets worse at each successive level of decomposition (while the frequency resolution improves). 40 x[n] Level 1 DWT g[n] he Coefficients | } Level 2 DWT HPF |Coefficierts Level 3 DWT Coefficients Figure 4.13: Filterbank Representation of the DWT [18]. 41 V. DENOISING METHODS USING WAVELET THRESHOLDING David Donoho and his colleagues at Stanford University have led the way in applying wavelet techniques to the theory of statistics. A particularly relevant outcome of this work is a denoising technique termed “wavelet shrinkage” by Donoho and his co-developer Iain Johnstone [26]. The wavelet transform is applied to a noisy data set producing a set of coefficients which are then “shrunk” towards zero using a Soft or hard statistical thresholding method to select the appropriate coefficients which represent the desired signal. The resulting coefficients are then back transformed to the time domain to produce a denoised signal. A brief and informal explanation of the process follows. A more detailed argument can be found in many of Donoho’s writings [26,27,28]. A. WAVELET COEFFICIENT THRESHOLDING Assuming a suitable wavelet basis function is chosen, the DWT decomposition of a signal will compress the energy of the signal into a small number of large magnitude wavelet coefficients. Gaussian white noise in any one orthogonal basis is transformed by the DWT into white wavelet coefficients of small magnitude. This property of the DWT allows the suppression of noise by applying a threshold which retains wavelet coefficients representing the signal and removes low magnitude coefficients which represent noise. Assume a finite signal s(k) of length N is corrupted by zero mean, additive white Gaussian noise n(k) with standard deviation o, leading to the noisy signal x(k) = s(k) + 6 n(k). (5.1) The objective is to recover the signal s(k) from the noisy observations x(k) using thresholding 43 of the DWT coefficients. The DWT of x(k) is performed using (4.18) and the estimate § is determined by thresholding the individual wavelet coefficients. The three step method for recovery of the desired signal, s(k), from the noisy data, x(k), follows (1) Perform a J-leve] DWT of the data yielding noisy wavelet coefficients w, , , j=j,, oJ , k=0, ... , 2/- 1, where j is the scale or decomposition level, and k is the length of the signal. (2) Apply thresholding in the wavelet domain, using either hard-thresholding of the coefficients defined as: ; Wig lw, | > 6 Wik 7 Tard W429) = 0 Iw < §’ (3.2) or soft-thresholding of the coefficients defined as: > = _ J sign(w,,)(w, 1-6) bw > 6 Mie = Trop Wiw) = 1s Te ISP (3.3) The threshold, 6, is usually determined in one of four ways to be described below. (3) Perform the inverse DWT and obtain the signal estimate S(k). Thus, thresholding cancels the wavelet coefficients which are below a certain threshold value defined in terms of the noise standard deviation o. The artificial signals used in this thesis employ the basic model (5.1) with the standard deviation set at one. The signal is then scaled to provide the desired SNR. In practice the noise standard deviation must be estimated. Standard statistical procedures estimate o as the median absolute deviation of the wavelet coefficients at the finest level J, divided by 0.6745 [8]. Results have shown that the empirical wavelet coefficients at the finest scale are, with few exceptions, essentially pure noise. Although subject to an upward bias due to the presence of some signal at the finest scale, this robust estimate has produced reliable results [26]. A level dependent determination of o is another alternative available if colored noise is suspected. All three options for determining the noise standard deviation are available in the MATLAB Wavelet Toolbox [25]. B. THRESHOLD SELECTION Four commonly used threshold determination methods are discussed next: Universal threshold, Stein’s Unbiased Risk Estimator (SURE), Hybrid threshold, and the Minimax threshold. 1. Universal Threshold The universal threshold, defined as 6,, = of 2logN, (5.4) where o is the standard deviation of the noise and N is the sample length. It uses a single threshold for all detail wavelet coefficients and ensures that asymptotically all detail coefficients are annihilated [29]. The universal threshold is a less conservative method which removes noise effectively but may remove small details of the signal which lie in the noise range [25]. It is easily implemented and is particularly well suited to minimal spectral 45 content signals which have sparse wavelet coefficients [26]. One drawback involves signals of short duration (several hundred samples in length), in which case other methods may provide more accurate results in a MSE sense. 2 SURE Threshold The SURE threshold, 6, , selects a threshold based on Charles Stein’s work on unbiased risk estimators [30]. First implemented by Donoho [27], this method minimizes a risk function to determine an optimum threshold. The unbiased risk estimate used for soft thresholding is defined as follows: RISK(X) = 1 + (X? - 2) I(IXl<8.) + 6; (XI > 6), (5.5) where X ~ N(@,1), J(-) is the indicator function, and 6, is the threshold [29]. Using the model defined by (5.1) the threshold becomes 6. = nin y » RISK “| } (5.6) fe where the wavelet coefficients w,,, divided by the standard deviation have replaced the variable X in the unbiased risk estimate. This technique suffers from the sparse wavelet coefficient problem mentioned for the universal threshold and is often combined with the universal threshold in a hybrid threshold method. 3. Hybrid Threshold The hybrid threshold method is a combination of the first two methods. Developed by Donoho [27], this method uses the SURE threshold unless a low SNR situation with sparse wavelet coefficients develops, in which case, the universal threshold method is utilized. 46 4. Minimax Principle Threshold The inarae principle, often used in statistics to design estimators, uses a fixed threshold, 6,,, selected to produce minimax performance in terms of the mean square error when compared to an ideal procedure [31]. Like the universal threshold, minimax uses a single threshold for all detail wavelet coefficients. The minimax threshold is a function of the sample size, the threshold type (soft or hard), and the oracle type (projection or shrinkage) and can be found in tabulated form in references [28] and [29]. MATLAB Wavelet Toolbox [25] uses the following expression which approximates the minimax threshold for the soft threshold nonlinearity, with comparison to a projection oracle: 6 (2) = 0.3936 + 0.1829 *(log(n)Aog(2)), (5) where n is the sample length. a Threshold Discussion The four thresholding methods are compared in Figure 5.1 where various wavelet thresholding schemes are applied to the same noisy sinusoidal signal used earlier in Chapter Il. The SURE and Hybrid thresholding methods represent the signal as a smoothed sinusoid after a four-level decomposition using the Daubechies-20 wavelet for the DWT. The Universal and Minimax thresholding methods applied to a four-level decomposition do not perform as well as expected, although this signal has minimal spectral content and therefore a Sparse representation in the wavelet domain as seen in Figure 5.2. This is explained by comparing Figures 5.2, 5.3, and 5.4 which show the approximation and detail coefficients before and after universal and Minimax thresholding. The Universal and Minimax methods 47 apply a constant threshold for a given decomposition level. This results in poor denoising performance at any scale other than the coarsest or highest scale otherwise known as the approximation. This is true because the high threshold value results in the signal as well as the noise being removed by the thresholding except for the approximation which is not thresholded in “wave shrink” methods. This effect 1s most noticeable at low SNR values where the signal and noise wavelet coefficients approach the same magnitude. The SURE and Hybrid thresholding methods employ scale dependent thresholds which are more effective in alow SNR environment as demonstrated in Figures 5.5 and 5.6 which show the wavelet coefficients after applying these two methods. The normalized spectral decomposition filters are pictured in Figure 5.7 for the Daubechies-20 wavelet. A MSE comparison of the four thresholding methods across the normalized frequency spectrum from 0 to 0.5 (Sampling Frequency = 1) is shown in Figures 5.8 and 5.9. The poor performance of the Universal and Minimax methods at ” SNR level is clearly evident. Another significant phenomenon noticeable in these figures is the improvement in normalized MSE with each successive scale and decomposition level. This occurs because with each successive scale there is less noise energy as represented by the wavelet coefficients, whereas the signal energy remains the same. The spectral bandwidth is reduced with each level of decomposition eliminating approximately half the remaining white noise which reduces the noise variance and improves the SNR at that particular scale. The energy difference between the signal and noise wavelet coefficients is greater at higher or coarser scales resulting in more effective removal of the noise by thresholding. 48 A closer inspection of the threshold values chosen for a particular thresholding method reveals that the threshold is dependent on the number of decomposition levels even though the signal length of 1024 data point remains the same. This behavior is due to the specific implementation employed in the MATLAB Wavelet Toolbox [25] which doesn’t keep the total number of wavelet coefficients constant. Threshold values for the Universal and Minimax methods are compared below for four levels of decomposition. Level Wavelet Coefficients Threshold Threshold sh 3.7220 22500 3.7607 22596 Since the threshold value for the Universal and Minimax methods is dependent on the number of wavelet coefficients, the threshold varies depending on the decomposition level resulting in independent MSE curves. The SURE and Hybrid methods calculate a different threshold for each scale allowing them to generally outperform the Universal and Minimax methods. This is particularly noticeable at the lower or finest scale where scale dependent SURE and Hybrid thresholding methods generate much better results. 49 Original Signal and Noisy Signal: SNR = 0dB, Signal Freq = .1 5 15 20 25 30 35 Universal Method: MSE = 1.416 fh © th On oO © 10 15 20 25 30 35 SURE Method: MSE = 0.2097 Be © aN On On © 10 ils, 20 25 30 35 40 Hybrid Method: MSE = 0.1939 & oO oO © 10 15 20 25 30 35 40 Minimax Method: MSE = 0.6907 f on Oo © Ampl kK WO 0 5 10 15 20 25 30 35 40 45 50 Time (sample number) Figure 5.1: Comparison of wavelet thresholding schemes; sinusoidal signal at Frequency =0.1, Sampling Frequency F=1, Original signal (dashed); Denoised signal (solid); Five- level DWT decomposition using Daubechies-20 wavelet. Approximation Coeft.: dice SNR=0dB 0 10 20 30 40 50 60 70 80 90 100 Detail Coef.: Level 4 4 2 ; @) ~2 0 10 20 30 40 50 60 70 80 90 100 Detail Coef.: Level 3 10 | apne — -10 0 20 40 60 80 100 120 140 160 180 Detail Coef.: Level 2 0 50 100 150 200 250 300 Detail Coef.: Level 1 0 100 200 300 400 500 600 Figure 5.2: Approximation and detail wavelet coefficients resulting from five-level DWT decomposition of noisy sinusoidal signal, Frequency=0.1 (F=1), SNR=0 dB using Daubechies-20 wavelet. 51 Universal Threshold, Approx. Coef.:Freq=0.1, SNR=0dB 4 2 eee -2 0 10 20 30 40 50 60 70 80 90 100 Detail Coef.: Level 4 thres = 3.761 1 0) —1 0 10 20 30 40 50 60 70 80 90 100 Detail Coef.: Level 3 thres = 3.761 S pee | =, 100 120 80 aa Pay Oe 2 thres = 3.761 1 ae = 100 200 300 Detail Coef.: — 1 thres = 3.761 - 0 100 200 300 400 500 600 Figure 5.3: Approximation and detail wavelet coefficients after applying Universal soft thresholding. Threshold values indicated for each level or scale. Minimax Threshold, Approx. Coef.:Freq=0.1, SNR=0dB 4 2 | ena —2 0 10 20 30 40 50 60 70 80 90 100 Detail Coef.: Level 4 thres = 2.26 0.2 0 0 10 20 30 40 50 60 70 80 90 100 Detail Coef.: Level 3 thres = 2.26 0 20 40 60 80 100 120 140 160 180 Detail Coef.: Level 2 thres = 2.26 0 50 100 150 200 250 300 Detail Coef.: Level 1 thres = 2.26 ; pape = 100 200 300 400 500 600 Figure 5.4: Approximation and detail wavelet coefficients after applying Minimax soft thresholding. Threshold values indicated for each level or scale. © SURE Threshold, Approx. Coef.:Freq=0.1, SNR=0dB 4 2 | AA, -2 0 10 20 30 40 50 60 70 80 90 100 Detail Coef.: Level 4 thres = 1.851 1 0.5 0 0 10 20 30 40 50 60 70 80 90 100. Detail Coef.: Level 3 thres = 0.1886 20 40 60 80 100 120 140 160 180 Detail Coef.: Level 2 thres = 2.26 1 | 1. —1 100 150 200 250 300 Detail Coef.: Level 1 thres = 1.682 8 felt lac rer ~2 100 200 300 400 500 600 Figure 5.5: Approximation and detail wavelet coefficients after applying SURE soft thresholding . Threshold values indicated for each level or scale. =) on © © Hybrid Threshold, Approx. Coef.:Freq=0.1, SNR=0dB 4 2 0 —-2 0 10 20 30 40 50 60 70 80 90 100 Detail Coef.: Level 4 thres = 3.035 1 0 —1 0 10 20 30 40 50 60 70 80 90 100 Detail Coef.: Level 3 thres = 0.1886 10 PS Te 100 120 80 an = — 2 thres = 3.362 1 ee We 100 150 300 Detail Coef.: Level 1 thres = Bec z| 0 100 200 300 400 500 600 Figure 5.6: Approximation and detail wavelet coefficients after applying Hybrid soft thresholding. Threshold values indicated at each level or scale. Four Level Spectral Decomposition Using Daubechies—20 Wavelet Normalized Gain ; \ ae : 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 Normalized Frequency Figure 5.7: Four-level Spectral Decomposition Filters Using Daubechies-20 Wavelet. 56 0.5 Universal SURE owl, 0.9 0.8 Oly a 2 2 9.6 i 2 2 ig) 3) 2.0.5 =. ” ” & 0.4 = ® ® = = 0.3 0.2 0.1 oki 0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 0.3 0.4 OD Normalized Freq Normalized Freq Hybrid Minimax 0.2 0.5 0.18 0.45 0.16 0.4 _ 0.14 =10-35 2 = Ww O.12 W 0.3 = 2 3 i) — 0.1 2.0.25 “) ” @ 0.08 Ss 0.2 ® ® = = 0.06 0.15 0.04 0.1 0.02 0.05 0 ob 0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 0.3 0.4 0.5 Normalized Freq Normalized Freq Figure 5.8: Soft thresholding methods applied to noisy sinusoidal signal: SNR = 10 dB. Four level DWT decomposition using Daubechies-20 wavelet. MSE curves for one through four levels (Lvl) of decomposition or scale pictured. Ten trial average. o/ Universal SURE Mean Square Error VN 78a See Mean Square Error 0 0.1 0.2 OS 0.4 0.5 @) 0.1 0.2 0.3 0.4 0:5 Normalized Freq Normalized Freq Hybrid Minimax Mean Square Error Mean Square Error 0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 0.3 0.4 0.5 Normalized Freq Normalized Freq Figure 5.9: Soft thresholding methods applied to noisy sinusoidal signal: SNR = 0 dB Four level DWT decomposition using Daubechies-20 wavelet. MSE curves for one through four levels (Lvl) of decomposition or scale pictured. Ten trial average. 58 VI. FIR WIENER WAVELET FILTER Wiener filtering and wavelet thresholding are both effective means of denoising acoustic signals. In this chapter the two methods are combined to form an enhanced FIR filter which will be referred to as the Wiener Wavelet filter. This method of denoising utilizes the discrete wavelet transform to represent the signal as a series of wavelet expansion coefficients. Wavelet analysis has several distinct properties which make it an ideal method of signal decomposition, analysis, and reconstruction. First, wavelets are an unconditional basis for a wide variety of signals which means that a signal can be represented by a small number of expansion coefficients, d;, in (4.20), since the magnitudes of these coefficients drop off rapidly for a wide class of signals [32]. Second, wavelet analysis provides a multi- resolution analysis in time and frequency allowing a more accurate description and separation of signal and noise [14]. Third, there are a variety of wavelet functions which makes wavelets adaptable to represent signals of differing characteristics. Finally, wavelet analysis and particularly the DWT is well suited to implementation on a digital computer since it involves only multiplications and additions and not derivatives or integrals. Once the signal has been transformed into the wavelet domain, the wavelet expansion coefficients are filtered using the optimal Wiener filter. In this process the wavelet coefficients are being shrunk towards zero, not by a nonlinear thresholding method but by an optimal linear filtering process. The shrunk coefficients are then inverse transformed to reconstruct the denoised signal. The Wiener filtering process introduces distortion and aliasing which affects the ability to perfectly reconstruct the denoised signal. These effects 3g will be discussed after introducing the necessary multirate system and filter bank theory. A. MULTIRATE SYSTEMS ANALYSIS The application of multirate systems and filter banks in signal analysis and reconstruction was discussed in Chapter IV. The necessary conditions for perfect reconstruction were derived. The optimal Wiener filter will now be applied following decimation and prior to interpolation in each stage of a filter bank. The two-channel filter bank is illustrated below in Figure 6.1. Input Analysis Decimators Wiener Filters Interpolators Synthesis Output Figure 6.1: FIR Wiener Wavelet Filter bank. The Wiener filters which have been inserted in the filter bank above must be provided with the statistical properties of the noise. This is accomplished by isolating a noise-only 60 portion of the noisy signal, x(n). A DWT of the noise, w(7), is performed using the same analysis filter bank structure, as was used for the noisy signal. The Wiener-Hopf equations (3.6) to (3.8) are applied to the noisy signal wavelet coefficients and the noise wavelet coefficients to produce the filter coefficients for the Wiener filters W,) and W,. This two- channel filter bank is repeated to produce a multi-level FIR Wiener wavelet filter bank. A four-level decomposition FIR Wiener wavelet filter bank is applied to the noisy sinusoidal signal from previous chapters and the results are displayed in Figure 6.2. The normalized analysis filters for the Daubechies-20 wavelet used in this analysis are pictured in Figure 6.3. Figure 6.4 and 6.5 show the MSE performance for this filter bank with the Wiener filter order varying from 4 to 20. In Figure 6.4 an independent white Gaussian noise sample with a zero mean was provided to the Wiener filter while in Figure 6.5 the actual noise sample used to produce the noisy signal was provided. The MSE difference between the independent and actual noise cases becomes greater at lower SNR levels. A common feature of all wavelet based techniques is the improved MSE performance with successive levels of decomposition. This was observed in wavelet thresholding and will occur in the DR Wiener wavelet methods described in the next chapter. The Wiener filters are applied on a level dependent basis, which is different from universal and minimax wavelet thresholding methods where the threshold was applied in a global manner across all of the scales of a particular level of decomposition resulting in independent MSE curves. Another characteristic observed in wavelet based denoising methods is the degraded MSE performance in the transition regions of the analysis and synthesis filters. The DWT 61 filter bank, as originally constructed, met the perfect reconstruction conditions. Consequently, aliasing and distortion effects which would normally occur in the transition eee of the filters are canceled. Denoising by filtering or thresholding the DWT coefficients between the decimation and interpolation operations results in a filter bank which no longer meets the perfect reconstruction conditions. This effect is most noticeable in the filter transition regions because this is where the perfect reconstruction property was so essential in canceling the effects of aliasing and distortion. This unfortunate result can be minimized by choosing wavelet analysis and synthesis filters of longer length to reduce the transition region. Some detrimental effects occur as a result, including more computational effort and longer filter transient periods. The transient period, during which the filter output is unreliable, is equal to the filter length. As the filter length increases relative to the signal length the transient period may become significant. 62 Original (dashed) and Nicisy Signal: SNR ~~ odB, Signal Freq — Oo. Criginal and FiR Wiener Wavelet Fiitered Signal: MSEe= 0.06392 —a & i Ko) 156 325 ao as6c So rire ECsemiiie muciess Figure 6.2: FIR Wiener Wavelet Filter (a) Original (dashed) and Noisy Sinusoidal Signal. Frequency=0.1, Sampling Frequency=1. (B) Original (dashed) and filtered signal. Three level decomposition using Daubechies-20 wavelet. Filter Order = 8. Four Level Spectral Decomposition Using Daubechies—20 Wavelet Normalized Gain 0.2 0.25 0.3 0.35 o.4 Normalized Frequency Figure 6.3: Gain normalized analysis filters using Daubechies-20 wavelet. 63 FIR Wiener Wavelet: Order 4 FIR Wiener Wavelet: Order 8 0.25 Oo o 0.2 iT im g o o = 3.0.15 7) 7) Cc o c @ @ = = 0.1 0.05 0.2 0.3 Normalized Freq 0 0.1 0.2 0.3 04 05 Normalized Freq FIR Wiener Wavelet: Order 12 FIR Wiener Wavelet: Order 20 O25 Guu 5 , , aT mT @ Fs 5 3.0.15 5 oO Y) ep) c= o 7 ® © = 0.1 Ss 0.05 0 01) 02060306040 (05 0 O01 O02 O38 04 O85 Normalized Freq Normalized Freq Figure 6.4: FIR Wiener Wavelet filter applied to noisy sinusoidal signal: SNR = 0 dB. Four level DWT decomposition using Daubechies-20 wavelet. Independent noise sample. Ten trial average. FIR Wiener Wavelet: Order 4 FIR Wiener Wavelet: Order 8 0.4 O55 0.25 0.3 5 5 0.2 i 0-2 AT © © wo rn) = 3.0.15 ep) (op) (<= = © 0.1 eS = = 0.1 0.1 0.05 0.05 O 0 0.1 0.2 0.3 0.4 0.5 O 0.1 0.2 0.3 0.4 0.5 Normalized Freq Normalized Freq FIR Wiener Wavelet: Order 12 FIR Wiener Wavelet: Order 20 0.25 © No Mean Square Error © oi, on Mean Square Error 0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 0.3 0.4 0.5 Normalized Freq Normalized Freq Figure 6.5: FIR Wiener Wavelet filter applied to noisy sinusoidal signal: SNR = 0 dB. Four level DWT decomposition using Daubechies-20 wavelet. Actual noise sample. Ten trial average. 65 B. ALIAS CANCELLATION The conditions for perfect reconstruction were specified in Chapter 4. A two-channel filter bank was analyzed in the Z-domain with the following result X(@) = YF y(2Ho@) +F@H,@)X@)+{Fo@H(-2)+F @H,(-2)X(-z). (6.1) Inserting the Wiener filters to denoise the data in the wavelet domain alters (6.1) as follows: R@ = YF )We (2) +F (W202) X@) (6.2) + UF o(z)Wo(2)Hy(-2) +F (2)W,(2?)H,(-2)] X(-2). Alias cancellation is obtained by choosing the following synthesis filters, F(z) = W(z2)H,(-z o(Z) 1( ; 1(-Z) 6.3) Leila) a o(Z )Hy(-2Z). which cancel the X(-z) alias term in (6.2) leading to X@ = LW, (e?)Wo(e HOH, (-2)-Hy(-H (IX) (6.4) As before, the alternating flip relationship between the analysis filters is chosen to eliminate distortion affects: EGA) = — vee 2) (6.5) Substituting (6.5) into (6.4) results in the following Z-domain expression for the reconstructed signal in terms of the orginal signal, X(z) = W,(@*)W,(z7)z “X(2). (6.6) Equation (6.6) is applied to the two-level decomposition pictured in Figure 6.6 with the following result for the reconstructed Z-domain signal, X(z) =W,(z7)W,(z*)W,(24)z *X(2). (6.7) 66 Equations (6.6) and (6.7) can be generalized for a J level decomposition as follows, X(z) =W,(z7)z *X(z) 1] Wz). (6.8) The reconstructed signal is now a function of a delay determined by filter length N, and the product of a series of alias cancellation filters given in (6.8). These alias cancellation filters are nothing more than interpolated Wiener Filters which are multiplied in the Z-domain or convolved in the time domain. Unfortunately, simulations show that the alias cancellation filters produce significant phase and amplitude distortion which is more detrimental than the W,(D) Analysis Filters Wiener Filters Synthesis Filters Alias Cancellation Filters Figure 6.6: Two-level FIR Wiener Wavelet Filter Bank with Alias Cancellation aliasing effects they are designed to remove. The FIR Wiener Wavelet filter with alias cancellation is applied to the same noisy sinusoidal signal as used before. The denoised signal is poorly reconstructed based on an evaluation of the MSE alone. The MSE values are high with erratic results across the 67 spectrum. A comparison of the power spectral densities of the noisy signal and the denoised signal for two-level FIR Wiener wavelet ftlters with and without alias cancellation is shown in Figure 6.7. Figures 6.7 (a) and (b) show the power spectral densities of the noisy and denoised signal without alias cancellation. The alias cancellation filter with frequency response as shown in Figure 6.7 (c) is used to produce the denoised signal shown in Figure 6.7 (d) which has been normalized to remove the effects of the gain attenuation. Although the aliasing frequency at a normalized frequency of 0.2 has been cancelled, significant distortion is occuring at this same frequency as well as at normalized frequencies of 0.05 and 0.45. The average gain attenuation for a two-level FIR Wiener wavelet decomposition with alias cancellation (but without energy normalization of signal spectra) is 41.5 dB measured at the signal frequency which is quite significant when compared to 0.49 dB for the same filter without alias cancellation (Attenuation values for signal frequencies from 0.01 to 0.49 were computed and then averaged). The FIR Wiener filter with alias cancellation might still be considered a suitable method if the denoised signal can be normalized. A final comparison is made between the two filtering methods based on the difference of two spectra, determined from energy normalized signals, on a log magnitude versus frequency scale defined by: Viw) = log S(w) - log S‘(w). (6.9) where S(w)is the energy normalized noise-free signal spectra and S(w)is the energy normalized denoised signal spectra. The Mean Absolute Log Spectrum Distortion (MALSD) is defined in reference [33] as: 68 MALSD = [Ivey (6.10) Expressed in terms of a summation of eee terms, (6.10) becomes MALSD = = lV(sti/N)I. (6.11) =0 where N is the number of frequency samples taken between zero and the half sampling frequency, 7. When applied to two-level decomposition filter banks with and without alias cancellation, the following results were obtained (ten trial average): FIR Wiener Wavelet Filter MALSD Without Alias Cancellation 1703 With Alias Cancellation 2584 The FIR Wiener Wavelet filter without alias cancellation exhibits a smaller magnitude MALSD indicating the addition of alias cancellation filters is detrimental to the denoising performance of this filter. For this reason any further results obtained using a FIR Wiener Wavelet filter will not include alias cancellation filters. 69 (a) Noisy Signal: SNR = OdB (b) FIR Wien. Wave. No Alias Canc. Amplitude (dB) Amplitude (dB) Freq ) Freq (c) Alias Canc. Filter Freq. Response (d) FIR Wien. Wave. With Alias Canc. -—20 —25 ran) fon 13) = ® @ 3 -35 E ro ie: £ = < < —40 —45 —50 0 0.1 0.2 0.3 0.4 OLS: 8) 0.1 0.2 0.3 0.4 0.5 Freq Freq Figure 6.7: Power spectral densities of: (a) noisy sinusoidal signal, Frequency=0.3, F,=1, SNR=0 dB, (b) denoised result as determined with no alias cancellation filter applied, (d) denoised result with alias cancellation filter in (c) applied. 70 VIL. WR WIENER WAVELET FILTER The denoising methods developed thus far have been restricted to the FIR Wiener filter, Wavelet thresholding methods, and a combination of the two. This chapter develops the noncausal IR Wiener filter and then introduces its application in the wavelet domain. Next, the results are compared to previously discussed methods. A. IIR WIENER FILTER The IIR Wiener filter 1s recursive in form and requires fewer parameters to determine the optimal filter weights than a comparable FIR form of the filter [12] . The IR causal filter problem was solved (originally in the continuous case) by Wiener in the transform domain using spectral factorization methods [11]. This section describes the transform domain solution of the noncausal Wiener-Hopf equation. The noncausal solution allows the filter to “look ahead” of real time and thus operates in an off-line application. This advantage results in improved performance in the mean-square error sense over the causal FIR Wiener filter. The noncausal IR Wiener filter estimate of the desired signal is of the form §(n) =y(n) = )) A *(@)x(n-k), . (7.1) k=-@ which differs from (3.1) in two respects. First, the upper limit of the sum extends to +o and the impulse response has nonzero values for n <0. Application of the orthogonality principle yields the noncausal IR form of the Wiener-Hopf equation x R (i-k)h*(k) = RQ); -o <i <o (7.2) k=-c fl which can be written in convolutional form as: Ri) *h (i) = Rid), (7.3) with mean-square error o. = R{0)- y h *(k)R,(k). (7.4) k= -0 After conjugating equation (7.3), the left-hand side is a discrete convolution of the filter coefficients h(k), and the auto-correlation function, R,(k). Expressing the Wiener-Hopf equation in the frequency domain leads to: S(@)H(w) = S_(o), C765.) where H(w) = > h(nje", -n<ws<t Sw). = > R (ne OW ee <0) Te (7.6) S(@) = > R (ne (On ee Thus the filter transfer function H(w) can be expressed as: S@) A A@) = CO) S(@) Assuming that the signal and noise are independent and have a zero mean results in: S(@)=S (@) and S(w)=S (@)+S (@), (7.8) allowing (7.7) to be expressed as S(@) - S(@ Hf, (@) = Sa Ss (7.9) 5 (w) Oe which may be rewritten as S.(w)/S (w) Aw) = S@y5,(a) #1 (7.10) The expression S,(@)/S,(@) is a measure of the signal-to-noise power ratio at frequency w. When S,(@)/S,,(@) >> 1 the filter coefficient magnitude, H,.(@), approaches unity, and when S,(q)/S,,(@) << 1 the filter coefficient magnitude approaches zero. Between these two limits the filter coefficients are chosen to pass the signal and eliminate noise with minimal distortion [34]. B. IIR WIENER FILTER APPLIED TO THE WAVELET DOMAIN A key property which underlies the success of wavelets is that they form an unconditional basis for a wide variety of signal classes [35]. Wavelet expansions can therefore concentrate the signal energy into a relatively small number of large coefficients. This signal compaction property makes the discrete wavelet transform an ideal tool in constructing an empirical IIR Wiener filter. The result will be referred to as the IIR Wiener Wavelet filter. The original time-series model given by x(n) = s(n) + w(n); n = 1,2,...,.N (7.11) becomes, in the wavelet domain, y, Up Ce eee! 3 ka 1,2,..,N, (7.12) where the terms are the DWT coefficients found from the following x(n) = » Yo ak) s(n) = y 6,27 y(2t-k) (7.13) w(n) = y Zi 2)? w(2t-k). 1S The indices j and k refer to level of DWT decomposition and signal sample number respectively. The signal estimate § is determined by applying the inverse wavelet transform to the coefficients, Oy The noise used to produce the noisy signal is provided in order to determine the optimal Wiener filter coefficients. Denoising an actual ocean acoustic signal would require an isolated noise only segment or some other noise estimate in order to determine the approprite noise statistical characteristics. A noise sample or estimate which is not highly correlated with the actual noise will degrade the denoising capability of this filter. The observation sequence wavelet coefficients, y,, , and noise wavelet coefficients, Zj,» are provided allowing (7.9) to be rewritten as follows for a multiple level DWT: > y2 - S23 N (7.14) h, = a aa eel ae i ee | m ; Dui DSH k=] where WN is the signal length, 7 is the level of the DWT decomposition, and J represents the maximum decomposition level. This filter is applied to the observation sequence wavelet coefficients at each decomposition level to determine the signal estimate coefficients as follows See 10 re alse eat k=1,2,...,N2 7. (Ties 74 The IIR Wiener Wavelet method is implemented in a one-level filter bank structure, as illustrated in Figure 7.1. This is the same filter bank structure as that used for wavelet thresholding and FIR Wiener wavelet filtering except that now at each scale the noisy wavelet coefficients are thresholded by the filter coefficients h, and h,. The structure may be expanded to additional levels or scales, limited only by the signal length. This denoising method is applied to the noisy sinusoid used in previous methods with results as shown if Figure 7.2. The filter represents the denoised signal as a smoothed sinusoid after a three level decomposition using the Daubechies-20 wavelet. The Daubechies-20 analysis filters Input Analysis | Decimators IR Wiener Filters Interpolators Synthesis Output Figure 7.1: One-level noncausal IR Wiener wavelet filter bank. of length forty are pictured in Figure 7.3 with gain normalized to the maximum gain filter 12 in order to show the spectral bands associated with the scales. This normalization is only used to allow a better spectral view of the filters and is not representative of the actual filter gains used in the DWT. Original (dashed) and Noisy Signal: SNR = OdB, Signal Freq = 0.1 Amplitude Amplitude oO 5 10 15 20 25 30 35 40 45 50 Time (Sample Number) Figure 7.2: IR Wiener Wavelet Filter (a) Original (dashed) and Noisy Sinusoidal Signal Frequency=0.1, F=1. (b) Original (dashed) and filtered signal. Three level decomposition using Daubechies-20 wavelet. A ten trial average normalized MSE is plotted across the normalized frequency range of 0 to 0.5 for a sampling frequency of 1 as shown in Figures 7.4 and 7.5. An independent noise sample of white Gaussian noise with a zero mean was used in Figure 7.4 while the actual noise sample was used in Figure 7.5. The MSE is reduced with each successive scale or decomposition level. This occurs because with each successive scale there is less noise energy as represented by the wavelet coefficients, whereas the signal energy remains the same. The spectral bandwidth is reduced with each level of decomposition eliminating approximately half the remaining white noise which reduces the noise variance and improves 76 the SNR at that particular scale. This leads to an IR Wiener filter which is more effective at removing noise as demonstrated by the improved normalized MSE. The IR Wiener filter is a level dependent filter which means that increasing the level of decomposition will only improve the MSE performance for the spectral band associated with the particular analysis filter. This concept is better understood by comparing Figures 7.3 and 7.4. Notice the degraded MSE that occurs in the analysis filter overlap regions. The spectral range over which this degradation occurs can be reduced by choosing analysis filters with smaller transition regions. Higher order filters will accomplish this at the expense of greater computational effort. The distortion can never be completely removed since the perfect reconstruction property of the filter bank has not been met due to the insertion of the IIR Wiener filter. Four Level Spectral Decomposition Using Daubechies—20 Wavelet Normalized Gain O 0.05 ‘ 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 Normalized Frequency Figure 7.3: Normalized Spectral Decomposition Filters Used For IIR Wiener Wavelet Denoising. (Daubechies-20 Wavelet). 77 MR Wiener Wavelet Filter: SNR @w CO daB Mean Square Encor o.35 o.4 oO 0.08 o.7 o.18 oC.2 0.25 eo. Normalized Freq Figure 7.4: Mean-square error comparison for four level decomposition of noisy sinusoid using Daubechies-20 wavelet. Independent noise sample. Average of 10 trials. uR Wiener Vaveilet Filter Reautlta: S~rRrm = OdBB o.1& o.2 o.265 4 NMorwnalizedcdt Freq Figure 7.5: Mean-square error comparison for four level decomposition of noisy sinusoid using Daubechies-20 wavelet. Actual noise sample. Average of 10 trials. 78 C3 WIENERSHRINK The IR Wiener wavelet filter just discussed requires an isolated noise segment in order to determine the statistical noise characteristics necessary to generate a Wiener filter. WienerShrink, developed by Ghael, Sayeed, and Baraniuk of Rice University [36], is a wavelet-based Wiener filtering technique which differs from previous methods in this regard. Initially, a wavelet thresholding technique similar to those discussed 1n Chapter V 1s applied to the noisy signal. The resulting denoised signal, referred to as the pilot signal, is then wavelet transformed a second time and used to construct an IR Wiener filter. This wavelet- based empirical Wiener filtering method is described by the block diagram shown in Figure MO. The noisy signal as represented by (7.11) 1s wavelet transformed by the DWT Figure 7.6: Wavelet-based empirical Wiener filtering [25]. 79 represented by W,. Any of the four wavelet thresholding methods, represented by block H, is applied to the wavelet coefficients y” ,,, to produce the wavelet coefficient signal estimate, Sines: This is then transformed by the inverse DWT represented by W,”' to form the signal estimate, 5,. The signal estimate is then wavelet transformed again in block W, to form the wavelet coefficient signal estimate cy This 1s used to design an empirical IIR Wiener filter given by Aan) >, (Or 5)” ” bh. =k ne= —e a 202 (7.16) ays Soul Vik k=1 where N is the signal length, 7 is the level of the DWT decomposition, and J represents the maximum decomposition level. This empirical IR Wiener filter is the ratio of the energy of the estimated noise-free signal wavelet coefficients for a particular scale divided by the energy of the noisy signal wavelet coefficients for the same scale. The IIR Wiener filter is applied to a DWT of the original noisy signal y”,,, formed by W;, to produce the wavelet coefficient signal estimate ome An inverse DWT, represented by W,’ is applied to produce the final signal estimate §. A Daubechies-8 wavelet is used for the DWT in block W, while a Daubechies-16 wavelet is used for the DWT in blocks W, and W,. Ghael et al provide the following explanation for the improved performance of their product compared to traditional thresholding methods [36]. The goal of any wavelet based denoising technique is to determine the noise-free signal wavelet coefficients from the noisy signal wavelet coefficients. The noise-free wavelet coefficients consist of a number of 80 trustworthy coefficients which are most assuredly due to the signal and then a number of less trustworthy coefficients which are more difficult to determine due to the presence of noise coefficients. WienerShrink uses the initial thresholding method to determine the trustworthy coefficients given by ee The untrustworthy or dubious coefficients are then predicted by the IIR Wiener filter formed from the trustworthy coefficients. Ghael et al indicate that classical wavelet thresholding methods are overly conservative in thresholding the coefficients resulting in the retention of only the trustworthy coefficients. WienerShrink is applied to the same noisy sinusoidal signal used previously with the MSE results displayed in Figure 7.7 for a SNR of 0 dB and Figure 7.8 fora SNR of 10 dB. The Heuristic SURE and SURE methods of thresholding were chosen due to their superior performance over minimax and universal thresholding methods at lower SNR levels. The hard threshold option outperforms the soft threshold option in some areas, especially at the lower SNR and the lowest scale. The difference is not substantial allowing either thresholding option to be utilized. A comparison of Figures 7.4 and 7.7 reveals that the minimum MSE obtained by both the IIR Wiener wavelet and WienerShmnnk at a particular scale are nearly identical. The IZR Wiener wavelet exhibits superior MSE performance in the finer scales for a particular decomposition level. For example, the single-level or scale decomposition for IR Wiener wavelet produces an average MSE value of 0.37 in the lowest scale (normalized frequency of 0.25 to 0.5) while the best WienerShrink method for a single- level or scale decomposition with Heuristic SURE hard thresholding achieves an average MSE value of only 0.6 at the same scale. 81 If the filter bank can be adapted by changing the level of decomposition (alternatively, the number of scales) to confine the signal of interest to the highest scale then this difference would have little effect. Otherwise, the IR Wiener wavelet filter may be preferable to WienerShrink. Thus the choice becomes a signal dependent decision. A significant improvement in MSE performance for Wienershrink is achieved when the SNR is increased to 10 dB (see Figure 7.8). This difference is particularly noticeable for denoising methods which employ wavelet thresholding. Comparing wavelet thresholding methods (Figure 5.9) to WienerShrink (Figure 7.7) reveals minimal difference between the two when applied to a noisy sinusoidal signal. WienerShrink is slightly better in terms of MSE performance. In Chapter VIII signals with more interesting characteristics will be denoised using the various methods and the results compared. 82 Heuristic SURE Soft Heuristic SURE Hard Mean Square Error 0.2 0.3 . , 03 04 «0.5 Normalized Freq a Normalized Freq SURE Soft SURE Hard Mean Square Error 0 0 ee C203 04 was 0 0.1 0.2 0.3 0.4 0.5 Normalized Freq Normalized Freq Figure 7.7: WienerShrink applied to noisy sinusoidal signal: SNR = 0 dB. Four level DWT decomposition using Daubechies-8 to form pilot signal and Daubechies-16 for IIR filtering. Heuristic SURE and SURE thresholds used with soft and hard options. 83 Heuristic SURE Soft Heuristic SURE Hard 0.12 0.1 Mean Square Error ©o So © © oO 10.6) ad oO BR ° re) © 0.2 0.3 : : ; 0.2 0.3 Normalized Freq Normalized Freq SURE Soft | SURE Hard OZ Mean Square Error od © © © Oo) 0O © © rs 0.02 ) 0.1 O02 03 O04 #405 0 01 02, 03° 0455 Normalized Freq Normalized Freq Figure 7.8: WienerShrink applied to noisy sinusoidal signal: SNR = 10 dB. Four level DWT decomposition using Daubechies-8 to form pilot signal and Daubechies-16 for IR filtering. Heuristic SURE and SURE thresholds used with soft and hard options. 84 VIII. COMPARISON OF DENOISING METHODS Various Wiener and wavelet filtering methods have been presented. Thus far the ability to denoise a noisy sinusoidal signal has been the basis for comparison between the denoising methods. In this chapter the denoising methods are applied to on synthetic signals which are commonly used as benchmarks for comparison because of their characteristics. These signals were originally chosen by the producers of Wavelab at the Stanford University Statistics Department [37]. The noise-free signals are shown below in Figure 8.1 and are referred to as Doppler, HeaviSine, Bumps, and Blocks. Doppler HeaviSine 0 900 1000 "0 900 1000 Bumps Blocks —20 0 900 1000 0 900 1000 Figure 8.1: Noise-free test signals: Doppler, HeaviSine, Bumps, and Blocks. 85 All of the denoising methods discussed in this thesis are applied to the four test signals and the average normalized MSE values are determined for twenty trials. The results are tabulated in the Appendix. Graphs depicting the MSE performance of the various methods are displayed as Figures 8.2 to 8.8. All methods, except for the Wiener filter, are compared in Figures 8.9 and 8.10. All Figures depicting the performance of Wiener wavelet denoising methods show two different assumptions concerning the noise estimation. Figures 8.4, 8.6, and 8.9 assume an ideal situation in which the exact noise sample that produced the noisy signal is supplied to the Wiener filter. In Figures 8.5, 8.7, and 8.10, an independent white Gaussian noise source with standard deviation equal to one and a zero mean is provided to the Wiener filter. The second option is realistic and better approximates the true performance of these methods. A visual comparison of four denoising methods is depicted in Figures 8.11 and 8.12. Noise-free and noisy (SNR=0 dB) doppler signals of length 1024 are shown in Figure 8.11. The noisy doppler signal has a MSE of 0.5774. Denoised doppler results are shown for wavelet thresholding, IIR Wiener wavelet, WaveShrink, and FIR Wiener wavelet noise removal techniques. These plots provide a visual interpretation of the MSE results. Specific observations regarding the denoising methods follow. A. WIENER FILTER The Wiener filter did not denoise these signals as well as the other methods except in the case of the HeaviSine signal. Doppler, Bumps, and Blocks are nonstationary and therefore are not effectively filtered by a standard Wiener filter. HeaviSine exhibits some 86 degree of stationarity, resulting in better denoising by the Wiener filter. A short-time Wiener filter could be designed with an optimal window size and filter length which would improve the Wiener filter MSE performance significantly. This type of filter 1s utilized in denoising nonstationary signals in the next chapter. B. WAVELET THRESHOLDING Wavelet thresholding or “Wave Shrink” performs as well as or better than all the other methods at SNR levels from -5 to 20 dB. In general the SURE and Hybrid thresholding methods are the best of the four thresholding methods. The appendix tables show results for both soft and hard threshold options, although only the soft option was used in the plots. In some cases better MSE results are obtained using the hard threshold, however, the soft threshold provided better overall results. The signals in this study have all been scaled in amplitude such that the desired SNR is achieved while maintaining the noise level at a standard deviation of one. Recall the method for determining the noise standard deviation in actual signals required determining the mean absolute deviation from the wavelet coefficients at the lowest or finest scale. In cases where some of the signal energy resides in this finest scale, the performance of the wavelet threshold denoising method will be degraded. c. FIR WIENER WAVELET FILTER The FIR Wiener wavelet filter is the best denoising method in terms of MSE performance in the ideal noise conditions described above and shown in Figures 8.4 and 8.9. Figures 8.5 and 8.10 show that this method is affected the most when provided with an 87 independent noise source. The MSE more than doubles for some of the test signals. An accurate noise estimate is critical to the performance of the FIR Wiener wavelet denoising method. Transient affects also must be considered when using a FIR filter. These shortcomings make this method the least preferred. D. IIR WIENER WAVELET FILTER The IR Wiener wavelet filter exhibits a MSE performance that is comparable to wavelet thresholding and WienerShrink. Although the MSE error performance is degraded when assuming an independent noise estimate as seen by comparing Figures 8.6 and 8.7, it still provides substantial denoising capability. Furthermore, it is the easiest to implement, requires the least computational effort, and has no transient effects. E. WIENERSHRINK WienerShrink outperformed all other methods but not by a substantial amount. It is likely that Ghael et al, the developers of this technique, may have applied an unpublished Statistical thresholding method in their original work which was not replicated in this thesis. 88 Doppler HeaviSine —5 0 5 10 15 —5 0 s) 10 15 SNR (dB) SNR (dB) Bumps Blocks 25 0 5 10 15 —5 0 5 10 15 SNR (dB) SNR (dB) Figure 8.2: MSE vs. SNR performance for Wiener Filter (Filter Order 8,12,16,20) on four test signals. Twenty trial average. 89 Doppler HeaviSine 0.25 0.05 £0 0 sa “ -5 0 S 10 15 -5 0 5 10 75 SNR (dB) SNR (dB) Bumps Blocks =5 0 5 10 15 -5 0 5 10 15 SNR (dB) SNR (dB) Figure 8.3: MSE vs. SNR performance for wavelet threshold denoising of four test signals. Threshold types shown include SURE, Heuristic SURE, Minimax, and Universal. Five level or scale decomposition using Daubechies-20 wavelet. Twenty trial average. 90 Doppler MSE 0 -5 0 5 10 15 SNR (dB) Bumps MSE o -5 0 5 Sie ets SNR (dB) HeaviSine 0.06 0.0 -5 0 5 10 15 SNR (dB) Blocks —5 0 5 10 15 SNR (dB) Figure 8.4: MSE vs. SNR performance for FIR Wiener wavelet denoising of four test signals. Wiener filter orders of 8 and 16 shown. Five level or scale decomposition using Daubechies-20 wavelet. Twenty trial average. Doppler HeaviSine ae 0 Sano Nani “5 0 5 10 15 SNR (aB) SNR (dB) Bumps Blocks 0.25 0.1 0.05 —5 0 5 10 15 -5 0 5 10 15 SNR (dB) SNR (dB) Figure 8.5: MSE vs. SNR performance for FIR Wiener wavelet denoising of four test signals. Wiener filter provided independent noise source. Filter orders of 8 and 16 shown. Five level or scale decomposition using Daubechies-20 wavelet. Twenty trial average. 22 Doppler HeaviSine 0.25 0.1 0.05 O Ra 2 ae 5 0 5 10 15 5 0 5 10 15 SNR (dB) | SNR (dB) Bumps Blocks 0.2 0.18 0.2 0.16} 0.14 0.2 O42 oD oD 2 0.15 2 0.1 0.08 0.1 0.06 04 0.05 oo 0.02 0 Ss 0 5 0 5 10 15 —5 0 5 10 15 SNR (dB) SNR (dB) Figure 8.6: MSE vs. SNR performance for IIR Wiener wavelet denoising of four test signals. Three, four, and five level or scale decomposition shown using Daubechies-20 wavelet. Twenty trial average. 93 Doppler HeaviSine MSE MSE 0 -5 0 5 10 15 SNR (dB) Bumps Blocks 0.2 0.18 0.25 0.16 0.14 0.2 0.1 QD D S 0.15 S 0.1 0.08 0.1 0.06 0.04 0.05 0.02 0 ————— 0 —5 0 5 10 15 -5 0 5 10 15 SNR (dB) SNR (dB) Figure 8.7: MSE vs. SNR performance for IIR Wiener wavelet denoising of four test signals. Wiener filter provided with independent noise source. Three, four, and five level or scale decomposition shown using Daubechies-20 wavelet. Twenty trial average. 94 Doppler HeaviSine MSE -5 0 5 10 15 -5 0 5 10 15 SNR (dB) SNR (dB) Bumps Blocks MSE -5 0 5 10 15 ) -5 0 5 10 15 SNR (dB) SNR (dB) Figure 8.8: MSE vs. SNR performance for WienerShrink denoising of four test signals. SURE, Heuristic SURE, Minimax, and Universal thresholding used with soft option. Five level or scale decomposition with Daubechies-20 wavelet. Twenty trial average. 5) Doppler HeaviSine —5 0 5 10 15 -5 0 5 10 15 SNR (dB) SNR (dB) Blocks 0 i ~ -5 0 5 10 «15 -5 0 5 10 15 SNR (dB) SNR (dB) Figure 8.9: MSE vs. SNR performance comparison of wavelet thresholding, WienerShrink, FIR Wiener wavelet, and IIR Wiener wavelet denoising of four test signals. 96 Doppler HeaviSine —5 0 5, 10 15 —5 ) 5 10 15 SNR (dB) SNR (dB) Blocks sso |S. 10mm. = 0 Ss MSCS SNR (dB) SNR (dB) Figure 8.10: MSE vs. SNR performance comparison of wavelet thresholding, WienerShrink, FIR Wiener wavelet, and IIR Wiener wavelet denoising of four test signals. Wiener filters provided with independent noise source. 97 Noise—Free Doppler 0 100 200 300 400 500 600 700 800 900 1000 Noisy Doppler: SNR= 0 dB, MSE = 0.5774 0) 100 200 300 400 500 600 700 800 900 1000 Figure 8.11: Noise-free and noisy doppler signal of length 1024 data points. SNR=0 dB. MSE of noisy doppler is 0.5774. 98 Wave. Thresh.: MSE=0.08859 IIR Wien/Wave: MSE=0.09048 2.5 2 eS 1 0.5 O —0.5 = -1.5 —2 Be" 200 400 600 800 1000 70 200 400 600 800 1000 WaveShrink: MSE=0.1021 FIR Wien/Wave: MSE=0.1449 2.5 2.5 2 2 ES) AS 1 1 0.5 0.5 O 0 | -0.5 —0.5 = 4 —1.5 —1.5 —2 -2 Fo 200 400 600 800 1000 -—< 200 400 600 800 1000 Figure 8.12: Comparison of four denoising techniques applied to noisy doppler with SNR=0 dB. MSE values are shown for each method. Independent noise sources provided to Wiener filters. Daubechies-20 wavelet used in five-level or scale decompositions. SURE (soft) threshold used for wavelet thresholding. be, IX. WAVELET PACKETS Wavelet analysis has been shown to be particularly effective at removing noise from low frequency signals such as narrowband tonals which might occur in the ocean environment. The constant-Q filters provide a logarithmic frequency resolution with narrow bandwidths at low frequencies and wide bandwidths at high frequencies. This high frequency resolution at low frequencies is produced by iterating the lowpass scaling branch of the Mallat algorithm tree. In addition to low frequency tonals, ocean acoustic signals often include high frequency signals which may be of interest as well. These signals are usually short duration transients or acoustic energy pulses. Wavelet based denoising techniques are not optimal for removing noise from this class of signals. A richer signal analysis tool, termed wavelet packets, was introduced by Ronald Coifman [38] to provide high resolution decomposition at high frequencies. This 1s accomplished by iterating the highpass wavelet branch of the Mallat algorithm tree. The detail coefficient vector is thus decomposed into two parts by splitting, filtering and decimating in the same manner as the approximation coefficient vector. The full binary tree is pictured in Figure 9.1 where the H and L blocks are the high and low-pass filters determined by the wavelet and scaling functions. The numerous signal expansions that are possible with wavelet packets come at a cost in computational complexity of O(N log,.()) compared to O(N) for the wavelet transform. 101 Hl W530 W354 W3. W333 W334 W355 W356 W37 Figure 9.1: Wavelet packet tree (Used with permission from Joshua Altmann, [14]). A. WAVELET PACKET THEORY Wavelet packets are formed by modifying (4.16 ) and (4.17 ), which determine the DWT filter coefficients in terms of the scaling and wavelet functions. A third index is required to describe the wavelet packet in terms of its position within the wavelet packet tree. This index describes the wavelet packet in terms of its frequency resolution. The wavelet packet equations are given by o,(2) (0) h(n) ¥2.,(2t-n), neZ 9.1 » a(n) ¥2.,(2t-n), neé€ Z. a The wavelet packets for a three-level wavelet packet decomposition using the familiar Daubechies-2 wavelet are shown in Figure 9.2. These wavelet packets correspond to the Figure 9.1 tree with the notation W, ,, with j = 3 denoting the scale parameter and p =0 to 7 the frequency parameter. Notice the first two wavelet packets correspond to the Daubechies-2 filters (Figure 4.11) developed from the scaling and wavelet functions for the 102 0 w2 W3 3 3 3 3 2 2 2 2 1 1 { 1 0 0 0 0 = ~ 4 | = 2, a: - xD Figure 9.2: Daubechies-2 wavelet packets for a three-level or scale decomposition. W, , for j = 3 and p = 0 to 7. Supported on the interval [0,3]. © ol, NO © ak NO © =i NO 3 3 3 3 2 2 2 2 1 1 1 1 0 0 0 0 i = a at &D 2 2 = ee wo 1 2 W4 = ol = Oo = ~] wavelet transform. Thus the wavelet transform is a subset of wavelet packets. The wavelet packet decomposition of a signal results in a binary tree composed of waveforms known as “atoms.” This collection of atoms forms an overcomplete “dictionary” or library from which a signal of length NV can be decomposed and reconstructed in at most 2” different ways. The Best Basis algorithm, developed by R. Coifman and V. Wickerhauser [38], and used in MATLAB [22], seeks to minimize an additive cost function in order to determine the best basis for representing the signal. The cost function chosen is 103 entropy given by E(s) = -}/ 5; log(s;) (9.3) where s is the signal and s, the coefficients of s in an orthonormal basis. Best Basis is an efficient algorithm which computes the entropy at each node in the binary tree and then finds the best wavelet packet representation of the signal based on minimizing this quantity. Reference [34] provides additional information. B. WAVELET PACKETS APPLIED TO TEST SIGNAL The test signal chosen for analysis 1s a series of three high frequency transient pulses of decreasing amplitude. This test signal, which will be referred to as transients, was used by Barsanti in reference [8]. The signal is produced from the following truncated exponentially decaying sinusoids: Sl = sin(21k1/4) xexp(-k1/200) Kigee2....,256 $2 = 1/2 sin(27k2/4) xexp(-k2/200) k2 = 1,2,...,200 (9.2) S3 = 1/3 sin(27k3/4) *exp(-k3/200) KS 2 clo This is a realistic acoustic signal which 1s artificially produced to allow MSE comparisons at differing SNR levels. The noisy version of transients is produced by adding white Gaussian noise with a zero mean. The standard deviation of the noise remains constant at a value of one while the noise-free signal is scaled to obtain the desired SNR level. Noise- free transients and noisy transients with an SNR of 0 dB are shown in Figure 9.3. The first transient is also pictured expanded in the noise-free and noisy form so that the details may be better observed. 104 Noise—free Transients Single Noise—free Transient 10 10 8 8 6 6 4 4 2 2 O 0 22 -—2 4 24 -6 6 -8 -8 ms 5000 10000 15000 ~105 100 200 300 Noisy Transients: SNR = 0 dB 73 Single Noisy Transient 8 6 4 2 O : -2 —4 —6 } 3 -6 -8 | -8 a: 5000 10000 15000 mS 100 200 300 Figure 9.3: Noise-free and noisy transients (SNR = 0 dB). Expanded view of first transient without and with noise (SNR = 0 dB). Seven different methods of denoising are applied to the transients. The results for the short-time Wiener filter is shown in Figure 9.4. The other six methods are compared in Figures 9.5 through 9.10 for SNR levels of 5, 0, and -5 dB. Specific comments concerning each method follow. 105 1. Short-time Wiener Filter The short-time Wiener filter of length twenty is applied to the 16384 point transient signal in segments of 256 points. The noise source provided to the Wiener filter is an independent noise source from that used to produce the noisy transients signal. It is white Gaussian noise with standard deviation equal to one and a zero mean. Windowing the data allows the Wiener filter to completely remove the noise in noise-only segments. This does not occur when the Wiener filter 1s applied to the entire signal in one block due to the non- stationary nature of the data. ae Wavelet Thresholding This signal does not meet the low frequency assumption of wavelet denoising and therefore the MSE is greater than that achieved by wavelet packet thresholding. When using wavelet thresholding for denoising, the approximation coefficients are not usually thresholded since the assumption is that these coefficients contain the signal for the most part. The algorithm was modified to support approximation coefficient thresholding which removed additional noise and reduced the MSE. 3. FIR Wiener Wavelet Filter The FIR Wiener wavelet filter does not perform well for this signal, again due to the high frequency characteristics of the signal. The MSE is higher than what could be optimally achieved due to the residual noise present in the noise-only segments. This is true for all of the Wiener/Wavelet methods. This is partially remedied by implementing a windowed or segmented filtering algorithm that uses a triangular window with fifty- percent overlap. A 106 window length of 4096 points with a Wiener filter length of four and decomposition level of two achieved the best results which are listed in the table at the end of this section. The short-time FIR Wiener wavelet filter approaches the MSE performance of the FIR Wiener Wavelet packet filter. 4. IIR Wiener Wavelet Filter The observations made for the FIR Wiener wavelet filter apply to this method as well. The IIR Wiener/Wavelet filter experiences the most severe degradation when the signal to be denoised doesn’t meet the low frequency assumption. Significant improvement is achieved by implementing the same triangular window short-time filtering algorithm described above for the short-time FIR Wiener wavelet filter. A window length of 256 points was used to produce the results shown in Figure 9.11. Additional improvement in terms of MSE is realized by applying a threshold to the WR Wiener wavelet filter coefficients. If the filter coefficient magnitude is below some threshold value, in this case 0.4 was chosen by trial and error, then the coefficient is set to zero. The resulting improvement is seen in column four of Figure 9.11. 5. Wavelet Packet Thresholding Wavelet Packet thresholding exhibits significant MSE improvement over wavelet thresholding. Greater frequency resolution at high frequencies allows this method to produce transient pulses with better fidelity. Wavelet packet thresholding uses the SURE threshold described previously in Chapter IV. The approximation coefficients, represented by the “atom” in the lower left corner of the binary tree, are thresholded as well, leading to a noise- 107 free result in the segments between the pulses which contained only noise. 6. FIR Wiener Wavelet Packet Filter In this method a separate FIR Wiener filter is applied to each atom in the terminal or lowest nodes of the binary tree. No Best Basis optimization is performed. Some MSE improvement is realized when comparing this method to the FIR Wiener wavelet method but, the short-time Wiener filter and wavelet packet thresholding still achieve better results. De IIR Wiener Wavelet Packet Filter In this method an IIR Wiener Wavelet filter is applied to each atom at the terminal nodes. The IIR Wiener Wavelet Packet filter displays some improvement over the IR Wiener wavelet filter but it is the poorest performer of all the wavelet packet methods. It would likely achieve better results in some type of windowed or segmented implementation. 8. Summary The table following shows a comparison of all the methods in increasing MSE order. An optimized ST Wiener filter still outperforms all other methods with the wavelet packet threshold method and the ST IIR Wiener/wavelet method following close behind. Methods which do not utilize wavelet packets or a short-time filter are unable to effectively denoise this high frequency transient signal. For this particular signal, a short-time implementation of wavelet packet based denoising methods produced no further improvement and was computationally cumbersome. 108 Method ts dB | op | sap Swine nae | oe | mes | Wavelet Threshold 0.8283 0.3832 0.1399 TR Wiener/wavelet 0.8555 0.4788 109 Short-time Wiener (SNR=—-5dB): MSE=.1468 0 2000 4000 6000 8000 10000 12000 14000 16000 Short-time Wiener (SNR=0dB): MSE=.0684 10 2000 4000 6000 8000 10000 12000 14000 16000 Short-time Wiener (SNR=5dB): MSE=.0363 -5 -10 -15 0 2000 4000 6000 8000 10000 12000 14000 16000 Figure 9.4: Denoising transients using short-time Wiener filter of length twenty. Window size: 256. SNR levels of -5, 0, 5 dB. 110 Noisy Transients: SNR=5dB 15 10 0 5000 10000 15000 FIR Wien/Wave: MSE=0.1768 15 10 -15 0 5000 10000 15000 Wavelet Thresh: MSE=0.1399 15 10 0 3000 10000 15000 IIR Wien/Wave: MSE=0.4788 15 10 -5 -10 -15 0 5000 10000 15000 Figure 9.5: Denoising transients using wavelet thresholding, FIR Wiener/Wavelet, and UR Wiener/Wavelet methods. SNR =5 dB. Noisy Transients: SNR=5dB 0 5000 10000 15000 WP FIR Wien/Wave: MSE=.0969 15 10 -15 0 9000 10000 15000 WP Thresh: MSE=.0439 15 10 0 5000 10000 15000 WP IIR Wien/Wave: MSE=.1283 15 10 —15 0 5000 10000 15000 Figure 9.6: Denoising transients using wavelet packet thresholding, FIR Wiener/Wavelet packet, and IR Wiener/Wavelet packet methods. SNR = 5 dB. Noisy Transients: SNR=0dB Wavelet Thresh: MSE=0.3832 10 10 0 5000 10000 15000 0 5000 10000 15000 FIR Wien/Wave: MSE=0.3558 IIR Wien/Wave: MSE=0.8555 10 10 8 8 6 6 4 4 Z 2 0 0 2 -2 fl =A 6 —6 -8 -8 m6 5000 10000 15000 6 5000 10000 15000 Figure 9.7: Denoising transients using wavelet thresholding, FIR Wiener/Wavelet, and IIR Wiener/Wavelet methods. SNR = 0 db. Ws Noisy Transients: SNR=0dB WP Thresh: MSE=.1231 10 0 5000 10000 15000 0 5000 10000 15000 WP FIR Wien/Wave: MSE=.1884 WP JiR Wien/Wave: MSE=.3029 10 10 -10 ~10 0 5000 10000 15000 0 5000 10000 15000 Figure 9.8: Denoising transients using wavelet packet thresholding, FIR Wiener/Wavelet packet, and IR Wiener/Wavelet packet methods. SNR = 0 dB. 114 Noisy Transients: SNR=-5dB Wavelet Thresh: MSE=0.8283 =¥) | —2 -4 4 6 5000 10000 15000 0 5000 10000 15000 p FIR Wien/Wave: MSE=0.6944 d IIR Wien/Wave: MSE=1.2640 4 4 2 2 0 0 2) —2 -4 4 sais 5000 10000 ~—-15000 "6 5000 10000 15000 Figure 9.9: Denoising transients using wavelet thresholding, FIR Wiener/Wavelet, and IR Wiener/Wavelet methods. SNR =- 5 dB. els Noisy Transients: SNR=—5dB WP Thresh: MSE=.2888 6 4 2 0 -2 =4 0 5000 10000 15000 “8 5000 10000 15000 ; WP FIR Wien/Wave: MSE=.4102 WP IIR Wien/Wave: MSE=.7055 4 2 0 =o -2 ~4 4 <0 5000 10000 15000 =6 5000 10000 15000 Figure 9.10: Denoising transients using wavelet packet thresholding, FIR Wiener/Wavelet packet, and IIR Wiener/Wavelet packet methods. SNR = -5 dB. 116 ST IIR: MSE=0.07569 Transients: SNR=5dB ro 10 5 5 , O 26 -5 “195 5000 10000 15000 ~ 0 5000 10000 15000 Transients: SNR=0dB ST IIR: MSE=0.2025 6 6 4} 4 2 2 0! O —2'7 —2 =4 = —6 —6 -8 “=e 5 -8 @) 5000 10000 15000 0 5000 10000 15000 Transients: SNR=—5dB ST IIR: MSE=0.5950 6 4 2 =6 6 0 5000 10000 15000 0 5000 10000 15000 ST IIR w/thr: MSE=0.06108 10 —10 O 5000 10000 15000 ST IIR w/éAthr: MSE=0.1562 Oo NM Ff ODO —6 —8 0 5000 10000 15000 ST 1IR w/thr: MSE=0.3751 6 4 2 0 —2 ~4 —6 ) 5000 10000 15000 Figure 9.11: Short-time IIR Wiener wavelet filter applied to transients with SNR levels of 5,0, and -5 dB. 1 column: noisy transients. 2" column: denoised transients. 3" column: denoised transients using ST IIR Wiener wavelet filter with hard threshold. 117 X. CONCLUSIONS This study was undertaken to determine the feasibility of combining Wiener filter and wavelet based denoising techniques. The two methods have been used independently for years to denoise many types of signals. Their development and performance were demonstrated in Chapters III through V. Before any attempt was made to combine the methods, legitimate issues arose concerning the aliasing effects produced by inserting a Wiener filter into the DWT filterbank and how this might affect the ability to perfectly reconstruct the denoised signal. This concern led to the development of an alias cancellation filter in Chapter VI. It was hoped that this would cancel the deleterious affects produced by aliasing. Unfortunately, the alias cancellation filter only produced greater distortion and was not found to be of any benefit. Though the original intent was to combine a FIR Wiener filter with wavelet analysis, the possibility of an IIR Wiener wavelet filter was motivated by the WienerShrink method of reference [36]. The simplicity of this method combined with its advertised performance led to the IR Wiener wavelet filter methods developed in Chapter VI. In the earlier chapters the Wiener filter, wavelet thresholding, the FIR Wiener wavelt filter and IR Wiener wavelet filters were developed and applied to a stationary noisy sinusoidal signal. Initial results indicated that the IIR and FIR Wiener wavelet filters might actually outperform, in a MSE sense, the more traditional Wiener filter and wavelet thresholding methods. These results proved that the Wiener wavelet combination could indeed denoise a signal. However, removing the noise from a single sinusoidal waveform 119 was not deemed to be a comprehensive test of the denoising capabilities of these techniques. Therefore, the robustness of the various methods was measured by applying them to four test signals. In these results, the IR Wiener wavelet methods compared favorable with wavelet thresholding. The FIR Wiener wavelet filter results which seemed promising in Chapter VI were rather disappointing compared to the other methods with MSE values as much as twice those of the other methods. Since many ocean acoustic signals of interest are high frequency transients, wavelet packet methods were developed in Chapter [IX to provide greater frequency resolution at high frequencies. The FIR and OR Wiener filtering techniques were applied to the terminal nodes of the wavelet packet decomposition oe providing an enhanced capability in denoising high frequency signals compared to the previously developed Wiener wavelet methods. The various denoising methods were applied to synthetically generated noisy transients. The results were similar in that the traditional methods of short-time Wiener filtering and wavelet packet thresholding outperformed all other methods. However, an IR Wiener wavelet filter used in a short-time implementation performed rather well with MSE values somewhat greater than the traditional methods. This study has proven the feasibility of combining Wiener and wavelet based techniques into a single filter. In general, it does not outperform either of the methods from which it is derived. The basis for comparison between the various methods is not by any means clear cut. Each method has numerous parameters which can be adjusted to optimize its denoising performance. To complicate matters further, many of these parameters are 120 signal dependent. Although the conclusions drawn are well supported, specific results in denoising a particular signal may vary depending on a variety of parameters such as Wiener filter length, wavelet filter length, window length, wavelet type, threshold type, hard or soft threshold option, thresholding of approximation coefficients, and the list continues. The variability makes this a fascinating topic which often leads to more questions than answers. No doubt this subject will stimulate interest for years to come. Future studies might involve the application of a median filter to the wavelet coefficients or some type of hybrid Wiener wavelet threshold method in which a minimum filter coefficient value is set based on Statistical analysis. 121 APPENDIX This appendix lists modified MSE results as defined by (3.11). The MSE results are tabulated for the four test signals Doppler, HeaviSine, Bumps, and Blocks which are pictured in Figure 8.1. Variable parameters include: SNR, Wiener filter order, level of decomposition, wavelet threshold type, and wavelet threshold option (s for soft, h for hard). Actual and independent noise sample refers to the noise provided to the Wiener filter. In the actual case the noise used to produce the noisy signal provided, while the independent case refers to an independent sample of white Gaussian noise with a zero mean and variance as required to meet the SNR level desired. 7 . cameaa . am oes cecaikntcaaete a pee x ag oe P Vig RET, ame s aise? or hanes AE om 4d 9 24 LOR Cyl PS = a ae 4 5 igs Tae o s te! a ME Ci BE Se Bie RIS: oes, anita “Glhe. Be Sipe 205 Re Bs EEG oh O a Fie ea. Pi gitar ee Gace Ne Se A Sn TD / oe yee ee BX ne ate ht 4’ Misa "4 . ig NR RE A S53 ML Gal sote bi YP ss | oss | oacs | oamzs | 02019 r re om ee Ti = 0.1543 0.06870 0.1415 0.08855 [0 [2% | oasi0 | o0m20 | o.v0e | o0sese 13 : yan Io = : = SE — =— , ee oj P at he ae eee ne ee FO ry ce Te eke a ean ae ae ee Senter 7 = - mm 3-2 "Se ef Ree — 6 i a ee bes (as, . Sed eee} ease Baad : =. 7 F Br. eS ee re eae ; rte : ee ~¢ 3 e b 54 - Cao es ae. ae Pe va stirs 5 f A”: tw : be ho Pe ae tue ~*~ ee i ae a Sas ee a Se Sheet ah eet ae Fe me at ta s 4 Rare od 3% me May? 3558 Rw? er, jn ite ak SOS : - tes - sagtes’>) > diay eS ~~ & ae ee Ke oe $ e ie; hae aes es ~% Y~ gi) — > ] ] I ] ] I [20 [8 | eons | c00asee | noises | 0005s Tart |b lmaice ones | oan ig py Ea D 8 0 0 5 15 5 3 ere eS Pere re SURE s 0.1005 s 0.04833 s 0.08631 s 0.06316 h 0.2228 h 0.17490 h 0.17420 h 0.13840 HE s 0.1326 s 0.03152 s 0.14560 s 0.07379 h 0.1241 h 0.04417 h 0.12580 h 0.07464 | 0 UNIV s 0.1470 s 0.03136 s 0.1638 s 0.10510 | | h 0.1357 h 0.03705 | h0O.1454 h 0.09166 124 Ss “WAVELET. ‘THRESHOLDING _ s 0.1190 s 0.03637 s 0.1268 s 0.08342 ——" ail 1560 h 0. =a h 0. — s 0. | §0.04490 s 0. | $0.02456 S | 50.03505 03505 S | 50.03082 03082 Se s 0.05948 s 0.01982 s 0.03499 s 0.03760 Es wav en errs wre ee ros MINI | s 0.05178 s 0.02136 s 0.04663 s 0.03469 SURE s 0.01707 s 0.006897 s 0.01335 s 0.01524 h 0.04011 h 0.02166 h 0.02437 h 0.02763 HEUR s 0.02027 s 0.005605 s 0.01629 s 0.02055 h 0.02043 h 0.007024 h 0.01555 h 0.02004 10 UNIV s 0.04507 s 0.005605 s 0.05424 s 0.04428 h 0.01955 h 0.006310 h 0.02529 h 0.02796 s 0.02485 s 0.006006 s 0.02930 s 0.02766 h 0. eg h 0. ene h 0. Fe — 0.02427 S / $0.006774 006774 | s | 0.003568 003568 | s | §0.005015 | 005015 | sO. } s 0.006461 ee ee es s 0.01571 s 0.003350 s 0.01632 s 0.01962 oe Patt ae tae 20 SURE | s0.002439 | s0.001517 | s 0.001941 | s 0.002452 h 0.005082 | h 0.003641 | h 0.024370 | h 0.003784 20 HEUR | s 0.003283 | s 0.002278 | s0.001953 | s 0.002452 h 0.003231 | h0.001900 | h 0.002258 | h 0.003784 UNIV s 0.008027 | s 0.002587 | s 0.009745 | s 0.012060 h 0.002472 | h0.002176 | h 0.002961 | h 0.004294 | ar s 0.003962 | s 0.001910 | s0.004805 | s 0.005942 | _h 0.002943 | h 0.002796 | h0.002950 | h 0.002862 | 125 > * » te « é Ss a ] eta al oS aa) SK ae . - - = “te Fe 8 eee , < a + Wer’ i * a { + “ i. we ‘ fy att = . + Bae! Fags = aes Mtoe, PIG I AOR TI GN Gg PONG RY J ' ie scr ae - nan , emi 3d : ? BR ae wy! MEE Lye seat “Tae Ge st ies” 1 a ES Sgr SS WS a OE gel = iat -— ga all a d Sue! ~*~ ) = + ay > ae a a SE Te th SiS oe TS oe Sees LINE cr aa ek ante ee Or a ee Pes ce Lae sae Fagen EE Set ee eR sat Sees ww PPR, s 0.1932 s 0.09354 s 0.1997 s 0.1221 oles) cease: | coon | woszes ue Naaman. | pale.) caren SURE s 0.09323 s 0.03300 s 0.08491 s 0.05615 h 0.1385 h 0.06595 h 0.10900 h 0.07381 HEUR s 0.11540 s 0.03226 s 0.12690 s 0.07499 h 0.09107 h 0.03374 h 0.08628 h 0.05957 UNIV s 0.11990 s 0.03226 s 0.13910 s 0.08629 h 0.09974 h 0.03250 h 0.09558 h 0.07248 MINI s 0.1100 s 0.03232 s 0.1212 s 0.07864 h 0.1119 h 0.05961 h 0.1003 h 0.06312 5004872 | s001210 | 003748 | 5003000 s 0.05846 s 0.01195 s 0.03773 s 0.03857 s 0.08133 s 0.01194 s 0.09507 s 0.06621 mI Sr SURE s 0.02533 s 0.005288 s 0.01500 s 0.01618 h 0.02893 h 0.009307 h 0.01791 h 0.02224 HEUR s 0.02903 s 0.005371 s 0.01529 s 0.02002 h 0.02278 h 0.005379 h 0.01392 h 0.01575 UNIV s 0.05147 s 0.005374 s 0.05208 s 0.0463 1 h 0.02309 h 0.005384 h 0.01507 h 0.01765 MINI s 0.03664 s 0.005355 s 0.03187 s 0.03208 h 0.02468 h 0.007919 h 0.01619 h 0.01599 =a SURE s 0.01221 s 0.002783 s 0.005936 | s 0.008363 T= ay HEUR s 0.01272 s 0.003171 s 0.006209 | s 0.009316 a5) | UNV s 0.02784 s 0.003214 s 0.02354 s 0.02822 126 2 ov yh oe we oem oe « ste ntti ton Sng at 25 PRET BRE CRAG OES Ne StF a = ne . rom 753 " TE are aan = ins Mp 20m tr ats > a ie Ne dh or yi zs ou Pt bee wipe yp wrt ps : . AE EOE f a. oa —_—. CHA : W WR dial BY. A _ ee sc 7 oS <a 7 = Tait Ygthenry - <4 i em POF eg ae vA ee ee Pr oT | 20 SURE 5 0.005398 | s0.001565 | s 0.002292 | s 0.003784 | h 0.005575 | h0.002023 | h0.002786 | h 0.003806 20 HEUR 5 0.005914 | 0.001863 | s0.002318 | s 0.003784 h 0.004698 | h0.001469 | 0.002182 | h 0.003806 20 UNIV 5 0.01312 | s 0.002333 | s 0.009679 | s 0.014820 h 0.004717 | h0.001546 | h0.002183 | h0.003411 20 MINI 5 0.008285 | s 0.001990 | s 0.005220 | s 0.008439 | h 0.004842 | h 0.001678 | h 0.002328 | h 0.003374 ae >" we . Ne a ae - nen es ih aes es, re C vies ove -. 7 ‘ aie ? SEN a Pb ie Whe fe, DY get ack Sees cS >. | gy more ee: Sf ery OO : 00458 0.04686 0.01820 0.03268 0.02902 I 2d gs 2 4 os aa Oe Ps [os [+ | cows Ts fs [5 | oniico o0essee [onsise [ 002003 ss [emia | omer Laois | 0010s rio fs | «| o0ares | o00ss50 [oor2s7 | consi Tio fs fs | one [0001s [001260 | 001460 ras fs | 3} oovsres | a00ses0 [ooossr0 | aooners ras || «| eowsrse | aco2see [ooossis | ooo ris [ss | o0wsre0 | cone [onmise | coor Tis [16 [soos | 0002730 [onnioas | coors 15 0 To [8 0.004269 0.001505 0.001942 | 0.003290 0.004268 0.001484 0.001943 | 0.003290 128 568 ¥ ios st —_ APSE SS ES EG RS REET 2 ine A fe SE ETS RE BON” AS % ¥ ‘ ag . pf % ~ FIR WIENER WAVELET : a Sai ee ee epeea™ hee eee a as A = “2X nn Sage SO BAe Lan = SS 3 Pe i a AR Cure < Wore t TOMES: rs fs [3 [ ose | ozs | ozs | osece ms fos f+ | oa oats | oat | ozs rs fs fs | oso [oz | o3iso | ose To fs fs | ons o20s | ooo | ones To fs | | cio | oor [eee | ose To fs fs | ons | ei | on | oars Ts fs [| oosesr | 003568 | onmse | o0sis7 Ts fs fs | oom [ons | 0csso | o0zsio Ts fs fs | oosras | oss [oor [ cosas Tie fs [3 | 0226s | ooorsis [oovst | ooisre Foams | 0005303 0.001608 0.002140 } 0.003333 001967 | 0.003295 Wy Tee A $e RS o WAT eS a ae BEES ERS SESS Be to si 5 oA ee SNR LEVEL DOPPLER | HEAVSINE BUMPS BLOCKS a eee Win & om Bese SS ARM CA gee £ ae a gk. 4 at ~ re 4 : 3 a ~ “ O30 hee’ We, es pe pe Pam i SE Se 0 [+ | eonese | ooseoe | oorres0 | o0s0ss sf 5] conta | oosie0 | oorrs | o00s70 0 [3 oeoe | oaoirs | cone | 0003s ef 0.004317 0.001540 0.002096 0.003328 0.004318 0.001499 0.002096 0.003328 “IR WIENER FILTER ‘(Independent Noise Sample) = 5 = era a a Mh ato Tes emma [ones | aos | coi 130 tye Ke = > as 4a Rey Te Pro edge re oer) ve es: yn, 3 ote Res Sp ,260- ys oe See shy: “ty “y @ at ~ ve a R WIENER FILTER (Independent Noise Sample) Se 2, on 5 4 : od Cr ay Some, pt AGE ae Ke a? ye eee np oa oa We Sage aero ee Re “S93. 0.009947 | 0.002998 | 0.005642 | 0.007588 0.004359 | 0.001572 | 0.002141 | 0.003333 on REFERENCES [1] D. Frieden, Principles of Naval Weapons Systems, Naval Institute Press, Annapolis, MD, 1985. [2] R. Urick, Principles of Underwater Sound, McGraw-Hill Book Company, New York, 1975. [3] R. Urick, Ambient Noise in the Sea, Peninsula Publishing, Los Altos, CA, 1986. [4] V. Knuden, R. Alford, and J. Emling, “Underwater Ambient Noise,” Journal of Mar. Research 7, 410, 1948. [5] W. Crouch and P. Burt, “Logarithmic Dependence of Surface-Generated Sea Noise Spectrum Load on Wind Speed,” JASA 51,1066,1972. [6] R.F.W Coates, Underwater Acoustic Systems, John Wiley & Sons, Inc., New York, 1989. [7] R. Mellen, “Thermal-Noise Limit in the Detection of Underwater Acoustic Signals,” JASA 24, 478, 1952. [8] R. Barsant, Jr., Denoising of Ocean Acoustic Signals Using Wavelet-Based Techniques, Master’s Thesis, Naval Postgraduate School, Monterey, CA, December, 1996. [9] K. Frack, Improving Transient Signal Synthesis Through Noise Modeling and Noise Removal, Master’s Thesis, Naval Postgraduate School, Monterey, CA, March, 1994. [10] S. Haykin, Adaptive Filter Theory, Prentice Hall, Englewood Cliffs, NJ, 1996. [11] N. Wiener, Extrapolation, Interpolation, and Smoothing of Stationary Time Series, The M.LT. Press (formerly Technology Press), Cambridge, MA, 1949. (Reprinted as a paperback in 1964). [12] C. Therrien, Discrete Random Signals and Statistical Signal Processing, Prentice Hall PTR, Englewood Cliffs, NJ, 1992. [13] B. Cipra, “Wavelet Applications Come to the Fore,” SIAM News, Vol.26, No.7, November 1993. 3 [14] C. Burrus, R. Gopinath, and H. Guo, Introduction to Wavelets and Wavelet Transforms, Prentice Hall, Englewood Cliffs, NJ, 1997. [15] S. Haykin, Communication Systems, John Wiley & Sons, New York, 1994. [16] D. Gabor, “Theory of Communications,” Journal of the Institute for Electrical Engineers, 93:429-439, 1946. [17] O. Rioul and M. Vetterli, “Wavelets and Signal Processing,” JEEE Signal Processing Magazine, October, 1991. [18] J. Altmann, “Surfing the Wavelets,” http://www.monash.edu.au/cmcm/wavelet/ wavelet.html, February 1998. [19] R. Polikar, “The Wavelet Tutorial: Part IV,” http://www.public.iastate.edu/~rpolikar/ WAVELETS/WTrtutorial.html, November 1997. [20] A. Croisier, D. Esteban, and C. Galand, “Perfect Channel Splitting by Use of Interpolation, Decimation, and Tree Decomposition Techniques,” Int. Conf. on Information Sciences/Systems, Patras, pp. 443-446, Aug. 1976. [21] P. Vaidyanathan and I. Djokovic, “Wavelet Transforms”. In Wai-Kai Chen, editor, The Circuits and Filters Handbook, chapter 6, pages 134-219, CRC Press and IEEE Press, Roca Raton, FL, 1995. [22] G. Strang and T. Nguyen, Wavelets and Filter Banks, Wellesley-Cambridge Press, Wellesley, MA, 1996. [23] P. Vaidyanathan, Multirate Systems and Filter Banks, Prentice Hall PTR, Englewood Cliffs, NJ, 1993. [24] I. Daubechies, “Orthonormal Bases of Compactly Supported Wavelets,” Communications on Pure and Applied Mathematics, Vol. XLI, 1988. [25] The Mathworks, Inc., Wavelet Toolbox, Natick, MA, 1996. [26] D. Donoho, I. Johnstone, “Ideal Spatial Adaptation by Wavelet Shrinkage,” Biometrika, 81:425-455, 1994. Also Stanford Statistics Dept. Report TR-400, July 1992. [27] D. Donoho, I. Johnstone, “Adapting to Unknown Smoothness via Wavelet Shrinkage,” Journal of American Statist. Assn., to appear 1995. Also Stanford 134 Statistics Dept. Report TR-425, June 1993. [28] D. Donoho, “De-noising by soft-Thresholding, ” JEEE Transactions on Information Theory, 41(3):613-627, May 1995. Also Stanford Statistics Dept. Report TR-409, November 1992. [29] H.Gao, “Threshold Selection in WaveShrink,” Technical Report, StatSci Division of MathSoft, Inc., Seattle, WA, 1995. [30] C. Stein, “Estimation of the Mean of a Multivariate Normal Distribution,” The Annals of Statistics, Vol. 9, pp.1135-1151, 1981. [31] M. Ledbetter, G. Lindgren, and H. Rootzen, Extremes and Related Properties of Random Sequences and Processes, Springer-Verlag, New York, 1983. [32] D. Donoho, “Wavelet Shrinkage and W.V.D. - A Ten Minute Tour,” Technical Report TR- 416, Statistics Department, Stanford University, Stanford, CA, January 1993. [33] L. Rabiner and B. Juang, Fundamentals of Speech Recognition, Prentice Hall PTR, Englewood Cliffs, NJ, 1993. [34] H.V. Poor, An Introduction to Signal Detection and Estimation, Springer-Verlag, New York, 1988. [35] D. Donoho, “Unconditional Bases are Optimal Bases for Data Compression and for Statistical Estimation,” Applied and Computational Harmonic Analysis, 1(1):100- 115, December 1993. Also Stanford Statistics Dept. Report TR-410, Nov. 1992. [36] S. Ghael, A. Sayeed, and R. Baranuik, “Improved Wavelet Denoising via Empirical Wiener Filtering,” Proceedings of SPIE, Mathematical Imaging, San Diego, CA, July, 1997. [37] J. Buckheit, S. Chen, D. Donoho, and J. Scargle, ‘““Wavelab.700," http://wwww.wavelab/playfair.stanford.edu, 1996. [38] R. Coifman and M. Wickerhauser, “Entropy-based algorithms for best-basis selection,” JEEE Transactions on Information Theory, vol. 38, pp. 713-718, March 2. ois INITIAL DISTRIBUTION LIST Defense Technical Information Center.......................08. 2 8725 John J. Kingman Rd., STE 0944 Ft. Belvoir, VA 22060-6218 Dudley Kioxer brary 2)<.<. 5: 210s ee Jee es: Zz Naval Postgraduate School 411 Dyer Rd. Monterey, CA 93943-5101 Clairol Od esis © cain, ea in. 4. nacre ene er reek evens ae eee ] Department of Electrical and Computer Engineering Naval Postgraduate School Monterey, CA 93943-5121 BroteNienique arenes Code BG(Ka. maw 2.24 ss feed ee Ys Department of Electrical and Computer Engineering Naval Postgraduate School Monterey, CA 93943-5121 Prot, Ralplptippenstiel, Code BC/H1 . a... .... 65 24 see te eee oo ] Department of Electrical and Computer Engineering Naval Postgraduate School Monterey, CA 93943-5121 INGW Ci3 ualamon 3200 Rim 387 ....%6i0< faweess ia seas He es ] att: Mr. Steve Greneider, Code 2121 Newport Division Newport, RI 02841 NUWC, Building 1320, Rm 387............. 2.00 eee eeeeeeeeee l att: Dr. Paul M. Baggenstoss, Code 2121 Newport Division Newport, RI 02841 OOD rednic IW Onomicye line... teeta es 6 5 dis ws ee one ee 3 12008 Green Court Glenn Dale, MD 20769 La ABNF i Miya, ae: i00 adm beter rr Ors a 04/30 buf nindsd. Seeesregs & pA ‘4 Vm orm —~eh ~~ Perr eek Ae A IRS Ec OR setesapreage Oe eo ieee Send a cetes -4:8 3 : i “0%. @$ 0° Obi: oa tees oa eves ee ll i Se... Lars + Kees ae a i - 7 _ ‘ a : ‘ — P o2 Li e : oo cone Prater s tate Swe ssagerene Soneta® wags * ’ a ee,t ly corny nn PAOD dn By Ae Ryo O~ 02% eat? Dott 0m. feces @ t KNOZA 3/5 ee Ci : . i . 2S a = = _ t , ‘ t : oe aoe’. « eA a perp oar ons bapa At re : : A ath Soft Or dp Snt. 4p & So Oo Oo Fo eo . we | . ee Pettetint ein <tttn te 2 co ; rotwr) by cud eyagte { 1 : “7 5 mo Meo 61 mete 4 as 8 raed ty * , « ; x 8° q ; corte ye) “ 'o Srey 2.4 ae 6 a : ; A : PE care pre CRO aR Rie wes 2 ae Oe. e00’ Sime ; eo a ies - Cat Biatn® Felie® x * waa, oF? eee §=§— 3 2768 00307936 6 i a e.atue + i ae k ete’ ‘ cant oe ry } ° 7 Pre reece. * ox : mee @e tf catat Oqteen OE" > 6.0.4 && ‘ elt ee0ed = erane A PAT bd e « <6 eta 8h ¥ E F ee 4 wngeete a q : ; i be 4 ° Were by ae ee te es Se ris a ah math 228. tots ote e e Reds fon. Letw set. @ OO ef Sac le . asiores eis ee ee fap trie et Ooft « nt ag PM Me oe ta tS Boe } " . at - "he oe 5 te ores &° oe aig a” sO ; ¥ - e « @ r) . ° estet «Oe. ay teste # : e ee? G20 fie : dene : "i atic oh 6 d e e« ° «t-@ : %e ‘ Pe id head ete Gis fot O. ; ei o? me ¥ a e © ee” e 8. ¢ s& 200 Dr 2 y foe 0 0 e UM © ® tai yeas %» a ewe re Le rn : ooe% ent Pad) aed toot ’ f, : ate @ be ph Goh 828. Oe 4 r , ® he etree tise ny 4 : ant . mene 6 te Mle alta epg ‘ are Ore hy Vast Tal Tar ae os] i er hn 8 oe ES ag aT bal : ar Rate al to? ae i Dotpled Petits ae ; wash, ’ . ‘ ‘ 4 . : ot 4 : mw 5 o ‘ a rece F ee . cape’ 3? ie aw oe ® 4 . ° e 4 he ‘ e . is = oe ) eo 8 mat 28 j 8 . . . ‘ e 7 ® an * 8 neta no Mi £ mt ep be8 & & borer . a6 @ @66 y em terme & ry et.er ¢ Pa - oe ts 28 "ee 10 & Ps BIe e, ape 8 4. * e beng te Gerywee% * = @ e.mee,? eu %@ es att eine or dos Rete ua te one rw-p Peery Se poi tem 9 gta A Rate bbe Behe 02g be tof Oe @ jet ant ot eGo? -e BG ote - . 4 Roh we On er, On keg dyes © eat er sf cS Ont ley ® de Fa Ge 00° Ota . ® a he 9 © 3 : ‘ om sh *et J ¢ . ‘ » . : ° ° Taina of to hela oh btee %- ‘ Ay rv 35-1.- Stiatnie on @ te ot riate®- : 4 : 20 6 ote alate © oe rife . . 8 atatett rte 29s eae ¢ ® ee . 5 wie! aight st toate » aber =o! ster meta & te? : ; © we Pee tem Pp. » - oe *o = Been katy . - her ew ? we oo ip? ‘ “es , - en fe . ° . . ° e . ° ' . 1 ° . e ® ® 8 ¢ e ® rire : ei ee ot we ( ‘I # iA ms ° _ j i eteaes werk ¢ te ' ; a j t ~~