滤波噪声信号fft

Question

我想创建两个随机噪声信号，它们在200kHz-20Mhz的频率范围内以2.5G sa / s采样，信号时间为5us，并对其进行fft计算，但是fft有问题。谢谢您的帮助，代码为

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()

Answer 1

为了运行它，我必须对您的代码进行一些更改。主要的是将fs1.fft（）更改为fft.fft（）。

另一个问题是您应该注意的'fft.fftshift（）'方法。您可以手动计算频率向量，但是由于所得fft向量中元素的顺序，这有些乏味。 fft的结果具有特殊的频率安排。从scipy.fft.fft() documentation：

频率项f = k / n在y [k]处发现。在y [n / 2]处，我们达到奈奎斯特频率，然后绕到负频率项。因此，对于8点变换，结果的频率为[0、1、2、3，-4，-3，-2，-1]。要重新安排fft输出以使零频分量居中，例如[-4，-3，-2，-1、0、1、2、3]，请使用fftshift。

因此，最简单的方法是使用scipy.fft.fftfreq()让scipy为您进行计算。如果要以自然方式绘制它，则应调用scipy.fft.fftshift()将零Hz频率移动到阵列的中心。

此外，由于您使用的是实信号，出于效率原因，您可能会考虑使用fft算法scipy.fft.rfft()的实数版本。输出不包括负频率，因为对于实数输入阵列，完整算法的输出始终是对称的。

请参阅下面的代码。

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) distribuited 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()