pro WDOSHRNK
;+
; NAME:
;	WDOSHRNK
;
; PURPOSE:
;	This procedure denoises a time series, using wavelet shrinkage,
;   and plots the result. The time series is the working dataset, 
;   the wavelet is the current one selected, and the shrinkage method
;   is the one chosen under the Wavelet Workbench Settings: DeNoising
;   Threshold.
;
; CATEGORY:
;	Wavelets.
;
; CALLING SEQUENCE:
;	WDOSHRNK
;
; NECESSARY INPUTS (via COMMON)
;	sig:	1-d signal array	
;	QMF: Quadrature Mirror Filter
;	LD: Coarseness level (set by WINTWAVE)
;	st: signal title
;	noise_lev: noise level (set by NoLev procedure, =0 by default)
;	noise_type: noise type (set by NoType procedure, 
;		1, 2 or 3 = Normal, Uniform, Poisson,  =1 by default)
;	shrinkage_type: shrinkage method (set by Threshold procedure, =1 by default)
;		1 = VisuThresh
;		2 = MultiHybrid
;		3 = MultiMAD
;		4 = MinMaxThresh
;		5 = WPDeNoise
;	WaveType: Wavelet Type: 'Haar', 'Daubechies', 'Coiflet', 'Symmlet'
;	ParVal: Parameter Value of the above wavelet, specific for the
;		wavelet. For:
;		WaveType = 'Haar'
;		ParVal = 1
;		WaveType = 'Daubechies'
;		ParVal = 4 or 6
;		WaveType = 'Coiflet'
;		ParVal = 2 or 3
;		WaveType = 'Symmlet'
;		ParVal = 6 or 8
;
; OUTPUTS:
;	The denoised time series is plotted to the screen.
;
; SEE ALSO:
;	wintwave (necessary for this procedure)
;
; EXAMPLE:
; 		IDL> ;Truncate sig array appropriately and get QMF
;		IDL> WINTWAVE, sig, WaveType, ParVal, QMF, LD
;		IDL> ;Need to normalize first, if working with natural data)
;		IDL> sig = WNNOISE(sig,QMF, tt)
;		IDL> If tt EQ `S' then WPLOTIT, st, 1, sig
;		IDL> WDOSHRNK
;
; MODIFICATION HISTORY:
; 	Written by:	Amara Graps		September-November 1995
;Translated from MatLab Wavelab routine: do_shrinkage.m

;-
COMMON WWB_COMMON, $
	TEXT_ANNOUNCE, wd, sig, len, st, noise_type, noise_lev, $
	QMF, LD, shrinkage_type, WaveType, ParVal
COMMON Wave_color, $
	max_color, max_image, $
	white, yellow, lavender, aqua, pink, green, red, orange, blue, $
	lt_gray, med_green, brown, olive, purple, dk_gray, black


x_work = sig
TT = WSIGTYPE(x_work, olen, error)

;Color of text strings
IF !p.background EQ 0 THEN cl = 255 ELSE cl = 0


CASE error OF
	1: BEGIN
		;Valid array or matrix
	
		CASE TT OF
			'I': BEGIN
				;Image
				WIDGET_CONTROL, TEXT_ANNOUNCE, SET_VALUE='Shrinkage Not Implemented for 2D Data.'				
				END	;image
				
			'S': BEGIN
				;Signal
				;Load basic color table and set the Wave_color common variables
				WVLOADCT, 2
				
				
				;++++++++++++++++++++++++++++++
				;Do Forward Wavelet Transform +
				;++++++++++++++++++++++++++++++
				
				stitle = st + ' (Before)'
				plotwin = 3
				WDODWT, x_work, QMF, LD, stitle, plotwin, wc
				
				;(NOTE: If we add more wavelets, make sure that we have a case 
				;statement here to perform the proper forward transform
				;operation, i.e.  FWT_CDJV for Daubechies, edge corrected etc.)
				
				nx = N_ELEMENTS( x_work ) 
				
				;+++++++++++++++++++++++++
				;Do Shrinkage Operation  +
				;+++++++++++++++++++++++++
				
				CASE shrinkage_type OF
				
					1: BEGIN
						;VisuThresh
						TThresh = 1	;Hard Thresholding (Soft Thresholding = 2)
						wc(2^LD:nx-1) = WVISUTHR(wc(2^LD:nx-1), TThresh)  
						END	;1
						
					2: BEGIN
						;MultiHybrid
						wc = WMHYBRID(wc,LD) 
						END	;2
					
					3: BEGIN
						;MultiMAD
						wc = WMULTMAD(wc,LD) 
						END	;3
						
					4: BEGIN
						;MinMaxThresh
						wc(2^LD:nx-1) = WMTHRESH(wc(2^LD:nx-1))  
						END	;4
						
					5: BEGIN
						;WPDeNoise
					  	D = 10 
						clean = WPDNOISE(x_work,D,QMF,btree,stree) 
						END	;5
						
						
					ELSE: print, 'shrinkage_type unknown!'
				
				ENDCASE 
				
				
				;++++++++++++++++++++++++++++++++++++++++++++++++++++++
				;Plot dyadic wave coeffs if not wavelet packet method
				;otherwise display best basis tree. Then, perform the 
				;inverse wavelet transform and create a denoised signal.
				;+++++++++++++++++++++++++++++++++++++++++++++++++++++++
				;Set up the plot window
				;dummy values
				plotwin = 4
				pt = '' & dd = fltarr(2)
				WPLOTIT, pt, plotwin, dd
				
				CASE 1 OF
				
					( shrinkage_type LT 6): BEGIN
					
						;Plot the dyadic coeffs
						scale_2 = 0 * MAX(ABS( wc ) )
						!C = 0	;reset from MAX function
						
						signaltitle = 'WT of ' + st+' (After)'
						stitle = 'WT (Dyadic coeffs) of ' + st+' (After)'
								
						;Stack two plots vertically on the page
						!p.multi = [0, 1, 2]
				
						;Plot the coeffs in a spike plot for each of the dyads
						WPLWCOEF,wc,LD,0., stitle
				
						;Plot the coefficients in bottom plot
						PLOT, wc, title = signaltitle, ystyle = 1, color = cl
				
						!p.multi = 0	;reset to no multiple plots
				
						;Do the inverse transform for the 4 Wavelet types:
						;Haar, Daubechies, Coiflet, Symmlet 
						;(NOTE: If we add more wavelets, make sure that
						;we have a case statement here to perform the proper inverse 
						;operation, i.e.  IWT_CDJV for Daubechies, edge corrected etc.)
						
						xhat = WIWTPO( wc, LD, QMF )
						
						END
					
					ELSE: BEGIN
	
						;Stack two plots vertically on the page
						!p.multi = [0, 1, 2]
							
						WPLBTREE, btree,10,stree,st
				
						WPHPLANE,'WP',btree,wp,titlestr=st, nTFR = 64

						!p.multi = 0	;reset to no multiple plots
										
						;We already have a cleaned dataset from the wavelet
						;packet synthesis
						xhat = clean
						END	;Else
				
				ENDCASE
				   
				
				;+++++++++++++++++++++++++++++++++++++++++++ 
				; Plot the Denoised (and Error) Data Arrays  + 
				;------------------------------------------- 
				;I've commented out the error array.. I'm not sure
				;people are very interested in seeing it.. (AG)
				;stitle='Denoised: '+ st+' (Blue),  Error (Red)'
				
				stitle='Denoised: '+ st
				ymin = MIN(xhat)
				ymax = MAX(x_work)
				!c = 0	;reset system cursor that min/max automatically sets
				
				;Set up the plot window
				;dummy values
				pt = '' & dd = FLTARR(2)
				plotwin = 5
				WPLOTIT, pt, plotwin, dd
				
				;Plot, handling colors correctly. First, draw axes
				;and title using default colors, then draw data points
				;using the common variable colors
				PLOT, xhat ,/nodata, title=stitle, yrange=[ymin,ymax], color = cl
				;plot denoised data
				OPLOT, xhat, color = blue
				
				;plot error
				;OPLOT, xhat - x_work, linestyle=1, color=red 
				
				 
				;+++++++++++++++++++++++++++++++++++++++++ 
				;  Compute and Plot the Power Spectra of  + 
				;    the data and de-noised version       + 
				;------------------------------------------ 
				 
				stitle = 'Power Spectra: Data (green) DeNoised (blue) of '+st
				xtitle = 'log(Frequency)'
				ytitle = 'log(Power)'
				
				;NOTE to myself (AG): ask Jeff why only half of the data is
				;spectrum analyzed.
				
				powersp = WSPECTRM( x_work, nx/2 )
				denoised = WSPECTRM( xhat, nx/2 )


				ymin = MIN([powersp, denoised])
				ymax = MAX([powersp, denoised])
					
				xmin = 1
				xmax = nx
				!c = 0	;reset system cursor that min/max automatically sets
				
				;Set up the plot window
				;dummy values
				pt = '' & dd = fltarr(2)
				plotwin = 6
				WPLOTIT, pt, plotwin, dd
				
				;Plot, handling colors correctly. First, draw axes
				;and title using default colors, then draw data points
				;using the common variable colors
				PLOT, xhat ,/nodata, /xlog, /ylog, title=stitle, yrange=[ymin,ymax],$
				xrange=[xmin,xmax], xtitle = xtitle, ytitle = ytitle, xstyle=1, $
				ystyle = 1, color=cl
				;plot power spectra data
				OPLOT, powersp, color = green
				;plot power spectra denoised data
				OPLOT, denoised,  color=blue, linestyle = 1 


				
				END	;signal/data array						
		ENDCASE	;TT

	
		END	;valid array or matrix
		
	ELSE: PRINT, 'Not valid data/signal array or matrix!'

ENDCASE		;error

;
; Algorithm Source:  WaveLab Version 0.600
; WaveLab WWW site: http://playfair.stanford.edu/ 
; WaveLab Questions? e-mail wavelab@playfair.stanford.edu
;
END  
;***************************************************