pro filter

device,decomposed=0
loadct,39

;b_1_1=-1.9587428340882587
;b_1_2=0.96134553442399129
;b_2_1=-1.9066292518523
;b_2_2=0.90916270571237567

b_1_1=double(-2.*32092./2.^15.)
b_1_2=double(2.*15750./2.^15.)
b_2_1=double(-2.*31238./2.^15.)
b_2_2=double(2.*14895./2.^15.)

f_samp=1./(0.00000002*100.*33.)

nf=30000.
freq=(1./(nf-1.))* findgen(nf)* f_samp

!p.multi=[0,1,2]

window,5,xs=850,ys=600,title='filter'

;H(Z) = [1 + 2*Z^(-1) + Z^(-2)] / [1+b1*Z^(-1) + b2*Z^(-2)]

omega=freq/max(freq)*2.*!pi

h1_omega=complex(n_elements(omega))
h2_omega=complex(n_elements(omega))
h_omega=complex(n_elements(omega))

h1_omega=(1.+2.*complex(cos(-omega),sin(-omega))+complex(cos(-2.*omega),sin(-2.*omega)))/(1.+b_1_1*complex(cos(-omega),sin(-omega))+b_1_2*complex(cos(-2.*omega),sin(-2.*omega)))
h2_omega=(1.+2.*complex(cos(-omega),sin(-omega))+complex(cos(-2.*omega),sin(-2.*omega)))/(1.+b_2_1*complex(cos(-omega),sin(-omega))+b_2_2*complex(cos(-2.*omega),sin(-2.*omega)))

h_omega=h1_omega*h2_omega/2048.;/max(h1_omega*h2_omega/2048.)

ind_f3db=where(abs(abs(h_omega(0:300))-max(abs(h_omega(0:300)))/sqrt(2)) eq min(abs(abs(h_omega(0:300))-max(abs(h_omega(0:300)))/sqrt(2.))))

phase=atan(imaginary(h_omega),float(h_omega))/2./!pi*360.

for i=0,n_elements(phase)-1 do if imaginary(h_omega(i)) gt 0 then phase(i)=phase(i)-360.

plot,freq,abs(h_omega),xr=[0.,500.];,/ylog
oplot,[0,n_elements(h_omega)-1],[max(abs(h_omega(0:500)))/sqrt(2.),max(abs(h_omega(0:500)))/sqrt(2.)]
oplot,[freq(ind_f3db),freq(ind_f3db)],[min(abs(h_omega(0:500))),max(abs(h_omega(0:500)))]

plot,freq,phase,xr=[0.,500.]

print,'DC amplification= '+string(max(abs(h_omega)))
print,'F3dB= '+string(freq(ind_f3db))+' Hz'
;print,'Omega= '+string(omega(ind_f3db))
print,'Gain 200Hz= '+string(abs(h_omega(396))/max(abs(h_omega)))

;set_plot,'ps'
;device, /encapsulate
;device, filename= 'C:\users\elia\Vancouver\MCE\sims\filter.eps', /portrait,ysize=18,xsize=26

;plot,freq,abs(h_omega),xr=[0.,500.],xtitle='frequency',ytitle='LPF magnitude';,/ylog
;oplot,[0,n_elements(h_omega)-1],[max(abs(h_omega(0:500)))/sqrt(2.),max(abs(h_omega(0:500)))/sqrt(2.)]
;oplot,[freq(ind_f3db),freq(ind_f3db)],[0,max(abs(h_omega(0:500)))]
;plot,freq,phase,xr=[0.,500.],xtitle='frequency',ytitle='LPF phase';,/ylog

;device,/close
;set_plot,'win'

;readcol,'C:\users\elia\Vancouver\MCE\sims\calc_impulse_p1000_4pole.dat',firm_sims
;output=dblarr(fix(f_samp))
;output(0:n_elements(firm_sims)-1)=firm_sims(0:n_elements(firm_sims)-1)
;fft_output=abs(FFT(output))
;ind_f3db_sims=where(abs(fft_output(0:300)-max(fft_output(0:300))/2.) eq min(abs(fft_output(0:300)-max(fft_output(0:300))/2.)))
;window,3,xs=850,ys=600,title='sims'
;plot,fft_output,xr=[0,500]
;oplot,[0,n_elements(fft_output)-1],[max(fft_output(0:300))/2.,max(fft_output(0:300))/2.]
;oplot,[ind_f3db_sims,ind_f3db_sims],[min(fft_output(0:300)),max(fft_output(0:300))]
;print,'F3dB_sims= '+string(ind_f3db_sims)+' Hz'

stop

end