pro smooth_5_by_convolution

;smooth_5 implements 5-point smoothing by multiplying in the frequency domain

;smoth_5_by_convolution also uses the direct convolution method.

;the program also demonstrates some quirks in the convolution and fourier calculations made by IDL

;They are not exactly compatible, in the sense that the convolution is defined in the

;engineering style, of a much shorter array, while for fourier you need same size arrays.

;i=complex(0,1.) ;the imaginary unit, not used in this program

n=100 & nm=n-1 ;n = points in sample, even integer.

& nh=n/2 ;

x=findgen(n) ;real variable space. Period is n times the interval of sampling.

;fundamental frequency is 1/n. Fundamental angular frequency is 2 pi / n.

;sampling frequency is 1. sampling angular frequency is 2 pi.

signal= fltarr (n) ;measured signal

for j=0,nm do begin

signal (j)= 10.*abs(nh-j)/nh;-5 ;a "V", changes at 0.2/point

endfor

signal(3*n/4)=0

plot,signal & wait,2

;the smoothing kernel in time

sm_5=fltarr(n)

sm_5(0:2)=1./5.

sm_5(nm-1:nm)=1./5. ;note that this has the same array size as the data

;suitable for fourier transforming, but not for the CONVOL operator below

;sm_5(nh)=1. ;translation operator, by half period

;use it to see how it works!

;note the smoothing array, needs to be centered in position 0!

;*************************

;fourier transforms

fsignal=fft(signal,/double) ;a n term, double precision, complex array.

fsm_5=fft(sm_5,/double)

fsm_5=fsm_5/abs(fsm_5(0)) ;the transfer function of the smoothing filter, note division

fsmsignal=fsm_5*fsignal ;multiplying the signal by the transfer function

smsignal=fft(fsmsignal,1,/double) ;the smoothed signal in time domain (inverse transform)

oplot,x,smsignal,line=2 &wait,2

;*************************

;now will do the same directly in the time domain

kernel=dblarr(7)

kernel([1,5])=0.1 ;convolutions use a short kernel, different from sm_5.

kernel([2,4])=0.2

kernel(3)=0.4

;kernel=[0,0.1,0.2,0.4,0.2,0.1,0]

smsig_conv=convol(signal,kernel,/edge_wrap)

oplot,x,smsig_conv, line = 3, thick=2

stop

end