pro inv_fourier

;to reverse the fourier transformation by using the definition

;that is, calculating the a's and b's (coefficients of cosine and sine) from the c's, coefficients

;of the complex fourier series.

;points in sample , even integer.

n=100 & nh=n/2 ;

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

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

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

signal= dblarr (n) ;measured signal

signal(n/4:3*n/4)=0.95 ;square pulse, centered (signal is symmetric)

;signal=cos(2*2*!pi*x/float(n)) ;just a cosine at freq 2xfundamental.

plot,x,signal,line=5 & wait,3

;fourier transform

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

sum=dblarr(n)

sum=sum +float(fsignal(0))

;sum will contain inverse transform

for k=1,nh-1 do begin ;Nyquist frequency is one half of sampling frequency, or 50/100.

angfreq=2.*!pi*k/n ;angular frequencies of the Fourier series

;print,angfreq

a=float(fsignal(k)) + float(fsignal(n-k))

;in IDL's fft, the kth. complex term occupies position k, the -kth term is

;in position n-k. Thus the -1 term is at fsignal(n-1), last position.

sum=sum+a*cos(angfreq*x)

plot,x,signal,line=5

oplot,x,sum & wait,0.1

b=imaginary(fsignal(n-k)-fsignal(k))

;the sine coefficient, b_k, is the imag. part of the -kth complex coefficient minus the kth. ;-kth

print,'a(',k,') =',a,' b=,',b

sum=sum+b*sin(angfreq*x)

;plot,x,sum & wait,0.7

;if k eq 2 then stop

endfor

a=float(fsignal(nh))

angfreq=!pi

sum=sum+a*cos(angfreq*x)

;sum = sum + float(fsignal(0))

plot,x,sum ;recovered original function by explicit sum of sines and cosines.

rec_signal=fft(fsignal,1,/double) ;recovered original function by inverse fourier transform

oplot,rec_signal,line=2,thick=2 ;compare recovered signal

oplot,signal,line=5 ;with the original

stop

end