```python
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
```

FT algorithm receives a trajectory, apply its filters to find the appropriate cycles, and outputs the full set of cyclic components. There are two algorithms:

- the Discrete Fourier Transform (DFT) which requires $O(n^2)$ operations (for n samples)
- the Fast Fourier Transform (FFT) which requires $O(nlog(n))$ operations

## DFT
\begin{equation}
X_k = \sum_{n=0}^{N-1} x_n e^{-i 2 \pi k n / N}
\end{equation}

\begin{equation}
x_n = \frac{1}{N} \sum_{k=0}^{N-1} X_k e^{i 2 \pi k n / N}
\end{equation}

where:
- $X_k$ amount of frequency $k$ in the signal; each $k$th value is a complex number including strength (amplitute) and phase shift
- $N$ number of samples
- $n$ current sample, $n\in{0\cdots N−1}$
- $k$ current frequency, between $0$ Hz to $N-1$ Hz
- $1/N$ not necessary but it gives the actual sizes of the time spikes
- $n/N$ is the percent of the time we’ve gone through
- $2\pi{k}$ the speed in radians/second
- $e^{−ix}$ the backwards-moving circular path. This last three tell how far we’ve moved, for this speed and time.


```python
# step by step dft
def dft_k(x, k):
    N = len(x)
    return sum((x[n]*np.e**(-1j*2*np.pi*k*n/N) for n in range(N)))

# step by step idft
def idft_n(X, n):
    N = len(x)
    return sum((1/N * X[k] * np.e**(1j*2*np.pi*k*n/N) for k in range(N)))
```


```python
def remove_frequencies(Xn, N, fs, fre_low, fre_high):
    # remove specific frequencies
    Yn = np.copy(Xn)
    fre_low_n = fre_low * N // fs
    fre_high_n = fre_high * N // fs
    # remove two side frequiencies
    Yn[fre_low_n:fre_high_n] = 0
    Yn[N-fre_high_n:N-fre_low_n] = 0
    return Yn


N = 1000
fs = 10000
T = 1/fs
t = np.linspace(0, N * T, N)
x = np.sin(2*np.pi*50*t) + 2 * np.sin(2*np.pi*150*t) + 0.5 * np.sin(2*np.pi*1000*t)

# dft: analysis
Xn = np.zeros(N, dtype=np.complex)
for i in range(N):
    X = dft_k(x, i)
    Xn[i] = X

# filter: remove specific requencies
Yn = remove_frequencies(Xn, N, fs, 800, 1200)
    
# idft: synthesis
Re = np.zeros(N, dtype=np.complex)
for i in range(N):
    Re[i] = idft_n(Yn, i)

# one side frequency range
xf = np.linspace(0.0, 1.0/(2.0*T), N//2)
xn = 2.0/N * np.abs(Xn[0:N//2])
yn = 2.0/N * np.abs(Yn[0:N//2])
```


```python
#------------------------------------------------------------
# Set up the plots
fig = plt.figure(figsize=(10, 10))
fig.subplots_adjust(left=0.09, bottom=0.09, right=0.95, top=0.95,
                    hspace=0.05, wspace=0.05)
#----------------------------------------
# plot the origional signal
ax1 = fig.add_subplot(221)
ax1.grid(color='#b7b7b7', linestyle='-', linewidth=0.5, alpha=0.5)
ax1.plot(t, x, '-k', label=r'data $D(x)$')
plt.setp(ax1.get_xticklabels(), visible=False)

#----------------------------------------
# plot the dft and area to remove
ax2 = fig.add_subplot(222)
ax2.plot(xf, xn, '-k')
ax2.grid(color='#b7b7b7', linestyle='-', linewidth=0.5, alpha=0.5)
ax2.axvspan(800, 1200, facecolor='#b7b7b7', alpha=0.5)
ax2.yaxis.tick_right()
plt.setp(ax2.get_xticklabels(), visible=False)

#----------------------------------------
# plot the left frequencies
ax3 = fig.add_subplot(224, sharex=ax2)
ax3.plot(xf, yn, '-k')
ax3.grid(color='#b7b7b7', linestyle='-', linewidth=0.5, alpha=0.5)
ax3.yaxis.tick_right()

#----------------------------------------
# plot the filtered signal
ax4 = fig.add_subplot(223, sharex=ax1)
ax4.plot(t, Re.real, '-k')
ax4.grid(color='#b7b7b7', linestyle='-', linewidth=0.5, alpha=0.5)

#------------------------------------------------------------
# Plot flow arrows
ax = fig.add_axes([0, 0, 1, 1], xticks=[], yticks=[], frameon=False)

arrowprops = dict(arrowstyle="simple",
                  color="#333333", alpha=0.5,
                  shrinkA=5, shrinkB=5,
                  patchA=None,
                  patchB=None,
                  connectionstyle="arc3,rad=-0.35")

ax.annotate('', [0.57, 0.57], [0.47, 0.57],
            arrowprops=arrowprops,
            transform=ax.transAxes)
ax.annotate('', [0.57, 0.47], [0.57, 0.57],
            arrowprops=arrowprops,
            transform=ax.transAxes)
ax.annotate('', [0.47, 0.47], [0.57, 0.47],
            arrowprops=arrowprops,
            transform=ax.transAxes)
plt.show()
```

## Reference:

- [Fourier Transform: A R Tutorial](http://www.di.fc.ul.pt/~jpn/r/fourier/fourier.html)
- [An Interactive Guide To The Fourier Transform](https://betterexplained.com/articles/an-interactive-guide-to-the-fourier-transform/)
- [An Interactive Introduction to Fourier Transforms](http://www.jezzamon.com/fourier/index.html)

2383e60b9e60024416a27166a4bb97d40ee08f1a
ipynb
Jupyter Notebook
003_fourier_transform.ipynb
ValleyZw/hill
51d4bc76d7c49b3dccf88659fade806906a9b4e2
MIT
2019-07-22T07:19:39.000Z
2019-07-22T07:19:39.000Z
003_fourier_transform.ipynb
ValleyZw/hill
51d4bc76d7c49b3dccf88659fade806906a9b4e2
MIT
003_fourier_transform.ipynb
ValleyZw/hill
51d4bc76d7c49b3dccf88659fade806906a9b4e2
MIT
2020-05-19T17:52:16.000Z
2020-05-19T17:52:16.000Z
Qwen/Qwen-72B
1. YES
2. YES