0

我想创建两个随机噪声信号，在 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()
4

1 回答 1

2

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

另一个问题是您应该注意的“fft.fftshift()”方法。您可以手动计算频率向量，但这有点繁琐，因为结果 fft 向量中元素的顺序。fft 的结果具有特殊的频率排列。从scipy.fft.fft() 文档：

频率项 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()将零赫兹频率移动到数组的中心。

此外，当您使用实数信号时，出于效率原因，您可能会考虑使用 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()
于 2020-05-29T16:00:49.650 回答