使用 python 计算频谱的 FFT
Computing FFT of a spectrum using python
频谱显示的纹波我们可以直观地量化为 ~50 MHz 纹波。我正在寻找一种方法来计算这些波纹的频率，而不是通过目视检查数千个光谱。由于该函数在频域中，采用 FFT 会使它回到时域（如果我是正确的，则使用时间反转）。我们怎样才能得到这些涟漪的频率？
要为蓝色图获得合适的光谱，您需要做两件事：
- 正确计算频谱图（红色）的频率
- 消除数据中的偏差，使光谱较少受到低污染
频率。那是因为你感兴趣的是涟漪，而不是缓慢的波动。
请注意，当您计算 fft 时，您会得到复数值，其中包含有关每个频率的振荡幅度和相位的信息。在您的情况下，红色图应该是振幅谱（与相位谱相比）。为此，我们采用的绝对值
fft 系数。
此外，您使用 fft 获得的频谱是两侧对称的（因为信号是真实的）。您真的只需要一侧就可以了解纹波峰值频率在哪里。我已经在代码中实现了它。
使用你的数据后，我得到了以下结果：
import pandas as pd
import numpy as np
import pylab as plt
import plotly.graph_objects as go
from scipy import signal as sig
df = pd.read_csv("ripple.csv")
f = df.Frequency.to_numpy()
data = df.Data
data = sig.medfilt(data) # median filter to remove the spikes
fig = go.Figure()
fig.add_trace(go.Scatter(x=f, y=(data - data.mean())))
fig.update_layout(
xaxis_title="Frequency in GHz", yaxis_title="dB"
) # the blue plot with ripples
fig.show()
# Remove bias to get rid of low frequency peak
data_fft = np.fft.fft(data - data.mean())
L = len(data) # number of samples
# Compute two-sided spectrum
tssp = abs(data_fft / L)
# Compute one-sided spectrum
ossp = tssp[0 : int(L / 2)]
ossp[1:-1] = 2 * ossp[1:-1]
delta_freq = f[1] - f[0] # without this freqs computation is incorrect
freqs = np.fft.fftfreq(f.shape[-1], delta_freq)
# Use first half of freqs since spectrum is one-sided
plt.plot(freqs[: int(L / 2)], ossp, "r-") # the red plot
plt.xlim([0, 50])
plt.xticks(np.arange(0, 50, 1))
plt.grid()
plt.xlabel("Oscillations per frequency")
plt.show()
所以你可以看到有两个峰值：低频。 1 和 2 赫兹之间的振荡
并且您的纹波频率约为每 GHz 17 次振荡。
问题是因为您混淆了您正在测量的术语 'frequency' 和数据频率。
你要的是纹波频率，其实就是你数据的周期
说完这些，让我们看看如何修复您的 fft。
正如 所指出的，您必须确定数据的采样频率，并去除 FFT 结果中的低频分量。
要确定采样频率，您可以通过每个样本减去其前一个样本并计算平均值来确定采样周期。平均采样频率正好是它的倒数。
fs = 1 / np.mean(freq[1:] - freq[:-1])
对于高通滤波器，你可以使用巴特沃斯滤波器，this是一个很好的实现。
# Defining a high pass filter
def butter_highpass(cutoff, fs, order=5):
nyq = 0.5 * fs
normal_cutoff = cutoff / nyq
b, a = signal.butter(order, normal_cutoff, btype='high', analog=False)
return b, a
def butter_highpass_filter(data, cutoff, fs, order=5):
b, a = butter_highpass(cutoff, fs, order=order)
y = signal.filtfilt(b, a, data)
return y
接下来，绘制fft时，需要取其绝对值，这就是你所追求的。此外，由于它同时提供了正面和负面的部分，因此您可以只使用正面的部分。就 x 轴而言，它将是从 0 到采样频率的一半。
对此进行了进一步探讨
fft_amp = np.abs(np.fft.fft(amp, amp.size))
fft_amp = fft_amp[0:fft_amp.size // 2]
fft_freq = np.linspace(0, fs / 2, fft_amp.size)
现在，要确定纹波频率，只需获取 FFT 的峰值即可。您正在寻找的值（大约 50MHz）将是纹波峰值的周期（以 GHz 为单位），因为您的原始数据以 GHz 为单位。对于这个例子，它实际上是 57MHz 左右。
peak = fft_freq[np.argmax(fft_amp)]
ripple_period = 1 / peak * 1000
print(f'The ripple period is {ripple_period} MHz')
这是完整的代码，其中还绘制了数据。
import numpy as np
import pylab as plt
from scipy import signal as signal
# Defining a high pass filter
def butter_highpass(cutoff, fs, order=5):
nyq = 0.5 * fs
normal_cutoff = cutoff / nyq
b, a = signal.butter(order, normal_cutoff, btype='high', analog=False)
return b, a
def butter_highpass_filter(data, cutoff, fs, order=5):
b, a = butter_highpass(cutoff, fs, order=order)
y = signal.filtfilt(b, a, data)
return y
with open('ripple.csv', 'r') as fil:
data = np.genfromtxt(fil, delimiter=',', skip_header=True)
amp = data[:, 0]
freq = data[:, 1]
# Determine the sampling frequency of the data (it is around 500 Hz)
fs = 1 / np.mean(freq[1:] - freq[:-1])
# Apply a median filter to remove the noise
amp = signal.medfilt(amp)
# Apply a highpass filter to remove the low frequency components 5 Hz was chosen
# as the cutoff fequency by visual inspection. Depending on the problem, you
# might want to choose a different value
cutoff_freq = 5
amp = butter_highpass_filter(amp, cutoff_freq, fs)
_, ax = plt.subplots(ncols=2, nrows=1)
ax[0].plot(freq, amp)
ax[0].set_xlabel('Frequency GHz')
ax[0].set_ylabel('Intensity dB')
ax[0].set_title('Filtered signal')
# The FFT part is as follows
fft_amp = np.abs(np.fft.fft(amp, amp.size))
fft_amp = fft_amp[0:fft_amp.size // 2]
fft_freq = np.linspace(0, fs / 2, fft_amp.size)
ax[1].plot(fft_freq, 2 / fft_amp.size * fft_amp, 'r-') # the red plot
ax[1].set_xlabel('FFT frequency')
ax[1].set_ylabel('Intensity dB')
plt.show()
peak = fft_freq[np.argmax(fft_amp)]
ripple_period = 1 / peak * 1000
print(f'The ripple period is {ripple_period} MHz')
剧情如下：
频谱显示的纹波我们可以直观地量化为 ~50 MHz 纹波。我正在寻找一种方法来计算这些波纹的频率，而不是通过目视检查数千个光谱。由于该函数在频域中，采用 FFT 会使它回到时域（如果我是正确的，则使用时间反转）。我们怎样才能得到这些涟漪的频率？
要为蓝色图获得合适的光谱，您需要做两件事：
- 正确计算频谱图（红色）的频率
- 消除数据中的偏差，使光谱较少受到低污染 频率。那是因为你感兴趣的是涟漪，而不是缓慢的波动。
请注意，当您计算 fft 时，您会得到复数值，其中包含有关每个频率的振荡幅度和相位的信息。在您的情况下，红色图应该是振幅谱（与相位谱相比）。为此，我们采用的绝对值 fft 系数。
此外，您使用 fft 获得的频谱是两侧对称的（因为信号是真实的）。您真的只需要一侧就可以了解纹波峰值频率在哪里。我已经在代码中实现了它。
使用你的数据后，我得到了以下结果：
import pandas as pd
import numpy as np
import pylab as plt
import plotly.graph_objects as go
from scipy import signal as sig
df = pd.read_csv("ripple.csv")
f = df.Frequency.to_numpy()
data = df.Data
data = sig.medfilt(data) # median filter to remove the spikes
fig = go.Figure()
fig.add_trace(go.Scatter(x=f, y=(data - data.mean())))
fig.update_layout(
xaxis_title="Frequency in GHz", yaxis_title="dB"
) # the blue plot with ripples
fig.show()
# Remove bias to get rid of low frequency peak
data_fft = np.fft.fft(data - data.mean())
L = len(data) # number of samples
# Compute two-sided spectrum
tssp = abs(data_fft / L)
# Compute one-sided spectrum
ossp = tssp[0 : int(L / 2)]
ossp[1:-1] = 2 * ossp[1:-1]
delta_freq = f[1] - f[0] # without this freqs computation is incorrect
freqs = np.fft.fftfreq(f.shape[-1], delta_freq)
# Use first half of freqs since spectrum is one-sided
plt.plot(freqs[: int(L / 2)], ossp, "r-") # the red plot
plt.xlim([0, 50])
plt.xticks(np.arange(0, 50, 1))
plt.grid()
plt.xlabel("Oscillations per frequency")
plt.show()
所以你可以看到有两个峰值：低频。 1 和 2 赫兹之间的振荡 并且您的纹波频率约为每 GHz 17 次振荡。
问题是因为您混淆了您正在测量的术语 'frequency' 和数据频率。
你要的是纹波频率，其实就是你数据的周期
说完这些，让我们看看如何修复您的 fft。
正如
要确定采样频率，您可以通过每个样本减去其前一个样本并计算平均值来确定采样周期。平均采样频率正好是它的倒数。
fs = 1 / np.mean(freq[1:] - freq[:-1])
对于高通滤波器，你可以使用巴特沃斯滤波器，this是一个很好的实现。
# Defining a high pass filter
def butter_highpass(cutoff, fs, order=5):
nyq = 0.5 * fs
normal_cutoff = cutoff / nyq
b, a = signal.butter(order, normal_cutoff, btype='high', analog=False)
return b, a
def butter_highpass_filter(data, cutoff, fs, order=5):
b, a = butter_highpass(cutoff, fs, order=order)
y = signal.filtfilt(b, a, data)
return y
接下来，绘制fft时，需要取其绝对值，这就是你所追求的。此外，由于它同时提供了正面和负面的部分，因此您可以只使用正面的部分。就 x 轴而言，它将是从 0 到采样频率的一半。
fft_amp = np.abs(np.fft.fft(amp, amp.size))
fft_amp = fft_amp[0:fft_amp.size // 2]
fft_freq = np.linspace(0, fs / 2, fft_amp.size)
现在，要确定纹波频率，只需获取 FFT 的峰值即可。您正在寻找的值（大约 50MHz）将是纹波峰值的周期（以 GHz 为单位），因为您的原始数据以 GHz 为单位。对于这个例子，它实际上是 57MHz 左右。
peak = fft_freq[np.argmax(fft_amp)]
ripple_period = 1 / peak * 1000
print(f'The ripple period is {ripple_period} MHz')
这是完整的代码，其中还绘制了数据。
import numpy as np
import pylab as plt
from scipy import signal as signal
# Defining a high pass filter
def butter_highpass(cutoff, fs, order=5):
nyq = 0.5 * fs
normal_cutoff = cutoff / nyq
b, a = signal.butter(order, normal_cutoff, btype='high', analog=False)
return b, a
def butter_highpass_filter(data, cutoff, fs, order=5):
b, a = butter_highpass(cutoff, fs, order=order)
y = signal.filtfilt(b, a, data)
return y
with open('ripple.csv', 'r') as fil:
data = np.genfromtxt(fil, delimiter=',', skip_header=True)
amp = data[:, 0]
freq = data[:, 1]
# Determine the sampling frequency of the data (it is around 500 Hz)
fs = 1 / np.mean(freq[1:] - freq[:-1])
# Apply a median filter to remove the noise
amp = signal.medfilt(amp)
# Apply a highpass filter to remove the low frequency components 5 Hz was chosen
# as the cutoff fequency by visual inspection. Depending on the problem, you
# might want to choose a different value
cutoff_freq = 5
amp = butter_highpass_filter(amp, cutoff_freq, fs)
_, ax = plt.subplots(ncols=2, nrows=1)
ax[0].plot(freq, amp)
ax[0].set_xlabel('Frequency GHz')
ax[0].set_ylabel('Intensity dB')
ax[0].set_title('Filtered signal')
# The FFT part is as follows
fft_amp = np.abs(np.fft.fft(amp, amp.size))
fft_amp = fft_amp[0:fft_amp.size // 2]
fft_freq = np.linspace(0, fs / 2, fft_amp.size)
ax[1].plot(fft_freq, 2 / fft_amp.size * fft_amp, 'r-') # the red plot
ax[1].set_xlabel('FFT frequency')
ax[1].set_ylabel('Intensity dB')
plt.show()
peak = fft_freq[np.argmax(fft_amp)]
ripple_period = 1 / peak * 1000
print(f'The ripple period is {ripple_period} MHz')
剧情如下：