1 from numpy import cos, sin, pi, absolute, arange
   2 from scipy.signal import kaiserord, lfilter, firwin, freqz
   3 from pylab import figure, clf, plot, xlabel, ylabel, xlim, ylim, title, grid, axes, show
   4 
   5 
   6 #------------------------------------------------
   7 # Create a signal for demonstration.
   8 #------------------------------------------------
   9 
  10 sample_rate = 100.0
  11 nsamples = 400
  12 t = arange(nsamples) / sample_rate
  13 x = cos(2*pi*0.5*t) + 0.2*sin(2*pi*2.5*t+0.1) + \
  14         0.2*sin(2*pi*15.3*t) + 0.1*sin(2*pi*16.7*t + 0.1) + \
  15             0.1*sin(2*pi*23.45*t+.8)
  16 
  17 
  18 #------------------------------------------------
  19 # Create a FIR filter and apply it to x.
  20 #------------------------------------------------
  21 
  22 # The Nyquist rate of the signal.
  23 nyq_rate = sample_rate / 2.0
  24 
  25 # The desired width of the transition from pass to stop,
  26 # relative to the Nyquist rate.  We'll design the filter
  27 # with a 5 Hz transition width.
  28 width = 5.0/nyq_rate
  29 
  30 # The desired attenuation in the stop band, in dB.
  31 ripple_db = 60.0
  32 
  33 # Compute the order and Kaiser parameter for the FIR filter.
  34 N, beta = kaiserord(ripple_db, width)
  35 
  36 # The cutoff frequency of the filter.
  37 cutoff_hz = 10.0
  38 
  39 # Use firwin with a Kaiser window to create a lowpass FIR filter.
  40 taps = firwin(N, cutoff_hz/nyq_rate, window=('kaiser', beta))
  41 
  42 # Use lfilter to filter x with the FIR filter.
  43 filtered_x = lfilter(taps, 1.0, x)
  44 
  45 #------------------------------------------------
  46 # Plot the FIR filter coefficients.
  47 #------------------------------------------------
  48 
  49 figure(1)
  50 plot(taps, 'bo-', linewidth=2)
  51 title('Filter Coefficients (%d taps)' % N)
  52 grid(True)
  53 
  54 #------------------------------------------------
  55 # Plot the magnitude response of the filter.
  56 #------------------------------------------------
  57 
  58 figure(2)
  59 clf()
  60 w, h = freqz(taps, worN=8000)
  61 plot((w/pi)*nyq_rate, absolute(h), linewidth=2)
  62 xlabel('Frequency (Hz)')
  63 ylabel('Gain')
  64 title('Frequency Response')
  65 ylim(-0.05, 1.05)
  66 grid(True)
  67 
  68 # Upper inset plot.
  69 ax1 = axes([0.42, 0.6, .45, .25])
  70 plot((w/pi)*nyq_rate, absolute(h), linewidth=2)
  71 xlim(0,8.0)
  72 ylim(0.9985, 1.001)
  73 grid(True)
  74 
  75 # Lower inset plot
  76 ax2 = axes([0.42, 0.25, .45, .25])
  77 plot((w/pi)*nyq_rate, absolute(h), linewidth=2)
  78 xlim(12.0, 20.0)
  79 ylim(0.0, 0.0025)
  80 grid(True)
  81 
  82 #------------------------------------------------
  83 # Plot the original and filtered signals.
  84 #------------------------------------------------
  85 
  86 # The phase delay of the filtered signal.
  87 delay = 0.5 * (N-1) / sample_rate
  88 
  89 figure(3)
  90 # Plot the original signal.
  91 plot(t, x)
  92 # Plot the filtered signal, shifted to compensate for the phase delay.
  93 plot(t-delay, filtered_x, 'r-')
  94 # Plot just the "good" part of the filtered signal.  The first N-1
  95 # samples are "corrupted" by the initial conditions.
  96 plot(t[N-1:]-delay, filtered_x[N-1:], 'g', linewidth=4)
  97 
  98 xlabel('t')
  99 grid(True)
 100 
 101 show()