I am trying to use FFT to filter out low frequency components from a signal and retain high frequency components in a real dataset (hourly electricity demand in California). I have tried this so far:
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()
I am not sure whether (a) my implementation is correct, and (b) the results obtained from inverse discrete fourier transform is the expected results. For instance, if I don't take abs(y_ifft.real)
I get negative values.
I tried the following two approaches on synthetic signal before applying the second approach to the real dataset.
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)
I compared the output from the above with the below code with a criterion that I want to remove only lowest four frequency component:
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)
Dataset used: link to electricity demand dataset used above
P.S.: Many answers on SO helped me to understand the concept on image denoising using FFT as well as on synthetic data but couldn't find much on applying FFT on real dataset
Ref: 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
Don't use the complex versions of the FFT; use the real ones.
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()