# Using Fourier transform for time series decomposition

In this chapter, we learn how to make use of Fast Fourier Transform (FFT) to deconstruct time series.

We start with an easy example. Let’s recap the example from the *Basic time series analysis* chapter, where we’ve generated an artificial signal. We modify it a bit in order to better account for the sampling frequency of our signal:

```
# Create some dummy data:
import numpy as np
import matplotlib.pyplot as plt
# imports most relevant Matplotlib commands
# define frequencies, amplitudes, and sampling rate and time array:
f1 = 2 # Frequency 1 in Hz
f2 = 10 # Frequency 2 in Hz
A1 = 6 # Amplitude 1
A2 = 2 # Amplitude 2
Fs = 100 # Sampling rate
t = np.arange(0,1,1/Fs)
# calculate prime signals:
A_sin = A1 * np.sin(2 * np.pi * f1 * t)
A_cos = A2 * np.cos(2 * np.pi * f2 * t)
A_signal = A_sin + A_cos
# add some noise:
np.random.seed(1)
A_Noise = 2
Noise = np.random.randn(len(t)) * A_Noise
A_signal_noisy = A_signal + Noise
# plots:
fig=plt.figure(3, figsize=(9,4))
plt.clf()
plt.plot(t, A_sin, label="sine", lw=5, alpha=0.7)
plt.plot(t, A_cos, label="cosine", lw=5, alpha=0.7)
plt.plot(t, A_signal, lw=5, c="mediumorchid",
label="superposition", alpha=0.75)
plt.plot(t, A_signal_noisy, lw=5, c="lime",
label="superposition + noise", alpha=0.5)
plt.legend(bbox_to_anchor=(1.0, 1.0), loc='upper left')
plt.xticks([0, 0.25, 0.5, 0.75, 1],
["0", r"$\frac{\pi}{2}$", r"$\pi$",
r"$\frac{3}{4}\pi$", r"$2\pi$"])
plt.tight_layout()
plt.show()
```

**Note**: The expressions within the `xtick`

command of the above plot are $\LaTeX$ commands. You can find further information about the usage of $\LaTeX$ in *Markdown* documents and, hence, also in *Jupyter* notebooks here. An overview of the most common $\LaTeX$ commands can be found in this guide.

We again see, that our signal is a superposition of a sine and cosine wave (`A_signal`

) plus some random noise (`A_signal_noisy`

). Now, we make use of Fourier Transform in order to deconstruct the signal.

## Fourier Analysis

**Fourier analysis**ꜛ, also know as harmonic analysis, is the mathematical field of *Fourier series* and *Fourier integrals*. A **Fourier series**ꜛ decomposes any periodic function (or signal) into the (possibly) infinite sum of a set of simple sine and cosine functions or, equivalently, complex exponentials. The **Fourier transform**ꜛ is a tool for decomposing functions depending on *space* or *time* into functions depending on their component *spatial* or *temporal frequency*.

## Fourier Transform in Python

For Python, where are several Fast Fourier Transform implementations availble. Here, we will use the `fft`

function from the `scipy.fft`

package:

```
import scipy.fft
A_signal_fft = scipy.fft.fft(A_signal)
frequencies = scipy.fft.fftfreq(np.size(t), 1/Fs)
fig=plt.figure(2, figsize=(15,6))
plt.clf()
plt.plot(frequencies, np.abs(A_signal_fft), lw=1.0, c='paleturquoise')
plt.stem(frequencies, np.abs(A_signal_fft))
plt.xlabel("frequency [Hz]")
plt.ylabel("amplitude [a.u.]")
plt.title(r"$|\mathcal{F}(A_{signal})|$")
```

```
frequency_eval_max = 100
A_signal_rfft = scipy.fft.rfft(A_signal, n=frequency_eval_max)
n = np.shape(A_signal_rfft)[0] # np.size(t)
frequencies_rel = n*Fs/frequency_eval_max * np.linspace(0,1,int(n))
fig=plt.figure(3, figsize=(15,6))
plt.clf()
plt.plot(frequencies_rel, np.abs(A_signal_rfft), lw=1.0, c='paleturquoise')
plt.stem(frequencies_rel, np.abs(A_signal_rfft))
plt.xlabel("frequency [Hz]")
plt.ylabel("amplitude [a.u.]")
plt.title(r"$|\mathcal{F}(A_{signal})|$")
```

The Fourier transformed signal $|\mathcal{F}(A_{signal})|$ can also be transformed back to the spatial domain, $|\mathcal{F}^{-1}(A_{signal})|$, by applying the fast inverse Fourier transform function `irfft`

:

```
A_signal_irfft = scipy.fft.irfft(A_signal_rfft)
fig=plt.figure(4, figsize=(9,4))
plt.clf()
plt.plot(t, A_sin, label="sine", lw=5, alpha=0.7)
plt.plot(t, A_cos, label="cosine", lw=5, alpha=0.7)
plt.plot(t, A_signal, lw=6, c="mediumorchid",
label="superposition", alpha=0.75)
plt.plot(t, A_signal_irfft, c='k',
label="$|\mathcal{F}^{-1}(A_{signal})|$")
plt.legend(bbox_to_anchor=(1.0, 1.0), loc='upper left')
plt.xticks([0, 0.25, 0.5, 0.75, 1],
["0", r"$\frac{\pi}{2}$", r"$\pi$",
r"$\frac{3}{4}\pi$", r"$2\pi$"])
plt.tight_layout()
plt.show()
```

