; http://www.wavemetrics.com/products/igorpro/dataanalysis/signalprocessing/powerspectra.htm x = (/ 1002, 1017, 1018, 1020, 1018, 1027, \ 1028, 1030, 1012, 1012, 982, 1012, \ 1001, 996, 995, 1011, 1027, 1025, \ 1030, 1016, 996, 1006, 1002, 982 /)*1.0 ; even ;;1030, 1016, 996, 1006, 1002, 982, 999 /)*1.0 ; odd N = dimsizes(x) df = 1.0/N xAvg = avg(x) xVar = variance(x)*(N-1.)/N ; biased estimate x = x - xAvg print ("N="+N+" df="+df+" avg(x)="+xAvg+" xVar="+xVar) print("---") ;****************************************** ; cfftf ;****************************************** cf = cfftf (x, 0.0, 0) ; imaginary part set to 0.0 printVarSummary(cf) print("---") cf = cf*(1.0/N) ; normalization for 2-sided fft print("cfftf: "+sprintf("%9.5f",cf@frq)+" "+sprintf("%9.3f", cf(0,:))+" "+sprintf("%9.3f", cf(1,:)) ) print("---") cf2 = cf(0,:)^2 + cf(1,:)^2 ; sum after normalization cf2sum= sum(cf2) print("---") print("cf2sum="+cf2sum+" ratio="+(cf2sum/xVar)) print("---") ;****************************************** ; ezfftf ;****************************************** cfez = ezfftf (x) printVarSummary(cfez) print("---") print("ezfftf: "+sprintf("%9.3f", cfez(0,:))+" "+sprintf("%9.3f", cfez(1,:)) ) print("---") cfez2 = cfez(0,:)^2 + cfez(1,:)^2 cfez2sum= sum(cfez2) print("cfez2sum="+cfez2sum) print("---") if (N%2 .eq. 0) then N2 = N/2 cfez_var = 0.5*sum(cfez2(0:N2-2)) + cfez2(N2-1) else cfez_var = 0.5*sum(cfez2) end if print("cfez_var="+cfez_var+" ratio="+(cfez_var/xVar)) print("---")