Довідничок для перетворення Фур'є

Період 2$\pi$

Фур'є

Припустімо, що $f(t)$ має період $2\pi$, тоді $$\begin{equation} f(t) = \frac{a_0}{2} + \sum_{n=1}^{\infty}\left(a_n\cos(nt) + b_n \sin(nt)\right)\label{eq:FourierPi} \end{equation}$$ де коефіцієнти $a_0, a_1, \dots, b_1, \dots$ обчислюються так: $$\begin{align} a_n &= \frac{1}{\pi}\int_{-\pi}^{\pi}f(t)\cos(nt)dt\label{eq:anpi}\\ b_n &= \frac{1}{\pi}\int_{-\pi}^{\pi}f(t)\sin(nt)dt\label{eq:bnpi}\\ \end{align}$$ Знак рівності може не виконуватись в місцях розривів функції.

Означення

В $\eqref{eq:FourierPi}$ частота $n=1$ називається фундаментальною, частоти з $n > 1$ називаються гармоніками.

Ортогональність базису

$$\int_{-\pi}^{\pi}\cos(nt)\cos(nt)dt = \frac{1}{2} \int_{-\pi}^{\pi}\cos^2(nt)+\sin^2(nt)dt = \frac{1}{2}n2\frac{\pi}{n}=\pi$$ Через це $\pi$ ми і ділимо в $\eqref{eq:anpi}$ і $\eqref{eq:bnpi}$ на $\pi$.

Тепер розглянемо з $n\neq m$: $$\begin{align} \int_{-\pi}^{\pi}\cos(nt)\cos(mt)dt &= \frac{1}{2}\int_{-\pi}^{\pi}\left(\cos((n+m)t) + \cos((n-m)t)\right)dt &= 0\\ \int_{-\pi}^{\pi}\sin(nt)\cos(mt)dt &= \frac{1}{2}\int_{-\pi}^{\pi}\left(\sin((n+m)t) + \sin((n-m)t)\right)dt &= 0\\ \end{align}$$

Період 2L

Фур'є

Тепер припустимо, що $f(t)$ має період $P = 2L$. Тоді ряд Фур'є набуває вигляду $$\begin{equation} f(t) = \frac{a_0}{2} + \sum_{n=1}^{\infty}\left(a_n\cos(n\frac{\pi}{L}t) + b_n\sin(n\frac{\pi}{L}t)\right). \end{equation}$$ А коефіцієнти такого: $$\begin{align} a_n &= \frac{L}{\pi}\int_{-L}^{L}f(t)\cos(n\frac{\pi}{L}t)dt\\ b_n &= \frac{L}{\pi}\int_{-L}^{L}f(t)\sin(n\frac{\pi}{L}t)dt\\ \end{align}$$

Дискретне перетворення Фур'є (ДПФ)

Базис в $\mathbb{C}^N$

Розглянемо вектори $\mathbf{w}^{(k)} = [\dots, e^{j\frac{2\pi}{N}ik}, \dots]'$. Фундаментальною частотою тут буде $\frac{2\pi}{N}k$. Можна побачити, що складові вектора $\mathbf{w}^{(1)}$ роблять один оборот навколо одиничного кола. Найшвидше ми рухаємось по колу коли $k = N/2$ далі, зі збільшенням $k$ швидкість падає, а ознакою того що фундаметальна частота перетнула $\pi$ буде протилежний знак уявної частини.

Базис в $\mathbb{C}^N$

$N$ сигналів $\mathbf{w}^{(k)}_n = e^{j\frac{2\pi}{N}nk}$, де $n,k \in [0, N)$ - це ортогональний базис в $\mathbb{C}^N$.

Для $k\ne h$ $$\begin{align*} \langle\mathbf{w}^{(k)}\mathbf{w}^{(h)}\rangle &= \sum_{n=0}^{N-1}e^{-j\frac{2\pi}{N}nk}e^{j\frac{2\pi}{N}nh}\\ &= \begin{cases} N,& k = h\\ \frac{1 - e^{j\frac{2\pi}{N}N(h-k)}}{1 - e^{j\frac{2\pi}{N}(h-k)}} = 0,& k \ne h\\ \end{cases} \end{align*}$$ Тут ми скористались тим, що $(1+x+\cdots+x^{N-1})(1-x) = 1-x^N$. А також фактом, що $e^{j2\pi(h-k)} = 1$ в результаті чого числівник стає нулем.

Отже, $\mathbf{w}^{(k)}$ формують ортогональний базис для $\mathbb{C}^N$. Але цей базис не ортонормальний. Коефіцієнтом ортонормалізації буде $\frac{1}{\sqrt{N}}$.

