1
$\begingroup$

In order to understand the behaviour of the butterworth (lowpass) filter I generated a "spectrum" consisting of 1's. I then went to time domain (ifft), applied the filter and transformed back (fft). I chose the window (maximum) frequency to be one. So i expected that the windowed spectrum should be equal to the actual spectrum of ones. Instead I got almost the spectrum, but at some point I got a huge line (I expect it to be spectral leakage?!). Also, the values are slightly moved up... Unfortunately I haven't understood how to solve this problem.

Here a sample source code from which I started:

from scipy import signal
import numpy as np
import matplotlib.pyplot as plt

wavelength = range(10000)
spectrum = [1]*len(wavelength)
filter_order = 3 #2
Wn = 1
inter = np.fft.ifft(spectrum) 
b,a = signal.butter(filter_order, Wn, btype='low', analog=False)
inter_w = signal.filtfilt(b, a, inter, method='gust') # filtered
spec_w = np.fft.fft(inter_w) # filtered spectrum

plt.plot(wavelength, spectrum, label='spec')
plt.plot(wavelength, spec_w, label='spec_w')
plt.legend()
plt.show()

And here a small plot. As stated I don't understand, why "spec_w" is (i) bigger than one and (ii) there is a peak that shouldn't be there.

Filter_order 3: spectrum_order3 Filter_order 2: spectrum_order2

$\endgroup$

1 Answer 1

2
$\begingroup$

Everything you see are numerical artefacts. A "low pass" filter with Nyquist as cut-off frequency is simply a connection. However, you're asking for a recursive filter of order $2$ or $3$, so what you get is a filter with identical numerator and denominator polynomials (i.e., pole/zero cancellation). Check your filter coefficients, a and b should be identical.

So what you have is an unnecessarily complicated implementation of the identity operation, which results in numerical errors. The numerical errors are caused by the "filtering" operation (made worse by the double filtering implemented in filtfilt), and by the combination of FFT and IFFT.

Also note that you should plot the absolute value of spec_w, because due to numerical inaccuracies spec_w might have a small imaginary part and I don't know what your plot function does with a complex-valued vector.

$\endgroup$
4
  • $\begingroup$ Yes, I wanted to test whether this algorithm returns the input. I am working on spectra and need to apply these filters. I made this "dummy-test" to analyse the discretization errors... do you know any way to reduce the error? I actually plot the absolute value, check the line with the fft. $\endgroup$
    – famfop
    Jul 7, 2016 at 11:56
  • $\begingroup$ @famfop: What are "these filters"? Low pass filters with cut-off of $1$ (Nyquist) are pointless, and they also cause numerical errors due to pole-zero cancellation. If you use an actual low pass filter (cut-off < Nyquist) you won't have any pole-zero cancellation. $\endgroup$
    – Matt L.
    Jul 7, 2016 at 12:00
  • $\begingroup$ Yes, a low pass digital filter with a cutoff frequency of 1 is nonsense. $\endgroup$
    – Peter K.
    Jul 7, 2016 at 12:00
  • $\begingroup$ With these filters I meant butterworth, chebyshev and some others. Although I know it's not useful I thought it might return kind of a frequency response or something to visualize the functionality. Anyway, by filtering one of two sinusoids I basically got the point. Thanks for your effort :) $\endgroup$
    – famfop
    Jul 7, 2016 at 13:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.