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



plot,signal & wait,2

;the smoothing kernel in time



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=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([1,5])=0.1 ;convolutions use a short kernel, different from sm_5.





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