Skip to content

Instantly share code, notes, and snippets.

@MartinNowak
Created April 28, 2022 10:31
Show Gist options
  • Save MartinNowak/445074f87f26776c5b50430da085b775 to your computer and use it in GitHub Desktop.
Save MartinNowak/445074f87f26776c5b50430da085b775 to your computer and use it in GitHub Desktop.
FIR low-pass filter
#!python
from numpy import cos, sin, pi, absolute, arange
from scipy.signal import windows, kaiserord, lfilter, firwin, freqz
from pylab import figure, clf, plot, xlabel, ylabel, xlim, ylim, title, grid, axes, show
N = window_size_in_days = 20
#------------------------------------------------
# Create a signal for demonstration.
#------------------------------------------------
sample_rate = 1.0
nsamples = 400
t = arange(nsamples) / sample_rate
x = cos(2*pi*0.5*t) + 0.2*sin(2*pi*2.5*t+0.1) + \
0.2*sin(2*pi*15.3*t) + 0.1*sin(2*pi*16.7*t + 0.1) + \
0.1*sin(2*pi*23.45*t+.8)
#------------------------------------------------
# Create a FIR filter and apply it to x.
#------------------------------------------------
# The Nyquist rate of the signal.
nyq_rate = sample_rate / 2.0
# The cutoff frequency of the filter.
tau = 10 # days
cutoff_hz = sample_rate / tau
# Use firwin with a ~~Kaiser~~ window to create a lowpass FIR filter.
window = 'hamming' # default sharper low-pass
window = 'exponential' # closer to Croston's method with single-pole IIR low-pass
taps = firwin(numtaps=window_size_in_days, cutoff=cutoff_hz, window=window)
# 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()
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.05, 1.05)
grid(True)
# Upper inset plot.
ax1 = axes([0.42, 0.6, .45, .25])
plot((w/pi)*nyq_rate, absolute(h), linewidth=2)
xlim(0,8.0)
ylim(0.9985, 1.001)
grid(True)
# Lower inset plot
ax2 = axes([0.42, 0.25, .45, .25])
plot((w/pi)*nyq_rate, absolute(h), linewidth=2)
xlim(12.0, 20.0)
ylim(0.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)
# 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)
show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment