Loading [MathJax]/extensions/MathMenu.js

Saturday, 9 May 2015

Filtering an ideal signal with FFT

Hi everyone! Remember that in the last article I wrote that you can use the FFT to clean a signal from background noise? Well here is an example of signal filtering.

We start by generating a signal and then add some random noise using the random number generator in numpy. Here is the noisy signal:

signal

Now, our signal is made up of two main frequencies: 20 and 30 Hz while the rest is mainly background noise. You can check this in the FFT of the signal:

FFT of the signal

Suppose we’d like to isolate the 20 Hz bit, that means we’d have to set to 0 each value in the Spectrum which is greater or lower than 20 Hz. Let’s block out everything outside 18 and 22 Hz and then apply the inverse FTT. Here is the signal we found (blue) compared to the exact signal we wanted to find (red). At the bottom of the graph there is the error between the two signals, notice that it is not that high, (at worst we missed 0.08, however on average we missed less).

final

Here is the python code I used to make this.

import matplotlib.pyplot as plt
import numpy as np
import scipy as sc
plt.style.use('ggplot')
np.random.seed(20)
#-------------------------------------------------------------------------------
# Set up
# Time
t = np.linspace(0,1,1000)
# Frequencies in the signal
f1 = 20
f2 = 30
# Some random noise to add to the signal
noise = np.random.random_sample(len(t))
# Complete signal
y = 2*np.sin(2*np.pi*f1*t+0.2) + 3*np.cos(2*np.pi*f2*t+0.3) + noise*5
# The part of the signal we want to isolate
y1 = 2*np.sin(2*np.pi*f1*t+0.2)
# FFT of the signal
F = sc.fft(y)
# Other specs
N = len(t) # number of samples
dt = 0.001 # inter-sample time difference
w = np.fft.fftfreq(N, dt) # list of frequencies for the FFT
pFrequency = np.where(w>=0)[0] # we only positive frequencies
magnitudeF = abs(F[:len(pFrequency)]) # magnitude of F for the positive frequencies
#-------------------------------------------------------------------------------
# Some functions we will need
# Plots the FFT
def pltfft():
plt.plot(pFrequency,magnitudeF)
plt.xlabel('Hz')
plt.ylabel('Magnitude')
plt.title('FFT of the full signal')
plt.grid(True)
plt.show()
# Plots the full signal
def pltCompleteSignal():
plt.plot(t,y,'b')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.title('Full signal')
plt.grid(True)
plt.show()
# Filter function:
# blocks higher frequency than fmax, lower than fmin and returns the cleaned FT
def blockHigherFreq(FT,fmin,fmax,plot=False):
for i in range(len(F)):
if (i>= fmax) or (i<=fmin):
FT[i] = 0
if plot:
plt.plot(pFrequency,abs(FT[:len(pFrequency)]))
plt.xlabel('Hz')
plt.ylabel('Magnitude')
plt.title('Cleaned FFT')
plt.grid(True)
plt.show()
return FT
# Normalising function (gets the signal in a scale from 0 to 1)
def normalise(signal):
M = max(signal)
normalised = signal/M
return normalised
#-------------------------------------------------------------------------------
# Processing
# Cleaning the FT by selecting only frequencies between 18 and 22
newFT = blockHigherFreq(F,18,22)
# Getting back the cleaned signal
cleanedSignal = sc.ifft(F)
# Error
error = normalise(y1) - normalise(cleanedSignal)
#-------------------------------------------------------------------------------
# Plot the findings
pltCompleteSignal() #Plot the full signal
pltfft() #Plot fft
plt.figure()
plt.subplot(3,1,1) #Subplot 1
plt.title('Original signal')
plt.plot(t,y,'g')
plt.subplot(3,1,2) #Subplot 2
plt.plot(t,normalise(cleanedSignal),label='Cleaned signal',color='b')
plt.plot(t,normalise(y1),label='Signal to find',ls='-',color='r')
plt.title('Cleaned signal and signal to find')
plt.legend()
plt.subplot(3,1,3) #Subplot 3
plt.plot(t,error,color='r',label='error')
plt.show()
view raw fft_py.py hosted with ❤ by GitHub

Have fun!

2 comments:

  1. Isn't this is what is called FIR filter? Only, you avoid the FFT -- > IFFT, just do convolution.

    ReplyDelete
    Replies
    1. Hi, the idea was to have a look at the signal in the frequency domain using the fast fourier transform. As far as I know, a FIR system/filter is a different thing.

      Delete