基于FFT对加州电力需求数据集保留高频成分的实现疑问

一、核心疑问解答 #

(a) 真实数据集的实现是否正确 #

你的代码整体流程符合频域高通滤波的标准逻辑：FFT变换→频谱移位→置零低频分量→逆移位→逆FFT，但存在几个细节需要优化：

1. 索引处理：m = len(X)/2如果遇到信号长度为奇数的情况会得到浮点数，建议改为m = len(X) // 2（整数除法）避免索引错误；同时要加边界判断，防止m-n出现负索引：

   start_idx = max(0, int(m - n))
   end_idx = min(len(y_fft_shift), int(m + n + 1))
   sig_fft_filtered_img[start_idx:end_idx] = 0
   

2. 低频分量选择：n=20是随机选取的，这种方式是固定移除中心附近的N对低频分量（fftshift后中心为0频，左右对应正负低频），和基于频率阈值的高通滤波逻辑不同。如果你的需求是去掉某一频率以下的所有分量，这种方式不够严谨；但如果只是想移除最“慢”的趋势分量，这个实现是可行的。

(b) 逆FFT结果的负值问题 #

• 不取abs(y_ifft.real)出现负值是完全正常的：高通滤波本质是提取原始信号减去低频趋势后的高频波动部分，这部分信号会围绕0上下波动（因为低频分量包含了信号的整体基准线）。
• 你之前取abs(y_ifft.real)是错误操作，会完全扭曲信号的真实形态，正确的做法是直接使用y_ifft.real（原始信号是实数，逆FFT结果的虚部是数值精度误差，可忽略）。

二、两种合成信号方案对比 #

方案1：基于截止频率的高通滤波 #

from scipy.fftpack import fft, ifft, fftfreq
sr = 2000 # 采样率
ts = 1.0/sr # 采样间隔
t = np.arange(0,1,ts)

# 生成合成信号
freq = 1. 
x = 3*np.sin(2*np.pi*freq*t)
freq = 4.
x += np.sin(2*np.pi*freq*t)
freq = 7.   
x += 0.5* np.sin(2*np.pi*freq*t)

y = fft(x, axis=0) # 原始信号FFT
freq = fftfreq(len(x), d=ts) # 计算实际频率（注意这里d应该用采样间隔ts，不是1.0/len(x)）
cut_off = 4.5 # 截止频率

# 高通滤波：置零低于截止频率的分量
sig_fft_filtered = y.copy()
sig_fft_filtered[np.abs(freq) < cut_off] = 0

# 逆FFT得到滤波后信号
filtered = ifft(sig_fft_filtered)

• 逻辑严谨：通过fftfreq计算每个FFT分量对应的实际频率，精准过滤低于阈值的低频，适合明确频率范围需求的场景。
• 注意修正：原代码中fftfreq的d参数应该用采样间隔ts，否则频率计算错误。

方案2：移除固定数量低频分量 #

y = fft(x, axis=0)
y_fft_shift = fftshift(y)
n = 4
m = len(x) // 2
sig_fft_filtered_img = y_fft_shift.copy()
sig_fft_filtered_img[int(m-n):int(m+n+1)] = 0
y_ifft_shift = ifftshift(sig_fft_filtered_img)
y_ifft = ifft(y_ifft_shift)

• 操作简单：直接移除频谱中心的N个低频分量，无需计算实际频率，适合快速去除整体趋势的场景，但无法对应到具体的频率阈值。

三、真实数据集优化建议 #

如果需要更贴合业务需求的高通滤波（比如去掉日周期以下的趋势），建议参考方案1的逻辑，结合电力数据的小时采样特性调整：

from scipy.fftpack import fft, ifft, fftfreq

# 采样间隔为1小时，所以d=1（单位：小时）
freq = fftfreq(len(df['demand']), d=1)
# 去掉周期大于24小时的分量（对应频率 < 1/24 Hz）
cut_off = 1/24

X = fft(df['demand'])
X_filtered = X.copy()
X_filtered[np.abs(freq) < cut_off] = 0
y_ifft = ifft(X_filtered)

# 可视化对比
plt.figure(figsize=(25, 6))
plt.plot(df['demand'], 'b', label='原始信号')
plt.plot(y_ifft.real, 'r', label='滤波后信号')
plt.xlabel('Datetime')
plt.ylabel('demand')
plt.title('原始信号 vs 高通滤波信号')
plt.xticks(rotation=25)
plt.legend()
plt.show()

内容的提问来源于stack exchange，提问作者Dhanya