Unexpected frequency components after applying bandpass filter in python

Question

When I apply a bandpass filter in python to a pure sinusoidal signal, it seems that after filtering I get different frequency components across the spectrum (see orange line below). It looks like leakage, however plotting of the FFT should be "clean", i.e. the frequency of the sinusoid lands on a single bin. Why does the output signal not also land on a single bin at 10MHz, since the filter is supposed to be linear? Can I somehow avoid this effect?

import numpy as np
import math
import matplotlib.pyplot as plt
import scipy
from scipy.fftpack import fft

# -----------------------------
# Parameters
# -----------------------------
f_samp = 100e6
f_in = 10e6
N = int((f_samp / f_in) * 128)

# -----------------------------
# Signals
# -----------------------------
t = np.linspace(0, float((1 / f_samp) * (N - 1)), num=int(N), dtype=float)
samples = np.sin(2 * math.pi * f_in * t) #+ np.random.normal(0, 1, size=N)

freq_axis = np.fft.fftfreq(N, 1/f_samp)
freq_axis = freq_axis[0:N//2]
freq_axis = np.divide(freq_axis, 1e6) # Convert to MHz

# -----------------------------
# Apply filter
# -----------------------------
[b_coeffs, a_coeffs] = scipy.signal.butter(1, np.array([9e6, 11e6]), fs=f_samp, btype='bandpass')
samples_filt = scipy.signal.lfilter(b_coeffs, a_coeffs, samples)

sig_lp_filt_imp = np.concatenate(([1], np.zeros(len(t) - 1)))
sig_lp_filt_imp_resp = scipy.signal.lfilter(b_coeffs, a_coeffs, sig_lp_filt_imp)

# -----------------------------
# Plotting FFT
# -----------------------------
fft_lp_filt = 20*np.log10(np.abs(fft(sig_lp_filt_imp_resp)))[0:N//2]
fft_samples_filt = 20*np.log10(np.divide(np.abs(fft(samples_filt)), N/2))[0:N//2]
fft_samples = 20*np.log10(np.divide(np.abs(fft(samples)), N/2))[0:N//2]

plt.figure()
plt.scatter(freq_axis, fft_samples, label = 'FFT Samples', marker = "o")
plt.scatter(freq_axis, fft_samples_filt, label = 'FFT Samples Filt', marker = "x")
plt.plot(freq_axis, fft_lp_filt, label = 'FFT BP Filter', color='green')
plt.grid(True)
plt.ylim(-80, 10)
plt.xlabel('Frequency [MHz]')
plt.ylabel('[dB]')

plt.show()

Answer 1

When you "turn on" filtering at time zero, the filtered output starts with a transient

plt.figure('signal filtered')
plt.plot(t, samples_filt)

To make the graph this transient more pronounced, set the number of waveform periods to a value less than 128:

N = int((f_samp / f_in) * 32)

Here the transient is noticeable up to 0.5us.

This transient contributes into frequency bins different than the central frequency 10MHz.

To suppress this contribution in the plot, you can increase the periodic steady state share of waveform periods:

N = int((f_samp / f_in) * 2048)

The resulting FFT graphs is

With 128 periods we have -55/-37/-70dB, with 2048, -77/-62/-95dB (low freq/near center/high freq end).

If you want to do a "pure" convolution of a pure sine signal and a filter transfer function "unmarred" by transients, you should use cyclic convolution.