# display phase and gian for a few first order low pass filters

set terminal  pdf color enhanced size 10cm,7cm font "Helvetica,10pt"
set output "firstlowpass.pdf"

# om is the product (omega Delta t), k is the shift coefficient 
h(k,om)=(2**(-k))/(1 -(1-2**(-k))*exp(-I*om))
# h(k,om)=(2**(-k))/(exp(I*om)-(1-2**(-k)))

decibel(x)=20*log10(abs(x))
phase(x)=180.0/pi*arg(x)

set multiplot

# magnitude plot
set origin 0,0.5
set size 1,0.5

set logscale x
set xrange [0.0001:1]
set yrange [-30:0]
set format x ""
set ylabel "|h({/Symbol w})| in dB"
set ytics 6
set grid x,y
set key left bottom reverse

plot decibel(h(1,2*pi*x))  title "k=1, {/Symbol w}_C=0.1103 {/Symbol w}_S", \
     decibel(h(2,2*pi*x))  title "k=2, {/Symbol w}_C=0.0458 {/Symbol w}_S", \
     decibel(h(4,2*pi*x))  title "k=4, {/Symbol w}_C=0.0103 {/Symbol w}_S", \
     decibel(h(7,2*pi*x))  title "k=7, {/Symbol w}_C=0.0012 {/Symbol w}_S"

# phase plot
set origin 0,0
set size 1,0.55

set xrange [0.0001:1]
set yrange [-90:90]
set ylabel "arg [h({/Symbol w})] / degree"
set ytics 30
set xlabel "{/Symbol w} / {/Symbol w}_S" offset 0,0.7
unset format x
set key off

plot phase(h(1,2*pi*x))  title "k=1", \
     phase(h(2,2*pi*x))  title "k=2", \
     phase(h(4,2*pi*x))  title "k=4", \
     phase(h(7,2*pi*x)) title "k=7"