I wanted to create two random noise signals sampled 2.5G sa/s in frequency range 200kHz - 20Mhz, time of signal 5us and calculate fft of it but I have a problem with fft. Thank you for help, the code is
import numpy as np
import matplotlib.pyplot as plot
from scipy import signal
from scipy import fft
import pandas as pd
t = np.arange(0, 5e-6, 4e-10)
s1 = 1e-8*np.random.normal(0, 1, 12500)
s2 = 1e-8*np.random.normal(0, 1, 12500)
sos1 = signal.butter(N=10, Wn=[200000, 20000000], btype='band', fs=2.5e9, output='sos')
sos2 = signal.butter(N=10, Wn=[200000, 20000000], btype='band', fs=2.5e9, output='sos')
fs1 = signal.sosfilt(sos1, s1)
fs2 = signal.sosfilt(sos2, s2)
f1 = abs(fs1.fft())
f2 = abs(fs2.fft())
ax1 = plot.subplot(311)
plot.plot(t, fs1, t, fs2)
#ax1.set_xlim([0, 5e-6])
plot.xlabel('Time (s)')
plot.ylabel('Current (A)')
ax2 = plot.subplot(312)
plot.plot(f1, f2)
plot.xlabel('Frequency (Hz)')
plot.ylabel('Current (A)')
plot.show()
I had to do some changes to your code in order to run it. The main one was to change fs1.fft() to fft.fft().
Another issue is the 'fft.fftshift()' method that you should be aware of. You can calculate the frequency vector by hand, but this is somewhat tedious because of the order of the elements in the resulting fft vector. The result of the fft has a peculiar frequency arrangement. From the scipy.fft.fft() documentation:
The frequency term f=k/n is found at y[k]. At y[n/2] we reach the Nyquist frequency and wrap around to the negative-frequency terms. So, for an 8-point transform, the frequencies of the result are [0, 1, 2, 3, -4, -3, -2, -1]. To rearrange the fft output so that the zero-frequency component is centered, like [-4, -3, -2, -1, 0, 1, 2, 3], use fftshift.
So, the easiest way is to use scipy.fft.fftfreq() to let scipy do the calculation for you. If you want to plot it in a natural way, then you should call scipy.fft.fftshift() to shift the zero Hz frequency to the center of the array.
Also, as you are using real signals, for efficiency reasons you might consider using the real number version of the fft algorithm scipy.fft.rfft(). The output does not include the negative frequencies, since for real input arrays the output of the full algorithm is always symmetric.
Please see the code below.
import matplotlib
matplotlib.use('Qt5Agg')
import numpy as np
import matplotlib.pyplot as plot
from scipy import signal
from scipy import fft
import pandas as pd
sampling_freq_Hz = 2.5e9
sampling_period_s = 1 / sampling_freq_Hz
signal_duration_s = 5.0e-6
wanted_number_of_points = signal_duration_s / sampling_period_s
f_low_Hz = 200.0e3
f_high_Hz = 20.0e6
msg = f'''
Sampling frequency: {sampling_freq_Hz} Hz
Sampling period: {sampling_period_s} s
Signal duration: {signal_duration_s} s
Wanted number of points: {wanted_number_of_points}
Lower frequency limit {f_low_Hz}
Upper frequency limit: {f_high_Hz}
'''
print(msg)
# Time axis
time_s = np.arange(0, signal_duration_s, sampling_period_s)
real_number_of_points = time_s.size
print(f'Real number of points: {real_number_of_points}')
# Normal(0,sigma^2) distributed random noise
sigma_2 = 1.0e-8
s1 = np.random.normal(0, sigma_2, real_number_of_points)
s2 = np.random.normal(0, sigma_2, real_number_of_points)
# Since both filters are equal, you only need one
sos1 = signal.butter(N=10, Wn=[f_low_Hz, f_high_Hz], btype='band', fs=sampling_freq_Hz, output='sos')
#sos2 = signal.butter(N=10, Wn=[f_low_Hz, f_high_Hz], btype='band', fs=sampling_freq_Hz, output='sos')
# Do the actual filtering
filtered_signal_1 = signal.sosfilt(sos1, s1)
filtered_signal_2 = signal.sosfilt(sos1, s2)
# Absolute value
f_1 = abs(fft.fft(filtered_signal_1))
f_2 = abs(fft.fft(filtered_signal_2))
freqs_Hz = fft.fftfreq(time_s.size, sampling_period_s)
# Shift the FFT for understandable plotting
f_1_shift = fft.fftshift(f_1)
f_2_shift = fft.fftshift(f_2)
freqs_Hz_shift = fft.fftshift(freqs_Hz)
# Plot
ax1 = plot.subplot(311)
ax1.plot(time_s, filtered_signal_1, time_s, filtered_signal_2)
#ax1.set_xlim([0, 5e-6])
ax1.set_xlabel('Time (s)')
ax1.set_ylabel('Current (A)')
ax2 = plot.subplot(313)
ax2.plot(freqs_Hz_shift, f_1_shift, freqs_Hz_shift, f_2_shift)
ax2.set_xlabel('Frequency (Hz)')
ax2.set_ylabel('Current (A)')
plot.show()