Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

FFT low-pass filter

Be warned, this is a newbie question. I acquired some noisy data (a 1x200 pixel sclice from a grayscale image), for which I am trying to build a simple FFT low-pass filter. I do understand the general principle of the Fourier Transform, but I ran into trouble trying to implement it. My script works well on example data, but behaves in a strange manner on my data.

I think I must be mixing dimensions at some point, but it's been quite a few long hours and I cannot find where! I suspect that, because the output (please see below) of print(signal.shape) is different between the example and real data. Furthermore, scipy.fftpack.rfft(signal) seems to do nothing to my signal instead of computing the function in the frequency domain, as it should.

My script:

(will run out-of-the-box using example data, just by copy-pasting everything below #example data)

import cv2
import numpy as np    
from scipy.fftpack import rfft, irfft, fftfreq, fft, ifft
import matplotlib as mpl
import matplotlib.pyplot as plt
#===========================================
#GETTING DATA AND SETTING CONSTANTS
#===========================================
REACH = 100
COURSE = 180
CENTER = (cx, cy)
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
gray = cv2.equalizeHist(gray)
gray2 = gray.copy()

#drawing initial vector
cv2.line(gray, (cx, cy + REACH), (cx, cy - REACH), 0, 5)
cv2.circle(gray, (cx, cy + REACH), 10, 0, -1)
cv2.circle(gray, (cx, cy), REACH, 0, 5)

#flooding contour with white
cv2.drawContours(gray2, contours, index, 255, -1)

#real data
signal = gray2[(cy - REACH):(cy + REACH), (cx-0.5):(cx+0.5)]
time = np.linspace(0, 2*REACH, num=200)

#example data
time   = np.linspace(0,10,2000)
signal = np.cos(5*np.pi*time) + np.cos(7*np.pi*time)

#=============================================
#THE FFT TRANSFORM & FILTERING
#=============================================
#signal filtering
f_signal = rfft(signal)
W = fftfreq(signal.size, d=time[1]-time[0])


cut_f_signal = f_signal.copy()
cut_f_signal[(W>5)] = 0

cut_signal = irfft(cut_f_signal)

#==================================
#FROM HERE ITS ONLY PLOTTING
#==================================

print(signal.shape)
plt.figure(figsize=(8,8))

ax1 = plt.subplot(321)
ax1.plot(signal)
ax1.set_title("Original Signal", color='green', fontsize=16)

ax2 = plt.subplot(322)
ax2.plot(np.abs(f_signal))
plt.xlim(0,100)
ax2.set_title("FFT Signal", color='green', fontsize=16)

ax3 = plt.subplot(323)
ax3.plot(cut_f_signal)
plt.xlim(0,100)
ax3.set_title("Filtered FFT Signal", color='green', fontsize=16)

ax4 = plt.subplot(324)
ax4.plot(cut_signal)
ax4.set_title("Filtered Signal", color='green', fontsize=16)

for i in [ax1,ax2,ax3,ax4]:
    i.tick_params(labelsize=16, labelcolor='green')

plt.tight_layout()
plt.show()

The result on real data:

parameters:

signal = gray2[(cy - REACH):(cy + REACH), (cx-0.5):(cx+0.5)]
time = np.linspace(0, 2*REACH, num=200)

filtering parameter:

cut_f_signal[(W<0.05)] = 0

Output:

output of signal.shape is (200L, 1L)

real data output

The result on example data:

parameters:

signal = np.cos(5*np.pi*time) + np.cos(7*np.pi*time)
time   = np.linspace(0,10,2000)

filtering parameter:

cut_f_signal[(W>5)] = 0

Output:

output of signal.shape is (2000L,)

example data output

like image 737
Raoul Avatar asked Feb 20 '26 03:02

Raoul


1 Answers

So I began to think about that, and after a time I realized the stupidity in my ways. My base data is an image, and I take a slice of it to generate the above signal. So instead of trying to implement a less-than-satisfying home-brewed FFT filter to smoooth the signal, it is in fact much better and easier to smooth the image with one of the numerous battle-tested filters (gaussian, bilateral, etc.) available in equally numerous image libs (in my case, OpenCV)...

Thanks to the people that took the time to try and help!

like image 102
Raoul Avatar answered Feb 21 '26 15:02

Raoul



Donate For Us

If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!