Skip to content

Instantly share code, notes, and snippets.

@moorepants
Last active December 23, 2015 06:49
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save moorepants/6597067 to your computer and use it in GitHub Desktop.
Save moorepants/6597067 to your computer and use it in GitHub Desktop.
Fixed Bode magnitude x values in the plot (corrected the units) and cleaned up the plot a bit.
"""This script makes a plot showing the the difference in filtering a signal
with a Butterworth filter if you select the digital cutoff frequency ratio
with respect to incorrect Nyquist frequency values. The cutoff frequency for
digital filters should be `cutoff/nyquist` where the `cutoff` is the desired
cutoff frequency for the filter and `nyquist` is the Nyquist frequency of
the data you are filtering, i.e. half the sample rate."""
from numpy import sin, cos, pi, array, absolute, ones, log10, linspace
from numpy.random import normal
from scipy.signal import butter, freqz, filtfilt
import matplotlib.pyplot as plt
# The noisy signal.
duration = 3.0 # sec
actual_sample_rate = 100.0 # Hz (samples / sec)
num_samples = int(duration * actual_sample_rate) + 1
t = linspace(0.0, duration, num=num_samples)
x = (1.0 * 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.5) +
0.1 * sin(2 * pi * 16.7 * t + 0.1) +
0.1 * sin(2 * pi * 23.45 * t + 0.8) +
normal(scale=0.1, size=t.shape))
fig, axes = plt.subplots(2, 1)
fig.set_size_inches(8, 6)
axes[1].plot(t, x, 'k.', label='Actual Signal')
butterworth_order = 2
desired_cutoff_hz = 6.0
# The correct Nyquist frequency is 50 Hz. But, let's choose some Nyquist
# frequencies other than and including the correct one.
sample_rates = array([50.0, 100.0, 200.0, 300.0])
nyquist_frequencies = sample_rates / 2.0
for sample_rate, nyquist in zip(sample_rates, nyquist_frequencies):
b, a = butter(butterworth_order, desired_cutoff_hz / nyquist)
w, h = freqz(b, a) # rad / sample, mag(x)
axes[0].semilogx(w * sample_rate / 2.0 / pi, # Hz
20.0 * log10(absolute(h)), # dB
label="Nyquist Frequency = {} Hz".format(nyquist))
xfiltered = filtfilt(b, a, x)
axes[1].plot(t, xfiltered,
label="Nyquist Frequency = {} Hz".format(nyquist))
axes[0].set_xlabel('Frequency [Hz]')
axes[0].set_ylabel('Amplitude [dB]')
axes[0].set_title('Magnitude plot of the Butterworth filter.')
axes[0].set_xlim([1.0, 10.0])
axes[0].set_ylim([-20.0, 10.0])
ylim = axes[0].get_ylim()
axes[0].plot(desired_cutoff_hz * ones(2),
[ylim[0] - 1.0, ylim[1] + 1.0], 'k--',
label='Cutoff Frequency: {} hz'.format(desired_cutoff_hz))
axes[0].plot([0.5, 10.5], -3.0 * ones(2), 'k', label='-3 dB')
axes[0].legend(fontsize=8, loc=3)
axes[1].set_xlabel('Time [s]')
axes[1].set_ylabel('Amplitude')
axes[1].set_title('Time series of the unfiltered and filtered data.')
axes[1].legend(fontsize=8)
plt.tight_layout()
fig.savefig('nyquist-error.png', dpi=300)
fig.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment