Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
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
You can’t perform that action at this time.