REM     FIR FILTER DESIGN      firdsn3.bas
REM     Bandpass, Kaiser Window
REM     RSL  27 Dec 98

dim g(150), w(150), b(301)

pi = 3.14159265359#
print
print "FIR Filter Design, Low-pass, Band-pass or High-pass"
input "Number of FIR coefficients? ", nf
input "Sample rate, Hz? ", fs
print
print "The following two entries define the pass band of the filter."
print "The first number is the start of the pass-band.  If it is 0 the filter"
print "is low-pass.  The second number (upper cutoff frequency) is the"
print "top of the pass-band.  If it is set to half of the sampling frequency"
print "the filter is high-pass."
print
input "Lower Cutoff Frequency, Hz, between 0 and half of sample rate? ", f1
input "Upper Cutoff Frequency, Hz, between 0 and half of sample rate? ", f2
input "Stop-band Attenuation, dB (eg 55.0)? ", db
print

f1 = f1 / fs
f2 = f2 / fs
n = int((nf + 1) / 2)

' The variable odd=1 when nf is an odd number, otherwise odd=0.
odd = 0
if 2 * INT(n + .1) - 1 = nf then odd = 1

' Calculate the FIR coefficients g(j) based on the Fourier transform of
' a rectangular frequency response.
c1 = f2 - f1
if odd = 1 then
    g(1) = 2! * c1
    j0 = 2
else
    j0 = 1
end if

for j = j0 to n
    xn = j - 1
    if odd = 0 then xn = xn + .5
    c = pi * xn
    c3 = c * c1
    g(j) = SIN(c3) / c
    g(j) = g(j) * 2! * COS(c * (f1 + f2))
next j

' Find beta, a parameter that determines the out-of-band attenuation.
if db > 50 then
    beta = .1102 * (db - 8.7)
elseif (db > 20.96 and db <= 50!) then
    beta = .58417 * (db - 20.96) ^ .4 + .07886 * (db - 20.96)
else
    beta = 0!
end if

' Find Kaiser window factors w(j) that limit the out-of-band response.
x = beta
gosub Bessi0
bes = in0
xind = (nf - 1) * (nf - 1)

for j = 1 to n
    xi = j - 1
    if odd = 0 then xi = xi + .5
    xi = 4! * xi * xi
    x = beta * SQR(1! - xi / xind)
    gosub Bessi0
    w(j) = in0 / bes
next j

' Window the ideal filter response:
for j = 1 to n
    g(j) = g(j) * w(j)
next j

' Rearrange to be causal filter coefficients:
for j = 1 to n
    b(j) = g(n + 1 - j)
next j

for j = n + 1 to nf
    if odd = 1 then
        b(j) = g(j - n + 1)
    else
        b(j) = g(j - n)
    end if
next j

for j = 1 to nf
    print "Coefficient "; j; " = "; b(j)
next j

input "Hit Enter to continue", Z$
stop

'  Bessi0  -  Returns in0, Modified Bessel function Io(x) for any real x.
'  From Press, et. al., Numerical Recepies in C, 2nd Ed. p237.
Bessi0:
ax = ABS(x)
if ax < 3.75 then
    y = x / 3.75
    y = y * y
    exp0 = .2659732# + y * (.0360768 + y * .0045813)
    in0 = 1! + y * (3.5156229# + y * (3.0899424# + y * (1.2067492# + y * exp0)))
else
    y = 3.75 / ax
    exp0 = y * (-.01647633# + y * 3.92377E-03)
    exp1 = y * (-.02057706# + y * (.02635537# + exp0))
    exp2 = y * (2.25319E-03 + y * (-1.57565E-03 + y * (9.16281E-03 + exp1)))
    in0 = (EXP(ax) / SQR(ax)) * (.39894228# + y * (.01328592# + exp2))
end if

return

end

' ## End of file ##