Частоти базису

Якщо розглянути як складові базисних векторів рухаються уздовж одиничного кола, то помітно, що спочатку зі сростанням індексу базисного вектора маркер, що позначає складову, рухається все швидше, досягаючі найбільшої швидкості коли $k = N/2 - 1$ і $k = N/2$ для парного $N$. Після цього ми спостерігаємо ефект колеса з фільмів, тобто коли колесо крутиться в зворотній бік, це пов'язано з частотою кадрів. Так і тут, коли маркер робить більше ніж півоберта за крок, то нам здається, що він обертається в зворотньому напрямку з меншою швидкістю.

Python-код генерації анімації.

import numpy as np
import matplotlib.pylab as plt
from matplotlib.animation import FuncAnimation

N = 32
t1 = np.linspace(0,N-1,N)
t10 = np.linspace(0,N-1,N*10)

fig, (ax1, ax2) = plt.subplots(2,1)
fig.set_tight_layout(True)

ax1.set_xlim(-1,N)
ax2.set_xlim(-1,N)

def update(k):
    ax1.cla()
    ax2.cla()
    ax1.set_ylabel("Re")
    ax2.set_ylabel("Im")
    if (k < N):
        if k < N/2:
            ax1.plot(t10, np.cos(2*np.pi/N*(k)*t10))
            ax2.plot(t10, np.sin(2*np.pi/N*(k)*t10))
        else:
            ax1.plot(t10, np.cos(2*np.pi/N*(N - k)*t10))
            ax2.plot(t10, -np.sin(2*np.pi/N*(N - k)*t10))
        ax1.set_title('Базисний вектор #' + str(k) + ' в С^' + str(N))
        markerline, stemlines, baseline = ax1.stem(t1, np.cos(2*np.pi/N*k*t1))
        plt.setp(baseline, color='black', linewidth=2)
        plt.setp(markerline, color='orange')
        plt.setp(stemlines, color='orange', linewidth=2)
        markerline, stemlines, baseline = ax2.stem(t1, np.sin(2*np.pi/N*k*t1), 'violet')
        plt.setp(baseline, color='black', linewidth=2)
        plt.setp(markerline, color='orange')
        plt.setp(stemlines, color='orange', linewidth=2)

anim = FuncAnimation(fig, update, frames=np.arange(0, N-1 + 5), interval=1000)
plt.show()

Зауважте, що уявна частина змінила знак в другій половині. Бо частота стала від'ємною.

До коефіцієнтів і назад

$$\begin{equation} X_k = \langle\mathbf{w}^{(k)},\mathbf{x}\rangle\label{eq:analysis}. \end{equation}$$ $$\begin{equation} \mathbf{x} = \frac{1}{N}\sum_{k=0}^{N-1}X_k\mathbf{w}^{(k)}\label{eq:synthesis}. \end{equation}$$

Зміна базису в матричній формі

Ми можемо записати зміну базису в матричній формі. Якщо позначити $W=e^{-j\frac{2\pi}{N}}$, то $$\begin{equation} \mathbf{W} = \begin{bmatrix} 1 & 1 & 1 & \dots & 1\\ 1 & W & W^2 & \dots & W^{N-1}\\ 1 & W^2 & W^4 & \dots & W^{2(N-1)}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 1 & W^{N-1} & W^{2(N-1)} & \dots & W^{(N-1)^2}\\ \end{bmatrix} \end{equation}$$ Фактично, в кожному рядку маємо спряжений до базового вектора. Тепер формули $\eqref{eq:analysis}$ і $\eqref{eq:synthesis}$ можна записати в матричній формі: $$\begin{equation} \mathbf{X} = \mathbf{W}\mathbf{x}. \end{equation}$$ $$\begin{equation} \mathbf{x} = \frac{1}{N}\mathbf{W}^H\mathbf{X}. \end{equation}$$

Перетворення Фур'є дискретного часу (ПФДЧ)

Якщо збільшувати кількість точок зростає, то фундаментальна частота $2\pi/N$ стає дуже малою і, чисто інтуїтивно, ми можемо казати, що $\lim_{N\to\infty}(2\pi/N)k = \omega$.

Формально

