Python: Lowpass Filter with only numpy

I see that the comments of @Cris Luengo have already developed your solution into the right direction. The last thing you're missing now is that the spectrum you obtain from np.fft.fft is composed of the positive frequency components in the first half and the 'mirrored' negative frequency components in the second half.

If you now set all components beyond your bandlimit_index to zero, you're erradicating these negative frequency components. That explains the drop in signal amplitude of 6dB, you're eliminating half the signal power (+ as you already noticed every real signal has to have conjugate symmetric frequency spectrum). The np.fft.ifft function documentation (ifft documentation) explains the expected format quite nicely. It states:

"The input should be ordered in the same way as is returned by fft, i.e.,"

• a[0] should contain the zero frequency term,
• a[1:n//2] should contain the positive-frequency terms,
• a[n//2 + 1:] should contain the negative-frequency terms, in increasing order starting from the most negative frequency.

That's essentially the symmetry you have to preserve. So in order to preserve these components just set the components between bandlimit_index + 1 -> (len(fsig) - bandlimit_index) to zero.

    def low_pass_filter(adata: np.ndarray, bandlimit: int = 1000, sampling_rate: int = 44100) -> np.ndarray:
        # translate bandlimit from Hz to dataindex according to sampling rate and data size
        bandlimit_index = int(bandlimit * adata.size / sampling_rate)
    
        fsig = np.fft.fft(adata)
        
        for i in range(bandlimit_index + 1, len(fsig) - bandlimit_index ):
            fsig[i] = 0
            
        adata_filtered = np.fft.ifft(fsig)
    
        return np.real(adata_filtered)

or maybe if you want to be slightly more 'pythonic', you can also set the components to zero like this:

fsig[bandlimit_index+1 : -bandlimit_index] = 0

Here's a colab to walk through: https://colab.research.google.com/drive/1RR_9EYlApDMg4jAS2HuJIpSqwg5RLzGW?usp=sharing

A full FFT of a real signal contains two symmetrical halves. Each half will contain reflected (purely imaginary) complex conjugates of the other half: everything past Nyquist is not real information. If you set everything to zero past some frequency, you have to respect the symmetry.

Here is a contrived signal, sampled at 1kHz, and its FFT:

import numpy as np
from matplotlib import pyplot as plt

t = np.linspace(0, 10, 10000, endpoint=False)
y = np.sin(2 * np.pi * 2 * t) * np.exp(-0.5 * ((t - 5) / 1.5)**2)
f = np.fft.fftfreq(t.size, 0.001)
F = np.fft.fft(y)

plt.figure(constrained_layout=True)

plt.subplot(1, 2, 1)
plt.plot(t, y)
plt.title('Time Domain')
plt.xlabel('time (s)')

plt.subplot(1, 2, 2)
plt.plot(np.fft.fftshift(f), np.fft.fftshift(F.imag), label='imag')
plt.xlim([-3, 3])
plt.title('Frequency Domain')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Imginatry Component')

The frequency axis looks like this:

>>> f
array([ 0. ,  0.1,  0.2, ..., -0.3, -0.2, -0.1])

Notice that aside from the DC component (bin 0), the axis is symmtrical about the midpoint, with the highest (Nyquist) frequency in the middle. This is why I called fftshift to draw the plot: it rearranges the array to go from smallest to largest.

You don't need to restrict your inputs to integers most likely. A fractional band_limit is totally acceptable. Keep in mind that to convert frequency to index, you multiply by size and sampling frequency (divide by ram, rather than divide:

def low_pass_filter(data, band_limit, sampling_rate):
    cutoff_index = int(band_limit * data.size / sampling_rate)
    F = np.fft.fft(data)
    F[cutoff_index + 1 : -cutoff_index] = 0
    return np.fft.ifft(F).real

You still need to return the real component, because the FFT will always have some imaginary roundoff errors in the lowest couple of bits.

Here is a plot of the sample signal cut off above 2Hz:

y2 = low_pass_filter(y, 2, 1000)
f2 = np.fft.fftfreq(t.size, 0.001)
F2 = np.fft.fft(y2)

plt.figure(constrained_layout=True)

plt.subplot(1, 2, 1)
plt.plot(t, y2)
plt.title('Time Domain')
plt.xlabel('time (s)')

plt.subplot(1, 2, 2)
plt.plot(np.fft.fftshift(f2), np.fft.fftshift(F2.imag), label='imag')
plt.xlim([-3, 3])
plt.title('Frequency Domain')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Imginatry Component')

Remember how we said that only half of the FFT of a purely real (or purely complex) signal contains non-redundant information? Numpy respects this, and provides np.fft.rfft, np.fft.irfft, np.fft.rfftfreq to work with real-values signals. You can use this to write a simpler version of the filter, since there is no longer a symmetry constraint.

def low_pass_filter(data, band_limit, sampling_rate):
     cutoff_index = int(band_limit * data.size / sampling_rate)
     F = np.fft.rfft(data)
     F[cutoff_index + 1:] = 0
     return np.fft.irfft(F, n=data.size).real

The only caveat is that we must explicitly pass in n to irfft, otherwise the output size will always be even, regardless of input size.