过滤后的噪声信号 fft

Question

我想创建两个采样频率为 2.5G sa/s 的随机噪声信号，频率范围为 200kHz - 20Mhz，信号时间为 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:

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.

因此，最简单的方法是使用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()将零赫兹频率移动到阵列的中心。

此外，由于您使用的是实数信号，出于效率原因，您可能会考虑使用 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) 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()

过滤后的噪声信号 fft

Flitered noise signal fft

numpy

filtering

fft

noise

noise-generator