• $x[n] \in l_2(\mathbb{Z})$ (сумовна з квадратом послідовність)
• означимо функцію $$\begin{equation} F(\omega) = \sum_{n=-\infty}^{\infty}x[n]e^{-j\omega n}\label{eq:dtft} \end{equation}$$
• обернена (коли $F(\omega)$ існує): $$\begin{equation} x[n] = \frac{1}{2\pi}\int_{-\pi}^{\pi}F(\omega)e^{j\omega n}d\omega,\quad n\in\mathbb{Z}\label{eq:idtft} \end{equation}$$

Зазвичай замість $F(\omega)$ пишуть $X(e^{j\omega})$. З поміж іншого, при такому записі краще видно $2\pi$ періодичність. Прийнято представляти $X(e^{j\omega})$ на $[-\pi, \pi]$.

Властивості

• лінійність: $$\begin{equation}\mbox{ПФДЧ}\{\alpha x[n] + \beta y[n]\} = \alpha X(e^{j\omega})+\beta Y(e^{j\omega})\end{equation}$$
• часовий зсув: $$\begin{equation}\mbox{ПФДЧ}\{x[n-M]\} = e^{-j\omega M}X(e^{j\omega})\label{eq:dtft_timeshift}\end{equation}$$
• модуляція: $$\begin{equation}\mbox{ПФДЧ}\{e^{j\omega_0 n} x[n]\} = X(e^{j(\omega - \omega_0)})\label{eq:dtft_modulation}\end{equation}$$
• зворотній час: $$\begin{equation}\mbox{ПФДЧ}\{x[-n]\} = X(e^{-j\omega})\end{equation}$$
• спряження: $$\begin{equation}\mbox{ПФДЧ}\{x^*[n]\} = X^*(e^{-j\omega})\end{equation}$$

Властивості $\eqref{eq:dtft_timeshift}$ і $\eqref{eq:dtft_modulation}$ двоїсті одна до одної.

ПФДЧ одиниці

Очевидно, що стала $1$ несумовна з квадратом. Але $\mbox{ДПФ}\{1\} = N\delta[n]$. Було б добре могти обчислити її ПФДЧ. Для цього нам знадобиться дельта-функція Дірака. В ній нас цікавить така властивість: $$\begin{equation} (\delta; f) = \int_{-\infty}^{+\infty}\delta(x-a)f(x)dx = f(a), \end{equation}$$ де $f(x)$ - це будь-яка неперервна функція. Таким чином ми можемо видобути значення функції в будь-якій точці.

Нам потрібно застосовувати дельта-функцію Дірака в просторі $2\pi$-періодичних функцій, тому нам потрібна така послідовність імпульсів:

$$\begin{equation} \tilde{\delta}(\omega) = 2\pi\sum_{k=-\infty}^{\infty}\delta(\omega-2\pi k) \end{equation}$$

Python-код генерації зображення

import matplotlib.pyplot as plt
import numpy as np
import matplotlib.ticker as tckr

f,(ax,ax2)=plt.subplots(2,1)
ax2.axis('off')
x=np.linspace(-4*np.pi,4*np.pi,5)
y=[1,1,1,1,1]
ax.stem(x/np.pi,y, markerfmt='^')
ax.xaxis.set_major_formatter(tckr.FormatStrFormatter('%g $\pi$'))
ax.xaxis.set_major_locator(tckr.MultipleLocator(base=1.0))
plt.show()

Добре, тепер давайте розглянемо обернене перетворення Фур'є дискретного часу цієї послідовності імпульсів. Згідно з $\eqref{eq:idtft}$: $$\begin{align} \mbox{ОПФДЧ}\{\tilde{\delta}(\omega)\} &= \frac{1}{2\pi}\int_{-\pi}^{\pi}\tilde{\delta}(\omega)e^{j\omega n}d\omega\nonumber\\ &= \int_{-\pi}^{\pi}\delta(\omega)e^{j\omega n}d\omega\nonumber\\ &= e^{j\omega n}|_{\omega=0}\nonumber\\ & = 1\label{eq:idtft_one} \end{align}$$ Давайте перевіримо це за допомогою чисельного аналізу. Розглянемо частинних сум $S_k = \sum_{n=-k}^k e^{-j\omega n}$.

Python-код для відображення частинних сум

import matplotlib.pyplot as plt
import numpy as np
from matplotlib.animation import FuncAnimation

fig, ax = plt.subplots(1,1)

N = 500
x = np.linspace(-np.pi, np.pi, N)

def update(k):
    re = np.zeros(N)
    im = np.zeros(N)
    for n in range(-k, k+1):
        re += np.cos(x*n)
        im += np.sin(x*n)
    y = np.multiply(re, re) + np.multiply(im, im)
    y = np.sqrt(y)
    ax.cla()
    ax.set_title('Частинні суми k = ' + str(k))
    ax.set_xlim(-np.pi, np.pi)
    ax.set_ylim(0,100)
    ax.plot(x, y);

