NOTE: If you don't have a copy of the Chapter 8 lecture notes, print it, and put in your binder. There are widget programming examples in the notes that will prove useful to you in this exercise and later on. 
The exercise has four main purposes:

Subsections  

CVS  
FT Aliasing and Orthogonality  
FT Wave experimentarium  
FT Power Spectra  
Wavelets experimentarium  
Wavelets Autocorreletion  
Wavelet and FT data comparison  
Home Work 
CVS 
[about 1 minute] 
To extract the exercise files for this weeks exercise, do
cd ~/ComputerPhysics cvs update d
In case of problems, see the CVS update help page.
FT Aliasing and Orthogonality 
[15 minutes] 
products.pro
program to explore if this constraint is
fulfilled. The procedure
shows why you need to integrate over the full length of the domain
to reach the correct result....

So, if you have a case with 8 intervals (as in aliasing.pro) and you choose a wave with 6 wavelengths ...

FT Wave experimentarium 
[about 30 minutes] 
play.pro
example in IDL.
"eta = 0.1"
first.
eta
and play with the parameters.
Do this either ...
damping
)
in an expression such as
exp(0.001*damping*t*omega)
, to make the
damping proportional to the
period for each component (the factor 0.001 is there to make
the damping small
enough  you may want to change it, depending on the range
of the slider).
exp(0.001*damping*t*omega)
,
to make the damping proportional to the phase for each
component (the factor 0.001 is
there to make the damping small enough  you may want to
change it).
.reset
(note
the dot).
Try answering the following question.
It is hard to see directly in the plot, but you don't have to; you
can
deduce the answer from what "overtones" means, counting on your
fingers
(careful now, this is a tricky question):

FT Power spectra and auto correlation 
[about 15 minutes] 
IDL> !x.style=1 ; exact yaxis IDL> !y.style=1 ; exact yaxis IDL> !p.multi=[0,3,2] ; 3 by 2 frames in a window IDL> window,xsize=1024,ysize=650 ; big window IDL> restore,'signal.sav' ; restore saved data IDL> contour,signal ; signal is a fltarr(1024,32) IDL> f=signal(*,10) ; one line of data IDL> plot,f ; plot it IDL> plot,abs(fft(f,1))^2 ; plot its power spectrum IDL> plot,abs(fft(f,1))^2,xrange=[0,200] ; lowest 200 frequencies IDL> plot,fft(abs(fft(f,1))^2,1),xr=[0,150] ; autocorrelation over 150 days


IDL> power=complexarr(1024,32) IDL> auto=complexarr(1024,32) IDL> for i=0,31 do power[*,i]=abs(fft(signal[*,i],1))^2 IDL> for i=0,31 do auto[*,i]=fft(power[*,i],1)
A pattern is obvious in surface and contour plots of power spectra and the autocorrelation function.
IDL> surface, power[0:100,*], xtitle='frequency', ax=70, az=10 IDL> surface, auto[0:100,*], xtitle='time', ax=70, az=10 IDL> contour, auto[0:100,*], xtitle='time' ; not so neat, but easier to read from IDL> contour, auto[20:40,*], xtitle='time' ; zoom in IDL> plot,fft(abs(fft(f,1))^2,1),xr=[0,150] ; autocorrelation over 150 days


Wavelet experimentarium 
[45 minutes] 
Choose the Morlet wavelet with the smallest order. This shows that the power peaks at time 180 with a signal length of order 64 and with a long tail extending into shorter scales for increasing time. This indicates the decrease in wavelength of the sinusoidal wave with increasing time.



