Filtering - multiplication in frequency domain

Question

I am trying to create a simple low pass filter, but I have got what I consider a surprising result when looking at the frequency response of a simple Butterworth filter.

I have copied much of the example below from this other post. I have added some code at the bottom of the script to compare the input and output spectra against the frequency response of the filter. I would expect that the output spectrum $\mathbf B$ should be the product of the input spectrum $\mathbf A$ and the frequency response $\mathbf H$: $$ \mathbf B = \mathbf H \mathbf A $$

However, the plot on below shows that the filter actually increases some low frequency components - see how the red line is above the green below around $4\textrm{ Hz}$.

Could anyone please explain why this is?

import numpy as np
from scipy.signal import butter, lfilter, freqz
import matplotlib.pyplot as plt
from scipy.fftpack import fft as fft
def butter_lowpass(cutoff, fs, order=5):
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = butter(order, normal_cutoff, btype='low', analog=False)
    return b, a

def butter_lowpass_filter(data, cutoff, fs, order=5):
    b, a = butter_lowpass(cutoff, fs, order=order)
    y = lfilter(b, a, data)
    return y


# Filter requirements.
order = 6
fs = 30.0       # sample rate, Hz
cutoff = 3.667  # desired cutoff frequency of the filter, Hz

# Get the filter coefficients so we can check its frequency response.
b, a = butter_lowpass(cutoff, fs, order)

# Plot the frequency response.
w, h = freqz(b, a, worN=8000)
plt.subplot(2, 1, 1)
plt.plot(0.5*fs*w/np.pi, np.abs(h), 'b')
plt.plot(cutoff, 0.5*np.sqrt(2), 'ko')
plt.axvline(cutoff, color='k')
plt.xlim(0, 0.5*fs)
plt.title("Lowpass Filter Frequency Response")
plt.xlabel('Frequency [Hz]')
plt.grid()


# Demonstrate the use of the filter.
# First make some data to be filtered.
T = 5.0         # seconds
n = int(T * fs) # total number of samples
t = np.linspace(0, T, n, endpoint=False)
# "Noisy" data.  We want to recover the 1.2 Hz signal from this.
data = np.sin(1.2*2*np.pi*t) + 1.5*np.cos(9*2*np.pi*t) + 0.5*np.sin(12.0*2*np.pi*t)

# Filter the data, and plot both the original and filtered signals.
y = butter_lowpass_filter(data, cutoff, fs, order)

plt.subplot(2, 1, 2)
plt.plot(t, data, 'b-', label='data')
plt.plot(t, y, 'g-', linewidth=2, label='filtered data')
plt.xlabel('Time [sec]')
plt.grid()
plt.legend()

plt.subplots_adjust(hspace=0.35)
plt.show()

def calculateFFT(time,signal):
    N=len(signal)      
    df=1/((time[-1]-time[0]))
    frequencies=[i*df for i in range(int(N/2.0))]
    fftValues = [2.0/N*abs(i) for i in fft(signal,N)[0:N/2.0] ]
    return frequencies,fftValues

plt.subplot(2, 1, 1)
originalfreqs,originalFFT=calculateFFT(t,data)
plt.plot(originalfreqs,originalFFT,"g",label="original")

filteredfreqs,filteredFFT=calculateFFT(t,y)
plt.plot(filteredfreqs,filteredFFT,"r",label="filtered")
plt.legend()

Answer 1

I gladly upvote this question, because this is how a good question should look like, and also, how I would expect students to verify their understanding of the stuff they learn. It is always good to be keen to understand the background of something that you want to use later on.

The problem you experience is the following: In principle you are right, that the convolution theorem implies that convolution in one domain is multiplication in the other domain. Hence you have

$$\mathcal{F}\{x(t)*h(t)\}=X(f)H(f)$$ where $\mathcal{F}$ denotes the Fourier Transform.

So, what is $x(t)$ and $h(t)$ in your system? $x(t)$ is your input signal (in blue), which is the sum of three sines, multiplied by a rectangular window. It is implicitely multiplied by a rectangular, because your time does not (in simulation, it cannot) reach from $-\infty$ to $\infty$. So, the spectrum of $X(f)$ is the convolution of a sinc function (from the rectangular window) with three Diracs (at the frequency of the sines). Clearly, this is not a purely discrete spectrum.

What is $h(t)$? It is the impulse response of the butterworth filter. The Butterworth filter has an infinitely long impulse response, hence your convolution product $x(t)*h(t)$ is also infinitely wide. So, in principle you cannot apply a finite (i.e. discrete) Fourier Transform (fft) and expect it be a similar to the continuous case.

So, what you see is reasonable. You can try zeropadding the input signal, such that you get (the main part of) the impulse response in your signal. However, even then, you will see the frequencies, since $X(f)$ is not discrete actually.

Answer 2

To expand on @Maximilian Matthé's answer, you can visualize the spectral leakage effects (convolution of by sinc function in the frequency-domain) by zero-padding the inputs used in calculateFFT. For example, the following function zero-pad the inputs to a length $k$ times longer than the original input (where in this instance $k=4$):

def calculateFFT(time,signal):
    k=4
    N=k*len(signal)      
    df=1/(k*(time[-1]-time[0]))
    frequencies=[i*df for i in range(int(N/2.0))]
    fftValues = [k*2.0/N*abs(i) for i in fft(signal,N)[0:int(N/2.0)] ]
    return frequencies,fftValues

Plotting the result you should see that the low-frequency components actually overlap (which indicates that the filtering does not actually increase the signal there):

So, where did those ripples around the spike came from? It turns out they were always there, but when computing the FFT you were obtaining the value of the spectrum at a discrete set of frequency values and those ripples happened to pass through exactly zero at those frequencies. If you had chosen a slightly different signal frequency that is not an exact multiple of the FFT frequency resolution (0.2Hz in your case), for example 1.25Hz instead of 1.2Hz, the sampling of the frequency spectrum would have looked quite different (due to the frequency sampling at different points in the oscillation of the ripples):