from numpy import cos, sin, pi, absolute, arange from scipy.signal import kaiserord, lfilter, firwin, freqz from pylab import figure, clf, plot, xlabel, ylabel, xlim, ylim, title, \ grid, show, subplot #------------------------------------------------ # Create a signal for demonstration. #------------------------------------------------ sample_rate = 500.0 nsamples = 800 t = arange(nsamples) / sample_rate x = 0.1*cos(2*pi*0.5*t) + 0.2*sin(2*pi*2.5*t+0.1) + \ 3.5*sin(2*pi*22*t) + 4*sin(2*pi*35*t + 0.1) + 1.1*cos(2*pi*40*t) + \ 1.25*sin(2*pi*60*t+.8) + 1.13*sin(2*pi*77.12*t + 0.12) +\ 0.3*cos(2*pi*95.1*t + .43) + 0.45*sin(2*pi*120*t+1) + \ 0.25*cos(2*pi*191*t + 2) #------------------------------------------------ # Create a FIR filter and apply it to x. #------------------------------------------------ # The Nyquist rate of the signal. nyq_rate = sample_rate / 2.0 # The desired width of the transition from pass to stop, # relative to the Nyquist rate. We'll design the filter # with a 5 Hz transition width. width = 10.0/nyq_rate # The desired attenuation in the stop band, in dB. ripple_db = 60.0 # Compute the order and Kaiser parameter for the FIR filter. N, beta = kaiserord(ripple_db, width) # The cutoff frequencies of the bandpass filter. low_cutoff_hz = 15.0 high_cutoff_hz = 50.0 #low_edge = low_cutoff_hz / nyq_rate #high_edge = high_cutoff_hz / nyq_rate # Use firwin with a Kaiser window to create a bandpassFIR filter. taps = firwin(N, [low_cutoff_hz, high_cutoff_hz], pass_zero=False, window=('kaiser', beta), scale=False, nyq=nyq_rate) # Use lfilter to filter x with the FIR filter. filtered_x = lfilter(taps, 1.0, x) #------------------------------------------------ # Plot the FIR filter coefficients. #------------------------------------------------ figure(1) plot(taps, 'bo-', linewidth=2) title('Filter Coefficients (%d taps)' % N) grid(True) #------------------------------------------------ # Plot the magnitude response of the filter. #------------------------------------------------ figure(2) clf() subplot(4,1,1) w, h = freqz(taps, worN=8000) plot((w/pi)*nyq_rate, absolute(h), linewidth=2) xlabel('Frequency (Hz)') ylabel('Gain') title('Frequency Response') ylim(-0.01, 1.01) grid(True) # Frequency response detail. subplot(4,1,2) plot((w/pi)*nyq_rate, absolute(h), linewidth=2) xlim(0, low_cutoff_hz) ylim(0, 0.0025) grid(True) # Frequency response detail. subplot(4,1,3) plot((w/pi)*nyq_rate, absolute(h), linewidth=2) xlim(low_cutoff_hz, high_cutoff_hz) ylim(0.9985, 1.0015) grid(True) # Frequency response detail. subplot(4,1,4) plot((w/pi)*nyq_rate, absolute(h), linewidth=2) xlim(high_cutoff_hz, high_cutoff_hz+10) ylim(0, 0.0025) grid(True) #------------------------------------------------ # Plot the original and filtered signals. #------------------------------------------------ # The phase delay of the filtered signal. delay = 0.5 * (N-1) / sample_rate figure(3) subplot(2, 1, 1) # Plot the original signal. plot(t, x) # Plot the filtered signal, shifted to compensate for the phase delay. plot(t-delay, filtered_x, 'r-') # Plot just the "good" part of the filtered signal. The first N-1 # samples are "corrupted" by the initial conditions. plot(t[N-1:]-delay, filtered_x[N-1:], 'g', linewidth=4) xlabel('t') grid(True) subplot(2, 1, 2) # Detail of the time series plot. # Plot the original signal. plot(t, x) # Plot the filtered signal, shifted to compensate for the phase delay. plot(t-delay, filtered_x, 'r-') # Plot just the "good" part of the filtered signal. The first N-1 # samples are "corrupted" by the initial conditions. plot(t[N-1:]-delay, filtered_x[N-1:], 'g', linewidth=4) xlim(0.75, 1.) xlabel('t') grid(True) show()