Wavelet Autocorrelation 
[20 minutes] 
where f is the function and T is the time interval over which the function is to be analysed. Integrating this with different values of tau, large values are only obtained when tau represents time differences equivalent to repetitions in the signal of f. From this we define the wavelet autocorrelation function:
where W in general is the complex wavelet transform of a function. Here we are only interested in the real part of the wavelet autocorrelation function. In the wavelet transform a represents the different "time" scales over which the wavelet kernel is applied (see fig 13.10.7 in NR). Numerically data are represented by discrete values along the "time"axis, and here the integrals transform into a sum.
To test it out lets start looking at a single sine function:
IDL> x=findgen(1024)/1024. ; define x [0,1[ IDL> fx=sin(x*!pi*2.) ; single sinewave IDL> plot,x,fx ; plot sinewave IDL> w = wtn(fx,4) ; make wavelet, Deubechies, 4th order IDL> cw = w ; array for containing AC function IDL> for i=1,9 do begin for j=0,2^i1 do $ ; calculate AC function IDL> cw(2^(i1)+j) = total(w(2^(i1):2^i1)*shift(w(2^(i1):2^i1),j)) IDL> cwcw=fltarr(512,10) ; array for expanding the information IDL> ; spread for easy plotting IDL> for i=1,9 do cwcw(*,i) = rebin(cw(2^(i1):2^i1),512) IDL> tvscl,rebin(cwcw,512,100) ; rebin in y direction for more pixelsNow redo the calculation with different wavelenghts and different orders of wavelenghts (limited choices available!). Notice how information about the different wavelenghts are located in different horizontal bands in the tv plot image.

Try a few more combinations to get more familiar with the information in the AC plot:
For this you want a different color table to make it easier to see small amplitude variations. Try the various ones to find one that fits the purpose, you can do this using xloadct.
IDL> w = wtn(signal(*,0),4) ; Wavelet transform IDL> cw = w ; array for containing AC function IDL> for i=1,? do begin for j=0,2^i1 do $ ; calculating the AC function IDL> cw(2^(i1)+j) = total(w(2^(i1):2^i1)*shift(w(2^(i1):2^i1),j)) IDL> cwcw = ??? ; array to contain the variable for plotting IDL> for i=1,? do begin cwcw(*,i) = rebin(cw(2^(i1):2^i1),512) ; stretching the information IDL> tvscl,rebin(cwcw,512,100) ; showing the AC in a positionwavelength plot

Wavelet and FT data comparison 
[20 minutes] 
IDL> na=512 ; Number of data points IDL> a=bytarr(na) ; Define an array to contain the data (saved in byte format) IDL> openr,10,'chirp.dat' ; Open the file for reading IDL> readu,10,a ; Read the data (binary format) IDL> close,10 ; close data file
IDL> window,0 ; opens window number 0 IDL> !p.charsize=2 ; doubles the character size IDL> !p.multi=[0,3,1] ; allows for three plot frames in the horizontal direction IDL> plot,a,xtitle='time',title='Function values',ytitle='Amplitude' IDL> .... ; calculate the Power spectrum for a IDL> ... ; calculate the autocorrelation function IDL> x=indgen(na)na/2+1 ; define the xaxis for the shifted plots IDL> plot_io,x,shift( ??, na/21),xtitle='Wave number',title='.... ; plot the power spectrum IDL> plot,x,????
IDL> window,2 ; opens a new window IDL> ff=fft(a,1) ; define the FFT of a IDL> level=.0100 ; define the threshold limit IDL> in=where(fa lt level, num) ; determine the k values in index space that needs to be zeroed IDL> plot_io,x,???? ; plot the power spectrum IDL> oplot,?? ; a horizontal line indicating the cutoff limit IDL> if num gt 0 then ff(in)=0. ; set the chosen amplitudes to zero IDL> print,??? ; the number of nonzero amplitudes IDL> b=?? ; The inverse FFT of the truncated amplitude information IDL> plot,b,xtitle='time',title='Function values' ; plot the truncated data set IDL> oplot,??,col=200 ; add the original data set to the plot IDL> plot,???,xtitle='time',title='Differences' ; plot the different between the two data sets IDL> print,rms(????) ; prints the rms value (real number) of the difference between the two data sets


Home Work 
[about 2 hours] 
Spend one of the reading hours of this week reading the chapter about Fourier Transforms and wavelet Transforms in Numerical Recipes. If you have a possibility to play with IDL at home, or by spending some odd hours at the Institute, then use the opportunity to play with the Xwidget interface; modifying the play.pro example, or using examples from the IDL manuals, or inventing examples of your own.