频谱显示波纹，我们可以直观地量化为 ~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 Hz之间 每个GHz的大约17个振荡约为你的涟漪。
【讨论】：
问题源于您混淆了您正在测量的术语“频率”和数据的频率。
你想要的是纹波频率，实际上就是你数据的周期。
解决了这个问题，让我们看看如何修复您的 fft。
正如Dmitrii's answer 所指出的，您必须确定数据的采样频率，并去除 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 到一半。这将在this answer
上进一步探讨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')
剧情如下：
【讨论】：