anim = FuncAnimation(fig, update, frames=np.arange(0, 50), interval=100)
plt.show()

З $\eqref{eq:idtft_one}$ за допомогою властивості $\eqref{eq:dtft_modulation}$ маємо: $$\begin{equation} \mbox{ПФДЧ}\{e^{j\omega_0 n}\} = \tilde{\delta}(\omega - \omega_0). \end{equation}$$ А використовуючи формулу Ейлера можна знайти ще дві: $$\mbox{ПФДЧ}\{\cos \omega_0 n\} = [\tilde\delta(\omega-\omega_0) + \tilde\delta(\omega+\omega_0)]/2$$ $$\mbox{ПФДЧ}\{\sin \omega_0 n\} = -j[\tilde\delta(\omega-\omega_0) - \tilde\delta(\omega+\omega_0)]/2$$

Вбудовування сигналу скінченної довжини

Природне спектральне представлення сигналу скінченної довжини це ДПФ ($X[k]$). Існує два способи вбудувати такий сигнал в нескіченну послідовність:

• періодичне продовження: $\tilde{x}[n] = n[n\mod N]$
• продовження зі скінченним носієм: $\bar{x}[n] = \begin{cases} x[n] & 0\le n \le N-1\\ 0 & \mbox{інакше}\\ \end{cases}$

Зараз ми розглянемо зв'язок між ДПФ скінченного сигналу $X[n]$ і ПФДЧ продовженого сигналу $\tilde{X}(e^{j\omega})$.

Випадок періодичного сигналу

Починаючи з формули $\eqref{eq:dtft}$ маємо: $$\begin{align*} \tilde{X}(e^{j\omega}) &= \sum_{n=-\infty}^{\infty}x[n]e^{-j\omega n}\\ &= \sum_{n=-\infty}^{\infty}\left(\frac{1}{N}\sum_{k=0}^{N-1}X[k]e^{j\frac{2\pi}{N} nk}\right)e^{-j\omega n}\\ &=\frac{1}{N}\sum_{k=0}^{N-1}X[k]\sum_{n=-\infty}^{\infty}e^{j\frac{2\pi}{N} nk}e^{-j\omega n}\\ &=\frac{1}{N}\sum_{k=0}^{N-1}X[k]~\mbox{ПФДЧ}\{e^{j\frac{2\pi}{N} k}\}\\ &=\frac{1}{N}\sum_{k=0}^{N-1}X[k]\tilde{\delta}(w - \frac{2\pi}{N} k)\\ \end{align*}$$

Випадок сигналу зі скінченним носієм

Спочатку розглянемо сигнал $$\bar{r}[n] = \begin{cases}1&0\le n \le N-1\\0& \mbox{інакше}\end{cases}$$ Його ПФДЧ дорівнює $\sum_{n=0}^{N-1}e^{-j\omega n}$, яке ми позначатимемо $\bar{R}(e^{j\omega})$. $$\begin{align*} \bar{X}(e^{j\omega}) &= \sum_{n=-\infty}^{\infty}\bar{x}[n]e^{-j\omega n}\\ &= \sum_{n=0}^{N-1}x[n]e^{-j\omega n}\\ &= \sum_{n=0}^{N-1}\left(\frac{1}{N}\sum_{k=0}^{N-1}X[k]e^{j\frac{2\pi}{N} nk}\right)e^{-j\omega n}\\ &= \frac{1}{N}\sum_{k=0}^{N-1}X[k]\sum_{n=0}^{N-1}e^{j\frac{2\pi}{N} nk}e^{-j\omega n}\\ &= \frac{1}{N}\sum_{k=0}^{N-1}X[k]\bar{R}(e^{\omega - \frac{2\pi}{N} k})\\ \end{align*}$$ Зараз докладніше розглянемо $\bar{R}(e^{j\omega})$: $$\begin{align*} \sum_{n=0}^{N-1}e^{-j\omega n} &= \frac{1 - e^{-j\omega N}}{1 - e^{-j\omega}}\\ &=\frac{e^{-j\omega N/2}\left(e^{j\omega N/2}-e^{-j\omega N/2}\right)}{e^{-j\omega N}\left(e^{j\omega/2}-e^{-j\omega/2}\right)}\\ &=\frac{\sin(\omega N/2)}{\sin(\omega / 2)}e^{-j\omega (N-1)/2}\\ \end{align*}$$