fir.py

   1 # -*- coding: utf-8 -*-
   2 """
   3 Python code for basis FIR filter design
   4 @author: Matti Pastell <matti.pastell@helsinki.fi>
   5 http://mpastell.com
   6 """
   7
   8 import sys
   9 from pylab import *
  10 import scipy.signal as signal
  11
  12 #Plot frequency and phase response
  13 def mfreqz(b,a=1):
  14     w,h = signal.freqz(b,a)
  15     h_dB = 20 * log10 (abs(h))
  16     subplot(211)
  17     plot(w/max(w),h_dB)
  18     ylim(-150, 5)
  19     ylabel('Magnitude (db)')
  20     xlabel(r'Normalized Frequency (x$\pi$rad/sample)')
  21     title(r'Frequency response')
  22     subplot(212)
  23     h_Phase = unwrap(arctan2(imag(h),real(h)))
  24     plot(w/max(w),h_Phase)
  25     ylabel('Phase (radians)')
  26     xlabel(r'Normalized Frequency (x$\pi$rad/sample)')
  27     title(r'Phase response')
  28     subplots_adjust(hspace=0.5)
  29
  30 #Plot step and impulse response
  31 def impz(b,a=1):
  32     l = len(b)
  33     impulse = repeat(0.,l); impulse[0] =1.
  34     x = arange(0,l)
  35     response = signal.lfilter(b,a,impulse)
  36     subplot(211)
  37     stem(x, response)
  38     ylabel('Amplitude')
  39     xlabel(r'n (samples)')
  40     title(r'Impulse response')
  41     subplot(212)
  42     step = cumsum(response)
  43     stem(x, step)
  44     ylabel('Amplitude')
  45     xlabel(r'n (samples)')
  46     title(r'Step response')
  47     subplots_adjust(hspace=0.5)
  48
  49
  50 Fe = 5000.
  51 Fniq = Fe/2
  52 transition = 1.001
  53
  54 n = 32
  55 Ffond = 500
  56 Fharm1 = 1500.
  57 Fother = 850.
  58
  59 def bandpass(n, freq):
  60     global Fe
  61     a = signal.firwin(n, cutoff = (freq/Fniq)/transition, window = 'boxcar')
  62     b = - signal.firwin(n, cutoff = (freq/Fniq)*transition, window = 'boxcar')
  63     b[n/2] = b[n/2] + 1
  64     d = - (a+b); d[n/2] = d[n/2] + 1
  65     w,h = signal.freqz(d,1)
  66     return d/max(abs(h))
  67
  68 def c_dump(fond, harm1, other):
  69     f = fond * 32768
  70     f = map(int, f)
  71     h = harm1 * 32768
  72     h = map(int, h)
  73     o = other * 32768
  74     o = map(int, o)
  75     print "/* Apply the 3 filters on the bitfield. The results are stored in"
  76     print " * fil_fond, fil_harm1 and fil_other. */"
  77     print "static void apply_filters(void)"
  78     print "{"
  79     print "     fil_fond = 0;"
  80     print "     fil_harm1 = 0;"
  81     print "     fil_other = 0;"
  82     print
  83     for i in range(len(f)):
  84         print " if (bitfield & (1UL << %d)) {"%(i)
  85         print "         fil_fond += %d;"%(f[i])
  86         print "         fil_harm1 += %d;"%(h[i])
  87         print "         fil_other += %d;"%(o[i])
  88         print " }"
  89         print
  90     print "#ifdef HOST_VERSION"
  91     print "     fil_fond = saturate16(fil_fond);"
  92     print "     fil_harm1 = saturate16(fil_harm1);"
  93     print "     fil_other = saturate16(fil_other);"
  94     print "#endif"
  95     print "}"
  96
  97
  98 figure(3)
  99 fond = bandpass(n, Ffond)
 100 harm1 = bandpass(n, Fharm1)
 101 other = bandpass(n, Fother)
 102 c_dump(fond, harm1, other)
 103
 104 #Frequency response
 105 mfreqz(fond)
 106 mfreqz(harm1)
 107 mfreqz(other)
 108
 109 if len(sys.argv) >= 2:
 110     savefig(sys.argv[1])
 111 else:
 112     show()