Last active
December 23, 2015 06:49
-
-
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 file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
"""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