###### tags: `Lesson 2`
# 13. Fourier Transform In Practice
{%youtube 3IWA8tMxxTA %}
We learned how to apply the theory and understanding of what a fourier transform is to actually solve some problems.
Let's make a signal that is composed of two sine waves at 2 Hz and 3 Hz plus some random noise. We will be using the `np.fft` module to compute the Fourier transform with the two main functions below:
* `rfft` - computes the actual Fourier transform coefficients
* `rfftfreq` - tells us the frequencies for which we are computing the Fourier transform
We then examined `freqs` to see that the FFT samples the Fourier transform uniformly from 0 Hz to the Nyquist frequency, which in this case is 25 Hz because our sampling rate is 50 Hz. We also saw the Fourier transform coefficients and only examined the magnitudes. Plotting the FFT, we see that the signal is composed primarily of two frequencies (2 and 3Hz) and a little bit of everything else (from the random noise).
We also saw how zero-padding could be used to visualize sinusoids with frequencies that are not present in `freqs`. To visualize at frequencies not in `freqs`, we need to sample twice as often and we can do this by adding 0s to the end of the signal. However, we also see this rippling, which is an artifact of padding the signal with 0s, which is the trade-off when doing zero-padding.
{%youtube W-5uoeNYFwM %}
We just learned about inverse Fourier Transform in which we used the `np.fft` module to compute the inverse Fourier transform with the function `irfft`.
We started with a noisy signal and removed all the frequency components not in the range, in this case, 2Hz and 3Hz. And we saw a recovered signal that looked very close to the signal we saw in the previous video.
But we did the process again from 2.15Hz and 2.95Hz. The recovered signal looked a bit distorted and not what we'd expect. This is because zeroing out Fourier coefficients is not the best way to filter a signal.
We then used `scipy` to **bandpass filter** our signal for us. A bandpass filter will remove all frequency components outside of a given passband. Let's bandpass filter our signal with a passband from 1 Hz to 4 Hz. This way, our desired frequencies of 2.15 Hz and 2.95 Hz are well within the passband. And now, our recovered signal looks very similar to what we want.
In this course, we will always bandpass our signals before processing them.
## Notebook Review
If you wanted to interact with the notebook in the video, you can access it [here](https://github.com/udacity/nd320-c4-wearable-data-starter/tree/master/intro-to-dsp/walkthroughs/fourier-tranform-II) in the repo `/intro-to-dsp/walkthroughs/fourier-transform-II/` or in the workspace below.
# Fourier Transform In Practice
Now that we know what the Fourier transform is, let's use it to solve some problems.
## Imports
Begin with our traditional imports.
```python
import numpy as np
from matplotlib import pyplot as plt
import mpld3
import scipy as sp
from scipy import io, signal
```
```python
%matplotlib inline
```
```python
import mpld3
mpld3.enable_notebook()
```
## Computing the Fourier Transform
Let's make a signal that is composed of two sine waves at 2 Hz and 3 Hz plus some random noise.
```python
fs = 50
ts = np.arange(0, 10, 1/fs)
s2 = np.sin(2 * np.pi * 2 * ts)
s3 = np.sin(2 * np.pi * 3 * ts)
s = s2 + s3 + np.random.randn(len(s2)) * 1
```
We use the `np.fft` module to help us compute the Fourier transform. `rfftfreq` tells us the frequencies for which we are computing the Fourier transform. `rfft` computes the actual Fourier transform coefficients. ==We use the `rfft` and `rfftfreq` functions instead of `fft` and `fftfreq` because we know our signals will only ever contain real numbers (as opposed to complex). You don't really need to know the difference for now.==
```python=
freqs = np.fft.rfftfreq(len(s), 1/fs)
fft = np.fft.rfft(s)
```
By examining `freqs`, we see that the FFT samples the Fourier transform uniformly from 0 Hz to the Nyquist frequency, which in this case is 25 Hz because our sampling rate is 50 Hz.
```python=
freqs
```
The actual Fourier transform coefficients are complex numbers because they describe not only the magnitude of the frequency component in the signal but also its phase. Most of the time, we're just interested in the magnitude, which we can get by using `np.abs`.
```python=
np.abs(fft)
```
array([ 57.23573065, 15.81111803, 4.04847168, 13.55639643,
22.64158495, 16.88440975, 8.11815656, 7.09311266,
17.31272093, 26.54198613, 13.03583761, 34.63634805,
3.25479059, 18.26256333, 24.66610087, 21.4483017 ,
9.50626426, 13.29291602, 41.79481683, 2.22190195,
221.81645738, 19.38819958, 44.15817548, 19.12564097,
32.66065773, 15.75070478, 12.2414837 , 27.87930597,
11.96215353, 24.36754827, 264.47441665, 12.70339451,
14.06123698, 1.6575093 , 24.26179513, 16.43748261,
22.97629983, 15.77443929, 1.88977759, 12.31987813,
4.10592962, 18.40480778, 20.3052585 , 7.8117142 ,
14.34117185, 8.44623743, 22.04743678, 25.48385078,
19.77189159, 21.02499256, 44.78890599, 19.8593788 ,
11.89292656, 15.91103787, 20.06163464, 9.49583923,
11.01839388, 24.72530113, 9.97948017, 2.41246779,
12.35690306, 43.50680104, 25.18858406, 25.65354494,
16.060977 , 36.3357762 , 57.26107685, 20.17138541,
41.97401414, 12.48764345, 19.56505247, 31.00880526,
21.47574229, 28.45953768, 8.47023415, 15.60446362,
5.40720479, 24.11965975, 22.52138059, 7.96563089,
5.79007533, 19.60162759, 7.18729822, 22.87463294,
20.5456545 , 28.93256436, 25.34238427, 25.00378292,
28.33844491, 35.81759994, 6.84924553, 3.99814796,
25.02090768, 6.21397642, 28.23073807, 23.29993058,
10.04764903, 16.69528057, 7.77220237, 8.15291002,
14.66699077, 21.88878539, 17.02476864, 19.52555584,
40.1169293 , 26.1242745 , 13.67944715, 18.48474865,
7.83988048, 26.06706412, 14.74317867, 21.11863107,
12.76041534, 16.16019246, 26.13792391, 10.35888236,
24.00071809, 27.49706053, 29.69543304, 2.53484552,
29.70441008, 31.50030439, 31.41777508, 36.04420734,
8.61069431, 3.90358183, 21.4456407 , 31.21340046,
17.38927261, 22.43331028, 34.2405972 , 12.22322696,
37.48714987, 6.68763065, 27.75094964, 18.28698861,
13.07879412, 6.26133067, 10.75650603, 29.78926695,
21.83782506, 19.62477635, 20.24230622, 19.29787633,
25.9687386 , 14.44015085, 17.51506284, 12.85151354,
27.43994841, 29.11267525, 32.5448301 , 6.68983444,
15.30859223, 19.80929442, 10.28165278, 16.26463274,
27.23387184, 35.22829508, 21.05117788, 17.186299 ,
12.71696473, 24.2914819 , 22.35454114, 24.12299966,
11.87196402, 29.47285742, 32.73585953, 23.16364369,
36.84252941, 28.87909979, 21.2028534 , 14.32646292,
12.63785451, 12.84729237, 8.85980606, 4.51603236,
26.77216015, 11.03725662, 18.84952343, 32.74223009,
26.69170059, 12.31486851, 20.74128832, 0.58772494,
19.18160917, 38.56186877, 25.29632365, 26.08965958,
9.76739424, 14.42766212, 30.23134488, 30.57138499,
26.30623779, 29.22991002, 27.1639673 , 40.21350608,
42.41884914, 4.36775801, 20.87095431, 17.40952762,
18.34151815, 28.42425873, 41.34668416, 36.82538021,
22.33787765, 21.23447509, 12.49464797, 4.06532782,
22.75834777, 36.26159089, 1.88132571, 18.14879415,
10.17472522, 15.94825425, 24.04556175, 11.17120224,
9.98513904, 23.6754109 , 20.35298906, 5.94729494,
22.58755169, 16.1293238 , 17.61756664, 17.89113194,
14.37757889, 47.55226292, 23.9687635 , 27.99620005,
30.11238718, 9.69166712, 28.38911805, 34.48080686,
13.63620987, 13.14677568, 19.28828614, 9.62600025,
23.37538339, 20.92706292, 9.87772897, 20.16128758,
6.22405061, 15.71870427, 28.40314764, 5.49247389,
18.95152861, 37.55525939, 5.38105759, 17.55971456,
15.78066845, 24.36638947, 18.03367689])
The output of the real FFT is half the length of the input signal (plus 1).
```python=
len(fft), len(freqs), len(s)
```
(251, 251, 500)
Plotting the FFT, we see that the signal is composed primarily of two frequencies (2 and 3Hz) and a little bit of everything else (from the random noise).
```python=
plt.figure(figsize=(12, 8))
plt.subplot(2,1,1)
plt.plot(ts, s)
plt.title('Time-Domain')
plt.xlabel('Time (sec)')
plt.subplot(2,1,2)
plt.plot(freqs, np.abs(fft))
plt.title('Frequency-Domain')
plt.xlabel('Frequency (Hz)')
plt.tight_layout()
```
## Zero Padding
Let's change the frequency slightly to be in between the `freqs` array above.
```python=
fs = 50
ts = np.arange(0, 10, 1/fs)
s2 = np.sin(2 * np.pi * 2.35 * ts)
s3 = np.sin(2 * np.pi * 3.65 * ts)
s = s2 + s3
```
The peaks in the Fourier transform aren't as sharp and there's no value for the true frequency 2.35 Hz. So with just this plot, it's hard to tell what the true frequency of the signal is. We're not sampling the Fourier transform with high enough resolution to see the power at 2.35 Hz.
```python=
freqs = np.fft.rfftfreq(len(s), 1/fs)
fft = np.fft.rfft(s)
plt.figure(figsize=(12, 8))
plt.subplot(2,1,1)
plt.plot(ts, s)
plt.title('Time-Domain')
plt.xlabel('Time (sec)')
plt.subplot(2,1,2)
plt.plot(freqs, np.abs(fft))
plt.title('Frequency-Domain')
plt.xlabel('Frequency (Hz)')
plt.tight_layout()
```
Remember, the real Fourier transform's output is half the length of the signal input and that it will always evenly sample all frequencies from 0 to the Nyquist frequency. So at a sampling rate of 50 Hz and an input signal of length 500, it will give us 251 samples between 0 and 25. This means we're sampling at 0 and 0.1 and 0.2 Hz as we saw in the freqs array, but we don't sample at 2.35 Hz, which is the actual frequency of the signal. To do that, we actually need to sample twice as often. We need 500 samples between 0 and 25 Hz. We can do this by adding 0s to the end of the signal.
```python=
padded_s = np.hstack((s, np.zeros(len(s))))
padded_ts = np.arange(0, 20, 1/fs)
freqs = np.fft.rfftfreq(len(padded_s), 1/fs)
fft = np.fft.rfft(padded_s)
plt.figure(figsize=(12, 8))
plt.subplot(2,1,1)
plt.plot(padded_ts, padded_s)
plt.title('Time-Domain')
plt.xlabel('Time (sec)')
plt.subplot(2,1,2)
plt.plot(freqs, np.abs(fft))
plt.xlim((1, 5))
plt.title('Frequency-Domain')
plt.xlabel('Frequency (Hz)')
plt.tight_layout()
```
Now we have all the points we need and see a nice sharp peak at 2.35 Hz and 3.65 Hz. However, we also see this rippling, which is an artifact of padding the signal with 0s, which is the trade-off when doing 0-padding.
Numpy makes this convenient for us and has a parameter to the function that lets us take a 0-padded FFT without actually having to pad the signal ourselves.
```python=
fftlen = len(s) * 2
freqs = np.fft.rfftfreq(fftlen, 1/fs)
fft = np.fft.rfft(s, fftlen)
```
```python=
plt.figure(figsize=(12, 8))
plt.subplot(2,1,1)
plt.plot(ts, s)
plt.title('Time-Domain')
plt.xlabel('Time (sec)')
plt.subplot(2,1,2)
plt.plot(freqs, np.abs(fft))
plt.xlim((1, 5))
plt.title('Frequency-Domain')
plt.xlabel('Frequency (Hz)')
plt.tight_layout()
```
## Inverse FFT
The Fourier transform is invertible. This means we can modify the signal in "frequency-domain" and then take the inverse FFT to recover our "time-domain" signal.
Recreate the original noisy signal.
```python=
fs = 50
ts = np.arange(0, 10, 1/fs)
s2 = np.sin(2 * np.pi * 2 * ts)
s3 = np.sin(2 * np.pi * 3 * ts)
s = s2 + s3 + np.random.randn(len(s2)) * 1
freqs = np.fft.rfftfreq(len(s), 1/fs)
fft = np.fft.rfft(s)
```
Let's remove all frequency components that are not at 2Hz and 3Hz.
```python=
new_fft = fft.copy()
new_fft[~np.isin(freqs, (2, 3))] = 0
```
```python
plt.figure(figsize=(12, 8))
plt.plot(freqs, np.abs(new_fft))
plt.title('Frequency-Domain')
plt.xlabel('Frequency (Hz)')
```
Text(0.5, 0, 'Frequency (Hz)')
We can get our signal back in the time-domain by calling `irfft` on the Fourier transform. After plotting the new signal, we can see that by removing all the other frequency components, we have significantly cleaned up our signal! Imagine doing this in the time-domain.
```python
%matplotlib inline
```
```python
new_s = np.fft.irfft(new_fft)
```
```python
plt.clf()
plt.subplot(3, 1, 1)
plt.plot(ts, s2 + s3, label='2Hz + 3Hz')
plt.ylim((-4, 4))
plt.legend()
plt.subplot(3, 1, 2)
plt.plot(ts, s, label='input signal')
plt.ylim((-4, 4))
plt.legend()
plt.subplot(3, 1, 3)
plt.plot(ts, new_s, label='output signal')
plt.ylim((-4, 4))
plt.legend()
plt.tight_layout()
```
Let's change the frequencies again slightly. Now the input signal is 2.15 Hz + 2.95 Hz + random noise
```python
s2 = np.sin(2 * np.pi * 2.15 * ts)
s3 = np.sin(2 * np.pi * 2.95 * ts)
sb = s2 + s3 + np.random.randn(len(s2)) * 1
fft = np.fft.rfft(sb)
new_fft = fft.copy()
new_fft[~np.isin(freqs, (2.1, 2.2, 2.9, 3.0))] = 0
new_sb = np.fft.irfft(new_fft)
plt.clf()
plt.subplot(3, 1, 1)
plt.plot(ts, s2 + s3, label='2.15Hz + 2.95Hz')
plt.ylim((-4, 4))
plt.legend()
plt.subplot(3, 1, 2)
plt.plot(ts, sb, label='input signal')
plt.ylim((-4, 4))
plt.legend()
plt.subplot(3, 1, 3)
plt.plot(ts, new_sb, label='output signal')
plt.ylim((-4, 4))
plt.legend()
plt.tight_layout()
```
The recovered signal looks a bit distorted and not what we'd expect. This is because zeroing out Fourier coefficients is not the best way to filter a signal.
## Bandpass Filter
We can use `scipy` to **bandpass filter** our signal for us. A bandpass filter will remove all frequency components outside of a given passband.
```python
def BandpassFilter(signal, pass_band):
b, a = sp.signal.butter(5, pass_band, btype='bandpass', fs=50)
return sp.signal.filtfilt(b, a, signal)
```
Let's bandpass filter our signal with a passband from 1 Hz to 4 Hz. This way, our desired frequencies of 2.15 Hz and 2.95 Hz are well within the passband.
```python
filtered_s = BandpassFilter(sb, (1, 4))
```
And now, our output signal looks very similar to what we want.
```python=
plt.clf()
plt.subplot(3, 1, 1)
plt.plot(ts, s2 + s3, label='2.15Hz + 2.95Hz')
plt.ylim((-4, 4))
plt.legend()
plt.subplot(3, 1, 2)
plt.plot(ts, sb, label='input signal')
plt.ylim((-4, 4))
plt.legend()
plt.subplot(3, 1, 3)
plt.plot(ts, filtered_s, label='output signal')
plt.ylim((-4, 4))
plt.legend()
plt.tight_layout()
```
In this course, we will always bandpass our signals before processing them.
Sign in with Wallet
Connect another wallet