scipy 从真实的数据集中确定和过滤低频分量

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

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)

型

的
使用的数据集：上面使用的link to electricity demand dataset
附言：关于SO的许多答案帮助我理解了使用FFT以及合成数据进行图像去噪的概念，但无法找到将FFT应用于真实的数据集的更多信息
参考号：High Pass Filter for image processing in python by using scipy/numpy
How to interpret the output of scipy.fftpack.fft?的
How to interpret the results of the Discrete Fourier Transform (FFT) in Python的
Understanding FFT output in python的

不要使用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()

字符串
x1c 0d1x的数据

scipy 从真实的数据集中确定和过滤低频分量

1条答案

相关问题

热门标签

最新问答