从真实数据集中确定并过滤掉低频成分

Question

我正在尝试使用 FFT 从信号中滤除低频分量，并在真实的数据集（加利福尼亚州每小时的电力需求）中保留高频分量。到目前为止我已经尝试过：

X = fft(df['demand'])
y_fft_shift = fftshift(X)
n = 20 # number of low freq components to be removed (selected randomly) 
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)

# compare original signal vs filtered signal 
plt.figure(figsize = (25, 6))
plt.plot(df['demand'],'b') #df['hour'], 
plt.plot(abs(y_ifft.real),'r')
plt.xlabel('Datetime')
plt.ylabel('demand')
plt.title('original vs filtered signal')
plt.xticks(rotation=25) 
plt.show()

我不确定（a）我的实现是否正确以及（b）从离散傅里叶逆变换获得的结果是否是预期的结果。例如，如果我不采取

abs(y_ifft.real)

，我会得到负值。

在将第二种方法应用于真实数据集之前，我在合成信号上尝试了以下两种方法。

from scipy.fftpack import fft, ifft, fftfreq
sr = 2000 # sampling rate
ts = 1.0/sr # sampling interval
t = np.arange(0,1,ts)

#generate a signal
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) #FT of original signal
freq = fftfreq(len(x), d=1.0/len(x)) #compute freq.
# define the cut-off frequency
cut_off = 4.5

# high-pass filter by assign zeros to the 
# FFT amplitudes where the absolute 
# frequencies smaller than the cut-off 
sig_fft_filtered[np.abs(freq) < cut_off] = 0

# get the filtered signal in time domain
filtered = ifft(sig_fft_filtered)

我将上面的输出与下面的代码进行了比较，其标准是我只想删除最低的四个频率分量：

y = fft(x, axis=0)
y_fft_shift = fftshift(y)
n = 4
m = len(x)/2
# y_fft_shift[m-n+1:m+n+1] = 0 
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)

P.S.：SO 上的许多答案帮助我理解使用 FFT 以及合成数据进行图像去噪的概念，但在真实数据集上应用 FFT 却找不到太多

参考：使用 scipy/numpy 在 python 中进行图像处理的高通滤波器

使用的数据集：链接到上面使用的电力需求数据集

Answer 1

不要使用复杂版本的FFT；使用正品。

import matplotlib
import numpy as np
from matplotlib import pyplot as plt

sr = 2_000  # sampling rate
ts = 1/sr  # sampling interval
t = np.arange(0, 1, ts)

# generate a signal
freq = 1
y = 3*np.sin(2*np.pi*freq*t)

freq = 4
y += np.sin(2*np.pi*freq*t)

freq = 7
y += 0.5*np.sin(2*np.pi*freq*t)

fft = np.fft.rfft(y)  # FT of original signal
f = np.fft.rfftfreq(len(y), d=ts)  # compute freq.

# define the cut-off frequency
f_cutoff = 4.5
i_cutoff = round(f_cutoff/sr * len(t))

# high-pass filter by assign zeros to the
# FFT amplitudes where the absolute
# frequencies smaller than the cut-off
fft_filtered = fft.copy()
fft_filtered[:1 + i_cutoff] = 0

# get the filtered signal in time domain
y_filtered = np.fft.irfft(fft_filtered)

matplotlib.use('TkAgg')

ax_t: plt.Axes
ax_f: plt.Axes
fig, (ax_t, ax_f) = plt.subplots(nrows=2, ncols=1)

ax_t.plot(t, y, label='unfiltered')
ax_t.plot(t, y_filtered, label='filtered')
ax_t.set_xlabel('t (s)')
ax_t.set_ylabel('y')
ax_t.legend()

ax_f.loglog(f, np.abs(fft), label='unfiltered')
ax_f.loglog(f, np.abs(fft_filtered), label='filtered')
ax_f.set_xlabel('f (Hz)')
ax_f.set_ylabel('|y|')
ax_f.legend()

plt.show()