Fourier Transform and Filtering

The focus of this post is on filtering frequencies and why we need it. And I’m sure you can find dozens of good tutorials on how to take Fourier transform of your data/signals

Introduction:

As a person who studies computational neuroscience, there’s a great deal of my time spending on brain signals and then processing those signals. If we can assume that signals are the sum of some pure sin waves with different frequencies (which in the majority of topics we can) then we can say that Fourier-transform is the building block of signal processing. Fourier-transform can tell from the signal that what sin waves with what frequencies are caused to create that particular signal. The noise in signals shows itself in the higher frequencies of sin waves. By these definitions, we can see that Fourier-transform can help us within at least two useful tools.

• It can give us the tool to study each oscillation in the signal individually.
• It can give us the tool to remove the noise from the original signal (The noise in signal shows itself with the higher frequencies in the frequency domain).

These above processings are categorized under the definition of filtering. To understand how Fourier-transform works first, we need to know the concepts of dot-product and convolutions.

Dot-product:

Dot-product which is the basic step of convolution is a very simple equation. To compute the dot product we simply multiply each element of one vector (array) by the corresponding element of the other vector:

Here n is the the length of the array and i is the index number of array a and b.

Convolution:

Convolution is an extension of a dot-product, in which the dot-product is computed repeatedly over time. We have two vectors in convolution, One is the signal itself, and the other is the kernel.

As you can see the convolution of the signal with the kernel, smooth out the sharp edges of the signal.

Fourier-transform:

Now that we understand how convolution works, we can understand how Fourier-transform works as well. Because Fourier transform is basically the same thing as convolution is, but with only a few differences.

• Instead of the simple line kernel, in Fourier transform the kernel is a sin wave with a specific frequency
• Instead of just only one kernel, in Fourier transform we use many kernels (sine waves with different frequencies) to cover a range of frequencies.

The above figure is just a simple convolution to show how Fourier-transform works, but to get precise result we have to consider technical details and calculate the convolution of a much wider range of sin waves with different frequencies.

The focus of this post is on filtering frequencies and why we need it. And I’m sure you can find dozens of good tutorials on how to take Fourier transform of your data/signals. here is a good source to find out more: https://towardsdatascience.com/fast-fourier-transform-937926e591cb

Filtering:

Having multiple sin waves creating one signal is common when you deal with signals. For instance, while your brain does a task, different regions in your brain working together to achieve that task. Those different areas of your brain are usually oscillating in different frequencies and the result of it creates a complicated signal which you see in the EEG signal.

As a neuroscience student, I have to consider these cooperations and to study only one area of the brain, I may need to filter out other sources. But these technic has a wide application in telecommunication, smoothing data by sub-pressing noises and in general studying any time-series which has periodic sources as well. Talking of noise, as we saw above the convolution of a simple line shape kernel with a signal also smooth the signal and remove the sharp edges of the signal. But convolution is not the best method of reducing the noise when it comes to perfection.

What we did is that we filtered the Fourier transform signal and then reverse it to the original signal. And once we do this process the noises from the main signal vanishes.

There is another kind of filtering in which the only thing remains is a narrow range in the frequency domain and all other frequencies are gone. So that you can focus on the exact oscillation you want.

Writing the Algorithm:

First you need to import these libraries:

import numpy as np
import pandas as pd
from scipy.fftpack import fft,ifft
from scipy.signal import find_peaks,blackman

• numpy and pandas libraries are really handy ones for dealing with arrays.
• fft, ifft modules imported from scipy.fftpack to get a Fast Fourier-transform and also to take a reverse signal from a Fourier-transform of a signal.
• find_peaks and blackman are also needed to first find the peak value of the signal, and second to create a bell-shape curve for the bandpass filtering process.

The next thing we need to do is to create a reverse sigmoid function to create the lowpass filtering.

def sigmoid(range):
    x = np.linspace(-6,+6,range)
    y = 1/(1 + np.exp(-x))
    return -y + 1

here, we reversed the sigmoid function and also shifted it to make a push curve for the lowpass filter.

The next thing would be importing the signal. As I’m using a simulated brain signal, the below part has been specified only for me. But you can import your own signal for your purpose (If you look for data-points I’m using, you can find it in the Github link at the end of this post).

data = pd.read_csv('lfp.csv')
lfp = np.array(data['ex_lfp'])

• I imported the data I’m dealing with from the same folder my python file is located in.
• I extracted the data I want from the pandas Dataframe and stored it into an array named “lfp”. here lfp is my main and raw signal.

N = lfp.size
T = 1.0/1000.0

• I get the length of the data and stored it in N.
• T which is the time resolution of your data-points is a very important constant and based on your data it can be different. I’m working with data-points that recorded every 1 ms so the time resolution of my data-points would be 1/1000 of a second (Notice that if you change this constant you’ll get different frequencies as result).

fourierT = fft(lfp)
fourierT = fourierT[1:N]

• In the first line, we took Fast Fourier-transform of the signal with fft function of scipy.fftpack.
• because of zero frequences in the frequency domain which is the average of the signal in the time domain. That’s why we only pick the [1:N] part of the array. You can find a further description of this link.(https://www.quora.com/What-is-a-zero-Hz-frequency-component)

So to this point, we get the Fourier-transform of the signal. Now it’s time to find the peaks in the frequency domain. To do so, numpy has a function which finds the peak in a very simple way.

peak = np.argmax(fourierT)

• np.argmax is a function to find the index where valued with the maximum number in an array.

After finding the peak in the array, it’s time to make a kernel and a sigmoid curve.

# Lowpass filter process
range = 300
y = sigmoid(range)
kernel = np.ones(len(fourierT))
kernel[peak:peak+len(y)] = y
kernel[peak+len(y):] = 0

• Defining a range of 300, this range can tell you how fast or slow your sigmoid curve will decay.
• We create a kernel consist of ones with the length of the Fourier-transformed signal.
• We start a decreasing sigmoid curve at the peak point and after that, the kernel values must be zero.

If we want to design a bandpass filter, instead of this part we had to follow a different path:

# Bandpass filter
blackman_ = blackman(2*range)
kernel = np.zeros(len(signal_fft))
kernel[fft_peak-range:fft_peak+range] = blackman_

• blackman function which imported from scipy.signal is a bell-shape curve.
• We create a kernel consist of zeros with the length of the Fourier-transform signal.
• We create our bell-shaped curve at the position of the Fourier-transform peak.

The final part of filtering is just simply multiplying each element of the kernel to each element of the Fourier-transform signal.

filtered = np.multiply(kernel,fourierT.real)

• np.multiply is a function that multiplies each element of two given arrays and returns an array as a result.
• because of the imaginary part of the Fourier-transform, we must do this multiplication in the real value domain. To achieve it we simply need to add .real term at the end of the arrays name. numpy made a convenient way of getting the real or imaginary part of an array with this instruction.

So the filtering part has finished. We only need to take an inverse Fourier-transform to get back to the time domain signal instead of the frequency domain.

inverse_signal = ifft(filtered)

• ifft function from scipy.fftpack would create a signal based on what frequencies exist.

The only step left is that get the real values of this transformation.

inverse_signal = inverse_signal.real

Here is the code is written in a clean way of all we discussed above.

And this is the result signals after bandpass and lowpass filtering.

I hope by reading this tutorial you can understand the Fourier-transform a bit more and start to design your own filters from scratch.

To look even further into the codes, especially the convolution part, I invite you to check out my Github page that contains these codes.

As this is my very first post in Medium, I would be happy to read your ideas about this post and the future posts. Peace ;).