### Exercise 1

- Apply and plot the Fast Fourier Transform to the noisy signal
`A_signal_noisy`

. - What do you notice?

```
# You solution 1 here:
```

**Toggle solution**

## Exercise 2: Frequency filtering by using Fourier Transform frequencies

- Use the function
`find_closest_within_array`

defined below to find the index`idx`

of the value within the`frequencies_rel`

array, that is closest to the frequency of 10 Hz. - Use the output index
`idx`

to set the corresponding amplutide of the Fourier transformed signal`A_signal_rfft`

to zero. Then, retransform the Fourier signal and plot the filtered signal. - Reapply your script for a filter frequency of 2 Hz.
- Reapply your script for the noisy signal.

```
def find_closest_within_array(array, value):
""" Finds closest value within an array.
array: input NumPy array (can also be a list)
value: value to search for
array[idx]: found closest value within array
idx: index of the found closest value
"""
numpy_array = np.asarray(array)
idx = (np.abs(numpy_array-value)).argmin()
return array[idx], idx
```

```
# Your solution 2 here:
```

**Toggle solution**

## Exercise 3: Frequency filtering by using Fourier Transform amplitudes

With the filter approach shown above we are able to specifically filter for a certain frequency. However, to filter, e.g., random noise that is present in almost every frequency (with low amplitudes though), this approach will not work. For such cases we need an alternative approach, that filters for amplitudes within the Fourier transformed signal instead for the corresponding frequency:

- Define a threshold
`A_pass_limit = 50`

, that sets the lower bound of the amplitudes we want to keep within the Fourier transform $|\mathcal{F}(A_{signal})|$ of our noisy signal (`A_signal_rfft`

). Filter`A_signal_rfft`

by setting all of its values, that are lower than`A_pass_limit`

, to zero. - Retransform the Fourier signal and plot the resulting filtered signal.

```
# Your solution 3 here:
```

**Toggle solution**

</script>
**Note**: ‘$|\mathcal{F}^{-1}(A_{signal, filtered})|$’ is $\LaTeX$.

## Exercise 4: Finalize your pipeline

- Put the filter solution from Exercise 2 into a function, e.g., called
`frequency_filter`

. Define the function in such a way, that it outputs the filtered signal, the Fourier transform of the signal, the filtered Fourier transform of the signal as well as the corresponding frequency array (`frequencies_rel`

). - Apply your function to
`A_signal`

as well as to`A_signal_noisy`

, both for a filter frequency of 10 and 2 Hz, respectively. - Repeat 1. and 2. for the solution from Exercise 3 (call the function, e.g.,
`amplitude_filter`

).

**Hint**: It is advisable to also put your plot commands into a plot function in order to avoid code repetitions within your script.

```
# Your solution 4.1 and 4.2 here:
"""Hints:
def frequency_filter(signal, filter_frequency,
frequency_eval_max):
...
return (e.g.) signal_filtered, signal_rfft_filtered,
signal_rfft, frequencies_rel
def plot_comparison(frequencies_rel, Current_signal_rfft,
Current_signal, fignum=1):
...
return NOTHING
"""
```

```
# Your solution 4.1 and 4.2 here:
```

**Toggle solution**

```
# Your solution 4.3 here:
```

**Toggle solution**

## Exercise 5: Compbine both filter approaches

So far, we have applied both filter functions separately. In order to, e.g., apply noise reduction and the exclusion of a certain frequency, you can combine both filters by applying them one after the other:

- Apply the
`amplitude_filter`

function to the noisy signal`A_noisy_signal`

in order to filter out the noise. - Re-use the resulting filtered signal and apply the
`frequency_filter`

function to it (filter for the frequency 2 or 10 Hz).

```
# Your solution 5 here:
```

**Toggle solution**

## Exercise 6: Application to real world data

Now, we apply our two filter functions to the data from the *Patch Clamp Analysis* chapter:

- Recap the
*Analyzing patch clamp recordings*chapter and write a script, that reads the IGOR file`ad1_12.ibw`

(lies within the`/Data1`

folder). - Apply a Fast Fourier Transform in order to assess the prominent frequencies and potential noise within the recording. What is the frequency of the most dominant signal within the recording?
- Apply your filter functions to filter out any available noise and any unwanted frequency.
- Again, recap the
*Analyzing patch clamp recordings*chapter and add a section to your script, in which you assess the prominent peaks within the recording by applying the`find_peaks`

function from the`scipy`

package. From these results, estimate the frequency of the dominant signal within the recording. Discuss differences (if any) between the detected peak frequency from this method and from the Fourier Analysis in 2.

```
# Your solution 6 here:
```