{{{id=159| %auto reset() forget() # # Copyright 2001, P. Lutus http://arachnoid.com # # Released under the GPL: http://www.gnu.org/licenses/gpl.html # # special equation rendering def render(x,name = "temp.png",size = "large"): if(type(x) != type("")): x = latex(x) latex.eval("\\" + size + " $" + x + "$",{},"",name) # NOTE: sample rate must be > 2x highest frequency of interest # just as in reality (Nyquist-Shannon) var('alpha beta cf g sr q df f ff dbGain a0 a1 a2 b0 b1 b2') # default sample rate def_sr = 40000 # default q required for some filters def_q = 1/sqrt(2) # constants that apply to all filters A = 10 ^ (dbGain / 20) beta = sqrt(A + A) omega = 2 * pi * df / sr sn = sin(omega) cs = cos(omega) alpha = sn / (2*q) phi = (sin(2 * pi * f/(2.0 * sr)))^2 # value function rf(sr,q,df,f,a1,a2,b0,b1,b2) = ((b0+b1+b2)^2 - 4*(b0*b1 + 4*b0*b2 + b1*b2)*phi + 16*b0*b2*phi*phi) / ((1+a1+a2)^2 - 4*(a1 + 4*a2 + a1*a2)*phi + 16*a2*phi*phi) # value to db conversion function db(x) = 10 * log(x) / log(10) # from db to value val(x) = 10 ^ (x/10) # vertical line drawing function def vl(x): return line(((x,-1),(x,1))) # lowpass a0 = 1+alpha b0 = ((1-cs) / 2) / a0 b1 = (1-cs) / a0 b2 = b0 a1 = (-2 * cs) / a0 a2 = (1-alpha) / a0 lpf(sr,q,df,f) = rf(sr,q,df,f,a1,a2,b0,b1,b2) # highpass a0 = 1 + alpha b0 = ((1 + cs) /2) / a0 b1 = -(1 + cs) / a0 b2 = b0 a1 = (-2 * cs) / a0 a2 = (1 - alpha) / a0 hpf(sr,q,df,f) = rf(sr,q,df,f,a1,a2,b0,b1,b2) # bandpass a0 = 1+alpha b0 = alpha / a0 b1 = 0 b2 = -b0 a1 = (-2 * cs) / a0 a2 = (1-alpha) / a0 bpf(sr,q,df,f) = rf(sr,q,df,f,a1,a2,b0,b1,b2) # notch a0 = 1 + alpha b0 = 1 / a0 b1 = (-2 * cs) / a0 b2 = b0 a1 = (-2 * cs) / a0 a2 = (1 - alpha) / a0 nf(sr,q,df,f) = rf(sr,q,df,f,a1,a2,b0,b1,b2) # peak a0 = 1 + (alpha /A) b0 = (1 + (alpha * A)) / a0 b1 = (-2 * cs) / a0 b2 = (1 - (alpha * A)) / a0 a1 = b1 a2 = (1 - (alpha /A)) / a0 pf(sr,q,df,f,dbGain) = rf(sr,q,df,f,a1,a2,b0,b1,b2) # high shelf a0 = (A + 1) - (A - 1) * cs + beta * sn b0 = (A * ((A + 1) + (A - 1) * cs + beta * sn)) / a0 b1 = (-2 * A * ((A - 1) + (A + 1) * cs)) / a0 b2 = (A * ((A + 1) + (A - 1) * cs - beta * sn)) / a0 a1 = (2 * ((A - 1) - (A + 1) * cs)) / a0 a2 = ((A + 1) - (A - 1) * cs - beta * sn) / a0 hsf(sr,q,df,f,dbGain) = rf(sr,q,df,f,a1,a2,b0,b1,b2) # low shelf a0 = (A + 1) + (A - 1) * cs + beta * sn b0 = (A * ((A + 1) - (A - 1) * cs + beta * sn)) / a0 b1 = (2 * A * ((A - 1) - (A + 1) * cs)) / a0 b2 = (A * ((A + 1) - (A - 1) * cs - beta * sn)) / a0 a1 = (-2 * ((A - 1) + (A + 1) * cs)) / a0 a2 = ((A + 1) + (A - 1) * cs - beta * sn) / a0 lsf(sr,q,df,f,dbGain) = rf(sr,q,df,f,a1,a2,b0,b1,b2) /// }}} {{{id=366| # bandpass example plot(bpf(def_sr,3,100,f),f,0,200, gridlines=True) ///

}}} {{{id=357| # lowpass example plot(lpf(def_sr,def_q,100,f),f,0,200, gridlines=True) ///

}}} {{{id=358| # highpass example plot(hpf(def_sr,def_q,100,f),f,0,200, gridlines=True) ///

}}} {{{id=359| # notch example plot(nf(def_sr,9,100,f),f,0,200, gridlines=True) ///

}}} {{{id=360| # peak example plot(pf(def_sr,9,100,f,9),f,0,200, gridlines=True) ///

}}} {{{id=361| # highshelf example plot(hsf(def_sr,def_q,100,f,1),f,0,200, gridlines=True) ///

}}} {{{id=362| # lowshelf example plot(lsf(def_sr,def_q,100,f,1),f,0,200, gridlines=True) ///

}}} {{{id=347| # example of modeling multi-filter FM detector cf = 1000 dev = 85 start = 500 end = 1500 @interact def curve( cf = slider([500..2000],default=1000), gain = slider([10..1000],default=200), bpq = slider([1..20],default=6), msq = slider([1..20],default=12), qdev = slider([50..300],default=90), qff = slider([-100..100],default=40), bw = slider([200..1200])): bpff = cf + (qff/10) * 1000 / cf bp(f) = bpf(def_sr,bpq,bpff,f) m(f) = bpf(def_sr,msq,bpff+qdev,f) s(f) = bpf(def_sr,msq,bpff-qdev,f) v(f) = (m(f) - s(f)) * bp(f) * gain * .01 mp = v(cf+qdev).n() zp = v(cf).n() sp = v(cf-qdev).n() print sp,zp,mp p = plot(v(f) ,f,cf-bw,cf+bw) ml = vl(cf+dev) sl = vl(cf-dev) show(p+ml+sl) ///

gain

bpq

msq

qdev

qff

}}} {{{id=365| /// }}}