When designing low/high/bandpass filters I encountered the problem that they do not work as expected (see code and output). I want to isolate one of the three frequency peaks (the one with the lowest frequency with lowpass filter, etc.). However, the wrong peak is isolated... Does anyone know what my mistake is?
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
def create_single_freq_state(N, w):
"""
Creates a state that only contains the frequencies in w. Outputs array with N entries.
"""
res = np.exp(2.0*np.pi*w[:,np.newaxis]*np.arange(N)/N*1j)
return np.sum(res, axis=0)
def band_pass_filter(f1, f2, numtaps, window="hamming"):
fil = signal.firwin(numtaps, [f1, f2], pass_zero=False, fs=2*numtaps, window=window, scale=False)
return fil
def high_pass_filter(f, numtaps, window="hamming"):
fil = signal.firwin(numtaps, f, pass_zero=False, fs=2*numtaps, window=window, scale=False)
return fil
def low_pass_filter(f, numtaps, window="hamming"):
fil = signal.firwin(numtaps, f, pass_zero=True, fs=2*numtaps, window=window, scale=False)
return fil
def plot_freq_response(fil):
fs = 2*len(fil)
w, h = signal.freqz(fil, 1, worN=2000)
plt.plot((fs * 0.5 / np.pi) * w, abs(h), label="Hamming window")
plt.show()
def plot_func(f):
x = np.arange(1, len(f)+1)
plt.plot(x, f, color ="red")
plt.show()
#create state with several frequencies (those contained in the variable ws)
N = 60
ws = np.array([15, 30, 50])
f = create_single_freq_state(N, ws)
numtaps = 60-1
freq1 = 19.9
freq2 = 40.0
#choose filter
fil = low_pass_filter(freq1, numtaps)
#fil = band_pass_filter(freq1, freq2, numtaps)
#fil = high_pass_filter(freq2, numtaps)
#plot signal in fourierspace
plot_func(np.fft.fft(f))
#plot frequency response
plot_freq_response(fil)
#plot result of the application of the filter
plot_func(np.absolute(np.fft.fft(np.convolve(fil, f, mode="same"))))
Output:The output of my code. The wrong frequency is isolated.
Two of your frequencies are so high that they appear as DC and a low frequency (this is called aliasing, and is a consequence of the Nyquist sampling theorem). The fact that you're working with a complex input signal, and that your signal is short relative to your filter, make this hard to see.
I changed ws
to be ws = np.array([7.5, 15, 25])
, and changed your first call to plot_func
to be plot_func(np.abs(np.fft.fft(f)))
(because one of the peaks is entirely in the complex part of the transform). Now the input and output graphs look as follows: