【机器学习笔记31】傅里叶变换

【参考资料】
【1】《复变函数与积分变换》
【2】《数字信号处理》

1. 傅里叶级数

定义: 设$f_T(t)$fT​(t)是以T为周期的实值函数，且在$[-\dfrac{T}{2}, \dfrac{T}{2}]$[−2T​,2T​]上满足狄利克雷条件，即：
1）连续或只有有限个第一类间断点；
备注：第一类间断点包括：
可去间断点：左右极限存在且相等，但不等于该点的值，或该点无定义；
跳跃间断点：左右极限存在但不相等
2）只有有限个极值点

那么$f_T(t)$fT​(t)在连续点处有

傅里叶级数的三角形式

$f_T(t)=\dfrac{a_0}{2} + \sum\limits_{n=1}^{+\infty}(a_n \cos n \omega_0 t + b_n \sin n \omega_0 t)$fT​(t)=2a0​​+n=1∑+∞​(an​cosnω0​t+bn​sinnω0​t)
$w_0 = \dfrac{2 \pi}{T}$w0​=T2π​

$a_n = \dfrac{2}{T}\int_{-T/2}^{T/2} f_T(t) \cos n \omega_0 tdt$an​=T2​∫−T/2T/2​fT​(t)cosnω0​tdt其中$n = 0, 1, 2, \dots$n=0,1,2,…

$b_n = \dfrac{2}{T}\int_{-T/2}^{T/2} f_T(t) \sin n \omega_0 tdt$bn​=T2​∫−T/2T/2​fT​(t)sinnω0​tdt其中$n = 0, 1, 2, \dots$n=0,1,2,…

傅里叶级数的复指数形式

$f_T(t) = \sum\limits_{-\infty}^{+\infty} c_n e^{jn \omega_0 t}$fT​(t)=−∞∑+∞​cn​ejnω0​t
$c_n = \dfrac{1}{T} \int_{-T/2}^{T/2} f_T(t) e^{-jn \omega_0 t}dt$cn​=T1​∫−T/2T/2​fT​(t)e−jnω0​tdt

---

2. 傅里叶变换

周期函数可以展开为傅里叶级数。从物理上将周期为T的函数可以由一系列$\omega_0 = \dfrac{2\pi}{T}$ω0​=T2π​为间隔的离散频率所形成的简谐波合成。当T越来越大时，取得频率间隔越来越小，当T变成无穷大时，周期函数就变成了非周期函数，此时傅里叶级数就转变称了傅里叶变换。

推导如下:

$f(t)=\lim\limits_{T \to +\infty}f_T(t)$f(t)=T→+∞lim​fT​(t)
$=\lim\limits_{T \to +\infty} \sum\limits_{n=-\infty}^{n=+\infty} [ \dfrac{1}{T} \int_{-T/2}^{T/2} f_T(\tau) e^{-jn \omega_0 \tau} d \tau ] e^{jn \omega_0 t}$=T→+∞lim​n=−∞∑n=+∞​[T1​∫−T/2T/2​fT​(τ)e−jnω0​τdτ]ejnω0​t

在周期函数下时域信号被分解为频率$\omega_0 = \dfrac{2 \pi}{T}$ω0​=T2π​的简谐波，此时T区域无穷。

因此有越来越小的$\omega_0$ω0​记作$\Delta \omega$Δω，而$n \omega_0$nω0​记作$\omega_n$ωn​, 并且$T=\dfrac{2 \pi}{\Delta \omega}$T=Δω2π​
于是有:

$f(t)=\dfrac{1}{2\pi} \lim\limits_{\Delta \omega \to 0} \sum\limits_{n=-\infty}^{n=+\infty} [ \int_{-\pi/\Delta \omega}^{\pi/\Delta \omega}f_T(\tau)e^{-j\omega_n \tau}d\tau . e^{j \omega_n t}] \Delta \omega$f(t)=2π1​Δω→0lim​n=−∞∑n=+∞​[∫−π/Δωπ/Δω​fT​(τ)e−jωn​τdτ.ejωn​t]Δω

$f(t)=\dfrac{1}{2\pi} \int\limits_{-\infty}^{+\infty}[\int\limits_{-\infty}^{+\infty} f(\tau) e^{-j\omega \tau}d\tau] e^{j\omega t}d\omega$f(t)=2π1​−∞∫+∞​[−∞∫+∞​f(τ)e−jωτdτ]ejωtdω

因此得到定理如下（傅里叶变换）：
$F(\omega)=\int\limits_{-\infty}^{+\infty} f(t) e^{-j\omega t} dt$F(ω)=−∞∫+∞​f(t)e−jωtdt
$f(t)=\dfrac{1}{2\pi}\int\limits_{-\infty}^{+\infty} F(\omega)e^{j\omega t}d \omega$f(t)=2π1​−∞∫+∞​F(ω)ejωtdω

3. 单位冲击函数

单位冲击函数定义:
1）当$t \ne 0$t̸​=0时，$\xi(t)=0$ξ(t)=0
2）$\int_{-\infty}^{+\infty} \xi(t) dt = 1$∫−∞+∞​ξ(t)dt=1

主要性质：
1 当实数域有界函数f(t)，在t=0点连续，则$\int_{-\infty}^{+\infty} \xi(t) f(t) dt = f(0)$∫−∞+∞​ξ(t)f(t)dt=f(0)
2 当实数域有界函数f(t)，在$t = t_0$t=t0​点连续，则$\int_{-\infty}^{+\infty} \xi(t - t_0) f(t) dt = f(t_0)$∫−∞+∞​ξ(t−t0​)f(t)dt=f(t0​)
3 $\xi(t)$ξ(t)是偶函数，有$\xi(t) = \xi(-t)$ξ(t)=ξ(−t)
4 $\xi(t)$ξ(t)是单位阶跃函数$\mu (t)$μ(t)的导数
$\mu(t) = \begin{cases} 1, &amp; t &gt; 0 \\ 0, &amp; t &lt; 0 \end{cases}$μ(t)={1,0,​t>0t<0​

5 $\xi(t)$ξ(t)函数的傅里叶变换

公式1：$F(\omega)=F[xi(t)] =\int_{-\infty}^{+\infty}\xi(t)e^{-j\omega t}dt = e^{-jwt}|_{t=0} = 1$F(ω)=F[xi(t)]=∫−∞+∞​ξ(t)e−jωtdt=e−jwt∣t=0​=1
公式2：$F^{-1}[1] = \dfrac{1}{2\pi} \int_{-\infty}^{+\infty} e^{-j \omega t} d \omega = \xi(t)$F−1[1]=2π1​∫−∞+∞​e−jωtdω=ξ(t) //根据逆变换可知
根据上面两个公式，我们可以作出如下推导:
对$f(t)=1$f(t)=1，其傅里叶变换（由公式2） $F(\omega)=\int_{-\infty}^{+\infty} e^{-j \omega t} dt = 2\pi \xi(t)$F(ω)=∫−∞+∞​e−jωtdt=2πξ(t)
对$f(t)=e^{j \omega_0 t}$f(t)=ejω0​t
$F(\omega) = \int_{-\infty}^{+\infty} e^{j \omega_0 t} e^{-j \omega t} dt = 2 \pi \xi(\omega_0 - \omega)$F(ω)=∫−∞+∞​ejω0​te−jωtdt=2πξ(ω0​−ω)
对$f(t)=cos \omega_0 t = 1/2 \int_{-\infty}^{+\infty}(e^{j \omega_0 t} + e^{-j \omega_0 t}) e^{-j \omega t} dt = \pi[\xi(\omega - \omega_0) + \xi(\omega + \omega_0)]$f(t)=cosω0​t=1/2∫−∞+∞​(ejω0​t+e−jω0​t)e−jωtdt=π[ξ(ω−ω0​)+ξ(ω+ω0​)]

因此在广义的傅里叶变换下周期函数仍然可以以傅里叶变换展开

4. 离散傅里叶变换

设x(n)是一个长度为M的有限长序列，则定义x(n)的N点离散傅里叶变换为:
$X(k) = DFT[x(n)]=\sum\limits_{n=0}^{N-1}W_{N}^{kn} \quad k = 0,1,2,\dots, N-1$X(k)=DFT[x(n)]=n=0∑N−1​WNkn​k=0,1,2,…,N−1

X(k)的傅里叶逆变换为:
$x(n) = IDFT[X(k)] = \dfrac{1}{N}\sum\limits_{n=0}^{N-1}X(k)W_{N}^{-kn} \quad n = 0,1,2,\dots, N-1$x(n)=IDFT[X(k)]=N1​n=0∑N−1​X(k)WN−kn​n=0,1,2,…,N−1
其中$W_N=e^{-j \dfrac{2\pi}{N}}$WN​=e−jN2π​，N是DFT变换区间长度，$N \ge M$N≥M

5. 代码实现(numpy和纯python)

#!/usr/bin/env python
# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
import matplotlib.ticker as tik
import pandas as pd
import numpy as np
import math


def _fft_code():

    plt.figure(figsize=(8, 4))

    df = pd.read_csv('./data/test.csv', sep=',')

    data = df['std_value'].astype(float)

    ax1 = plt.subplot(1,2,1)
    ax1.set_title('time series')
    plt.gca().yaxis.set_major_formatter(tik.FormatStrFormatter('%1.2f'))

    plt.plot(data,'r')

    list_len = len(data)

    y = data

    m = np.floor(list_len) #DFT中的N,要求N >= M，这里字节用序列长度

    a = np.zeros(int(m))
    b = np.zeros(int(m))
    c = np.zeros(int(m))

    N = int(list_len - 1)

    for i in range(0, int(m)- 1):
       for t in range(0, N - 1):

           a[i + 1] = a[i + 1] + y[t + 1] * np.cos(2 * np.pi * i * t/N)
           b[i + 1] = b[i + 1] + y[t + 1] * np.sin(2 * np.pi * i * t/N)

       c[i + 1] = a[i + 1] + b[i + 1]

    ax2 = plt.subplot(1,2,2)
    ax2.set_title('Fourier Result')
    plt.plot(c,'b')
    plt.show()


def _fft_np():

    df = pd.read_csv('./data/test.csv', sep=',')

    data = df['std_value']

    fig = plt.figure(figsize=(8, 4))

    

    ax1 = plt.subplot(1,2,1)
    ax1.set_title('time series')
    plt.gca().yaxis.set_major_formatter(tik.FormatStrFormatter('%1.2f'))


    plt.plot(data,'r')

    transformed = np.fft.fft(data)

    d = abs(transformed) 
    x = np.arange(len(d))

    ax2 = plt.subplot(1,2,2)
    ax2.set_title('Fourier Result')
    plt.plot(d,'b')

    plt.show()

"""
说明：

傅里叶变换代码实现，对应的笔记《傅里叶变换》

作者：fredric

日期：2018-11-18

_fft_np: 采用numpy自带傅里叶变换函数

_fft_code: 纯python代码，依据基本的DFT推导，未采用FFT


"""

if __name__ == '__main__':

    _fft_code()

    _fft_np()