filter_c.pro

IDL语言编写的用于天文自适应光学仿真的软件CAOS V6.0的第一部分。

; $Id: filter_c.pro,v 1.2 2002/03/14 11:49:12 riccardi Exp $

pro filter_c, c_num_p, c_den_p, dt, f_min=f_min, f_max=f_max, n_freq=n_freq,$
filter_tf=filter_tf, hold_tf=hold_tf, LOG=log, DISCRETE=discrete, $
HI2LOW=hi2low, C_NUM_P_2=c_num_p_2, C_DEN_P_2=c_den_p_2


if n_elements(f_min) eq 0 then f_min=1.0
if n_elements(f_max) eq 0 then f_max=10000.0
if n_elements(n_freq) eq 0 then n_freq=1000.0

if n_elements(c_num_p_2) eq 3 and n_elements(c_den_p_2) eq 3 then begin
	if keyword_set(hi2low) then begin
		c_num_2=rotate(c_num_p_2,2) 
		c_den_2=rotate(c_den_p_2,2) 
	endif else begin
		c_num_2=c_num_p_2
		c_den_2=c_den_p_2
	endelse
endif

if keyword_set(hi2low) then begin
	c_num=rotate(c_num_p,2) 
	c_den=rotate(c_den_p,2) 
endif else begin
	c_num=c_num_p
	c_den=c_den_p
endelse

if keyword_set(discrete) then begin
	cost=c_den(2)+c_den(1)*dt+c_den(0)*dt^2
	c_o=[2*c_den(2)+c_den(1)*dt,-c_den(2)]/cost
	c_i=[c_num(2)+c_num(1)*dt+c_num(0)*dt^2,-2*c_num(2)+c_num(1)*dt,c_num(2)]/cost
	print,'First biquadratic discrete filter coefficients'
	print,c_o
	print,c_i
endif

freq = mk_vector(n_freq, f_min, f_max, LOG=log)
hold_tf = sin(!pi*freq*dt)^2/(!pi*freq*dt)^2*exp(-2*complex(0,1)*!pi*freq*dt)
filter_tf = poly(2*!pi*complex(0,1)*freq, c_num)/ $
			poly(2*!pi*complex(0,1)*freq, c_den)*113e-6
if n_elements(c_num_p_2) eq 3 and n_elements(c_den_p_2) eq 3 then begin
	filter_tf_2 = poly(2*!pi*complex(0,1)*freq, c_num_2)/ $
			poly(2*!pi*complex(0,1)*freq, c_den_2)
endif else begin
	filter_tf_2=1.0
endelse

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

filter_tf=filter_tf*filter_tf_2

plot_oo, freq, abs(filter_tf), xtitle='Frequency [Hz]', ytitle='bitC/bitP'
plot_oi, freq, imaginary(alog(filter_tf))*180/!pi, xtitle='Frequency [Hz]',$
								ytitle='Phase [deg]'
plot_oo, freq, abs(filter_tf*hold_tf), xtitle='Frequency [Hz]', ytitle='bitC/bitP'
plot_oi, freq, imaginary(alog(filter_tf*hold_tf))*180/!pi, xtitle='Frequency [Hz]',$
								ytitle='Phase [deg]'
!p.multi=pmulti_save
!p.multi(0)=0
return
end