#band-pass for pulse monitor
fs = 30
A0= 64
A1= -108
A2= 48
B0 = A0
B1 = 0
B2 = -B0

A0_2 = 64
A1_2 = -97
A2_2 = 38

#band-pass for pulse monitor
fs = 240
A0= 256
A1= -489
A2= 234
B0 = A0
B1 = 0
B2 = -B0

A0_2 = 256
A1_2 = -479
A2_2 = 224

#band-pass for pulse monitor
fs = 60
A0= 256
A1= -388
A2= 141
B0 = A0
B1 = 0
B2 = -B0

A0_2 = 256
A1_2 = -449
A2_2 = 199

unset arrow
set arrow nohead from 40000,-20 to 40000,20

unset label
pole1a = ( -A1 + sqrt(A1*A1 - 4 * A0 * A2)) / (2*A0)
pole1b = ( -A1 - sqrt(A1*A1 - 4 * A0 * A2)) / (2*A0)
set label sprintf("pole magnitude %.3f", abs(pole1a)) at 1,-3
set label sprintf("pole1a at %.3f + %.3f j", real(pole1a), imag(pole1a)) at 1,-5
set label sprintf("pole1b at %.3f + %.3f j", real(pole1b), imag(pole1b)) at 1,-7


pole2a = ( -A1_2 + sqrt(A1_2*A1_2 - 4 * A0_2 * A2_2)) / (2*A0_2)
pole2b = ( -A1_2 - sqrt(A1_2*A1_2 - 4 * A0_2 * A2_2)) / (2*A0_2)
set label sprintf("pole2b at %.3f + %.3f j", real(pole2b), imag(pole2b)) at 1,-17
set label sprintf("pole2a at %.3f + %.3f j", real(pole2a), imag(pole2a)) at 1,-15
set label sprintf("pole magnitude %.3f", abs(pole2a)) at 1,-13

set title sprintf("Design of biquad filter, fs=%3g Hz",fs)

set key bottom right
set ylabel "gain [dB]"
unset logscale y
set yrange [-30:*]

set xlabel "frequency [Hz]"
set logscale x
set xrange [0.01:0.5*fs]
# set xrange [1000:0.5*fs]

#set xtics add (0.43, 2.15)
set grid xtics
set mxtics 10

j=sqrt(-1)
biquad(zinv,b0,b1,b2,a0,a1,a2) = (b0+zinv*(b1+zinv*b2))/(a0+zinv*(a1+zinv*a2))
gain(f,b0,b1,b2,a0,a1,a2) = abs( biquad(exp(-j*2*pi*f/fs),b0,b1,b2,a0,a1,a2))
phase(f,b0,b1,b2,a0,a1,a2) = imag(log( biquad(exp(-j*2*pi*f/fs),b0,b1,b2,a0,a1,a2)))

set samples 5000

plot \
    20*log10(gain(x,B0,B1,B2, A0,A1,A2)) \
	title sprintf("%.0f (1 + %.0f z^-1 + %.0f z^-2)/(%.2f+ %.3f z^-1 + %.3f z^-2)", \
			B0, B1/B0, B2/B0, A0, A1, A2), \
    20*log10(gain(x,B0,B1,B2, A0_2,A1_2,A2_2)) \
	title sprintf("%.0f (1 + %.0f z^-1 + %.0f z^-2)/(%.2f+ %.3f z^-1 + %.3f z^-2)", \
			B0, B1/B0, B2/B0, A0_2, A1_2, A2_2) lt 3