fft()函数详解

fft()函数简单到发指，一般使用时就两个参数fft(nparray,n),n还可以缺省。上代码：

import numpy as np
from scipy.fftpack import fft,ifft

fft_y=fft(y)
print(fft_y)

执行结果：

[180444.84         -0.j          -1764.15187386-6325.24578909j
   2924.15553426-4550.04232005j ...    733.09073339+3427.79110955j
   2924.15553426+4550.04232005j  -1764.15187386+6325.24578909j]

这是啥？怎么用？
之前的文章里已经讲解了傅里叶变换中复数的含义。在这篇文章里，将回答以下问题：

1. n是什么？该怎么选择？
2. 为什么生成的fft序列和原始信号长度相同？它们中的每个数字是什么含义？
3. 生成的fft序列准确吗？
4. 怎么用fft来滤波？为什么滤波后ifft还原回来的值含有虚部？

一、DFT的基本原理

理解了傅里叶变换之后，DFT理解起来就不是太难了。但是DFT确实有着不同于傅里叶变换的性质。
下面是傅里叶变换的公式，这里$\lambda$是频率，t是时间：
$$\hat{f}(\lambda )=\frac{1}{\sqrt{2\pi }}{\int }_{-\infty }^{\infty }f(t)e^{-i\lambda x}d\lambda$$
当理想照进现实，我们实际面对的数据序列，只是一个个的数据采样值，没有$f(t)$函数，而且非常关键的是：它的长度都是有限的。因此我们不可能实现从$-\infty$到$\infty$的积分。所以实现离散序列的傅里叶变换的基本思路是：1）面对一个定义在有限长区间里的函数$y=f(t)$;2)使用有限的采样值作为$y_k$的值，来近似计算积分。
当分析到这里的时候，我们赫然发现，离散傅里叶变换更像是傅里叶级数，而不是傅里叶变换。我们再来看傅里叶级数的定义：
如果在区间$-\pi <=t<=\pi$上$f(t)={\sum }_{n=-\infty }^{\infty }{a}_n{e}^{int}$,则：
$$a_n=\frac{1}{2\pi }{\int }_{-\pi }^{\pi }f(t)e^{-int}dt$$
那位问了，为什么要从$-\pi$到$\pi$,还有为什么傅里叶变换的系数是$\frac{1}{\sqrt{2\pi }}$,这我也不知道，大概是当时数学家的恶趣味吧。实际上这些参数不重要，只要在逆傅里叶变换中采用对应的系数，就不会影响傅里叶变换的使用。
敲黑板，现在说重要的事：

1. 傅里叶级数定义在有限的时间长度上，这和离散傅里叶变换面对的情况是一致的，前面说了，时间长度是可以变换出来的。
2. 傅里叶级数中的n是整数 { ⋯ − 2 , − 1 , 0 , 1 , 2 , …   } \{\dots-2,-1,0,1,2,\dots\} {⋯−2,−1,0,1,2,…}。n为0的时候代表直流分量，n为自然数的时候，代表在$2\pi$时间里有n个周期的简谐振荡分量，就是频率为$\frac{n}{2\pi }$。n为负数的时候是负频率，没有意义，此时$a_n$的值为$a_{-n}$的共轭。离散傅里叶变换中的采样频率不一定为整数，但是这个也是可以转换出来的。
3. 现在说最重要的事，和傅里叶变换完全不同，傅里叶级数的对象$f(t)$必须是周期函数。对于在一定时间段内定义的函数，要求傅里叶级数，必须进行周期延拓。对于离散傅里叶变换，所获得的采样数据无论长短，一定是有限长的，因此，要以傅里叶级数的方式计算，它也必须进行周期延拓。就是在计算的时候，认为$f(t)$是周期为给出的时间段长度的周期函数。
   在讨论离散傅里叶变换的计算之前，我们得先看看积分近似计算的基本公式——梯形法则。
   
   函数$f(x)$的积分可以近似地写成很多个分段梯形面积的和。上图只画出了一个分段。面积:$$S=(b-a)\frac{f(a)+f(b)}{2}$$
   现在我们来讨论根据采样值近似计算离散傅里叶变换。我们来讨论一个在$(0,2\pi )$上的周期函数（为什么又成了$(0,2\pi )$了？你给我拿负时间上的采样结果？别忘了周期延拓，$(0,2\pi )$和$(-\pi ,\pi )$上的积分结果是一样的），在整个周期中，采样次数为n，则梯形法则中的$(b-a)=\frac{2\pi }{n}$,然后用$Y_j$来表示第j个采样点的结果，就是$f(\frac{2j\pi }{n})$,则有：
   $$\frac{1}{2\pi }{\int }_0^{2\pi }F(t)dt\approx \frac{1}{2\pi }\frac{2\pi }{n}[\frac{Y_0}{2}+Y_1+\cdots +Y_{n-1}+\frac{Y_n}{2}]$$
   注意我们讨论的$F(t)$是以$2\pi$为周期的，所以有$Y_0=Y_n$,所以上式可以进一步化简为：
   $$\frac{1}{2\pi }{\int }_0^{2\pi }F(t)dt\approx \frac{1}{n}\sum _{j=0}^{n-1}{Y}_j$$
   注意这里那个$2\pi$被悄悄地约掉了。
   我们用这个式子来计算第k个傅里叶系数的近似值，有：
   $$a_k=\frac{1}{2\pi }{\int }_0^{2\pi }f(t)e^{-ikt}dt\approx \frac{1}{n}\sum _{j=0}^{n-1}f(\frac{2\pi j}{n})e^{\frac{-2\pi ijk}{n}}$$
   注意式中那个$f(\frac{2\pi j}{n})$不就是第j个采样点的值$y_j$吗？
   所以，离散傅里叶变换的定义就隆重出场了：
   $${\hat{y}}_k=\sum _{j=0}^{n-1}{y}_j{e}^{\frac{-2\pi ijk}{n}}$$
   累死我了。
   那位又说了，$\frac{1}{n}$,$\frac{1}{n}$跑哪里去了。佛系的笑容，再说一遍，系数什么的不重要，逆变换的时候改一下系数就行了。这些数学家把系数变来变去，是为了玩我们吗？
   很完美的公式是吗？错！里面蕴含着巨大的问题！
   我们来看里面的k和j。k是频率参数，是周期为$\frac{T}{k}$的频率分量，j是第j个时间点。它们在实际应用中是什么关系？
   实际应用中，我们面对的是一个序列$y_0,y_1,\dots ,y_{n-1}$。一般我们根本就不知道面前的数据是不是周期函数，一般也不会是，所以我们假定它不是周期函数，必须用周期延拓成周期函数。那么它对应的频率分量有多少呢？傅里叶级数告诉我们有无限个。但我们只有n个采样值，奈奎斯特定理又告诉我们，我们只能求出$\frac{n}{2}$个频率分量，换句话说只有前$\frac{n}{2}$个频率分量是有意义的。
   一般而言，通过fft()函数计算离散傅里叶变换的时候，给出的结果中的频率分量也是n个，和时域上是一样的。但是，只有前$\frac{n}{2}$个是有意义的。那么，后面的是什么呢？用离散傅里叶变换的定义公式，我们可以证明：
   $${\hat{y}}_{n-k}=\overline{{\hat{y}}_k},0<=k<=n-1$$
   很可怕是不是，这些看起来一模一样的浓眉大眼的家伙，竟然只是前$\frac{n}{2}$个分量的共轭！要是用它们来滤波什么的，肯定是要掉坑里的。
   等等，这个共轭什么的，怎么和傅里叶变换中负频率的部分有点像？不是有点像，它们其实就是从$-\frac{n}{2}$到-1上的频率分量。它们的作用也是一样的，没有它们，你调用ifft()逆离散傅里叶变换就会出现一堆的复数值。
   需要注意的是，离散傅里叶变换的结果和傅里叶变换的值是有差别的，一是高频分量的缺失，二是系数的不同：<，第三，还有个Gibbs效应的问题，这个在本文的最后一部分讨论。
   这一节的最后，说一下fft算法。离散傅里叶变换叫做dft算法，fft的名称是快速傅里叶变换，它是dft算法在计算方法上的一个改进，通过把dft中的奇数项和偶数项分别计算，它可以把$n^2$次的乘法运算减少到大约$5n{\log }_{2}n$次。它要求n是偶数，同时它是一种可以迭代的算法，所以在n是2的整数幂的时候效率是最高的，就是1024,2048什么的。

二、fft()代码的注解

关于fft()具体使用的文章很少。这里把fft()函数注解的原文翻译放在这里，肯定比我之前写的要准确得多。

def fft(x, n=None, axis=-1, overwrite_x=False):
    """
	返回实数或者复数序列的离散傅里叶变换结果。
	返回的复数数组包括``y(0), y(1),..., y(n-1)``，其中
    ``y(j) = (x * exp(-2*pi*sqrt(-1)*j*np.arange(n)/n)).sum()``.
    参数
    ----------
    x : 数组
        用来做傅里叶变换的数组。
    n : 整数，可选
        傅里叶变换的长度。如果``n < x.shape[axis]``，`x`被截断。如果
         ``n > x.shape[axis]``，x被填充0。默认值是``n = x.shape[axis]``。
    axis : 整数，可选
        fft沿哪个轴计算；默认值是沿最后一个轴 (i.e., ``axis=-1``)。
    overwrite_x : 布尔量，可选
        如果为True,x的内容可能被摧毁；默认值False。

    返回值
    -------
    z : 复数ndarray
		成员::
            [y(0),y(1),..,y(n/2),y(1-n/2),...,y(-1)]        如果n是偶数
            [y(0),y(1),..,y((n-1)/2),y(-(n-1)/2),...,y(-1)]  如果n是奇数
        其中::
            y(j) = sum[k=0..n-1] x[k] * exp(-sqrt(-1)*j*k* 2*pi/n), j = 0..n-1

    See Also
    --------
    ifft : Inverse FFT
    rfft : FFT of a real sequence

    注意
    -----
    结果的组成形式是这样的“标准”：如果``A = fft(a, n)``，那么``A[0]``是
    直流分量，``A[1:n/2]``包含正频率分量，``A[n/2:]``包含负频率分量而且是
    用负频率减小的顺序。举一个8个点的例子，结果中的频率是 [0, 1, 2, 3, -4, 
    -3, -2, -1]。如果想把fft的输出调整成直流分量居中的形式，像 [-4, -3, -2,
     -1,  0,  1,  2,  3]这样，使用`fftshift`。

    single和double精度都支持。半精度输入被转换成single精度。非浮点输入被
    转换成double精度。long-double精度输入不支持。（这是啥意思？）
    
    这个函数在n是2的整数幂的时候效率最高，n是质数时候效率最低。

    注意如果x是实值的，那么``A[j] == A[n-j].conjugate()``（这个conjugate应
    该是取共轭的意思）。如果x是实值而且n是偶数，那么``A[n/2]``也是实数。

    如果x是实数，那个实数FFT逻辑会被自动应用，这大概能减半计算时间。使用
    rfft能稍微提高一点效率，它做同样的计算，但只输出对称频谱的一半。如果数
    据既是实数又是对称的，那么dct能够进一步提高一倍的效率，它只使用一半的
    数据生成一半的频谱。

    Examples
    例子
    --------
    >>> from scipy.fftpack import fft, ifft
    >>> x = np.arange(5)
    >>> np.allclose(fft(ifft(x)), x, atol=1e-15)  # within numerical accuracy.
    True

    """

三、使用fft()来滤波

fft()主要的用途大概就是分析频谱和滤波了。分析频谱就不说了，画出来就得。我们来看看滤波。从前文我们已经知道，从Y(0)到Y(n/2)是我们需要的频谱分量。其中Y(0)是直流分量，Y(1)是最低的频率分量，指在整个时间段内只有一个周期的分量，然后Y(2)Y(3)等依次增加，Y(n/2)是最高的频率。
我们先来看看傅里叶变换的逆变换能不能得到原来的结果，这里y是个(14000,)的序列：

fft_y=fft(y)
y_i = ifft(fft_y)
print(y[:10])
print(y_i[:10])
[14.69 14.67 14.69 14.64 14.68 14.71 14.71 14.74 14.74 14.75]
[14.69+1.13306190e-16j 14.67-4.38696698e-17j 14.69+3.30275489e-16j
 14.64+5.13218941e-16j 14.68-8.50113630e-17j 14.71-7.32429990e-17j
 14.71+2.86627864e-16j 14.74+4.43565819e-16j 14.74-6.77553252e-17j
 14.75-1.23678845e-16j]

可以看到数据基本被恢复了。逆变换的实部都对，虚部都小于1e-15，这应该是计算误差造成的。所以逆变换后需要舍弃过小的虚部。
现在我们来滤波，把高频部分去掉，只保留前100个频率分量（1个直流分量，99个正频率分量）。

fft_y=fft(y)
fft_y[100:]=0
y_i = ifft(fft_y)
print(y[:10])
print(y_i[:10])
[14.69 14.67 14.69 14.64 14.68 14.71 14.71 14.74 14.74 14.75]
[12.94888776-3.05519145j 12.97335191-3.05470399j 12.99777097-3.0536972j
 13.02212976-3.05217328j 13.04641322-3.05013494j 13.07060638-3.04758537j
 13.09469443-3.04452829j 13.11866268-3.04096789j 13.14249661-3.03690883j
 13.16618189-3.03235628j]

结果，一堆的复数。为什么？因为滤波的时候把负频率分量全滤掉了。这样，ifft的时候就不能生成正确的结果。
正确的滤波方法是：同时滤出正频率分量和负频率分量。

fft_y=fft(y)
fft_y[100:len(fft_y)-99] = 0 
y_i = ifft(fft_y)
print(y[:10])
print(y_i[:10])
[14.69 14.67 14.69 14.64 14.68 14.71 14.71 14.74 14.74 14.75]
[13.00885838+4.87229305e-17j 13.05778668+1.38048303e-16j
 13.10662479+2.73051423e-16j 13.15534237+4.44596740e-16j
 13.20390929+1.59364585e-16j 13.25229562+3.55271368e-17j
 13.30047172+1.07596471e-16j 13.34840821+2.91322522e-16j
 13.39607608-1.42108547e-17j 13.44344665+1.64439890e-16j]

可以看到，虚部都小于1e-15。舍弃后就是实函数了。

四、Gibbs现象

理论书籍告诉我们要注意Gibbs现象。Gibbs现象是指当用有限的傅里叶级数来复原函数时，在末端会出现频率震荡的情况。

在靠近函数的不连续点的时候，傅里叶级数的合成会产生一个尖峰，随着频率分量的增加尖峰的高度会稳定在振幅的9%左右，频率分量的增加只会压缩震荡的时间宽度，但不影响尖峰的高度。
Gibbs现象产生的原因是信号的不连续。对于离散傅里叶变换而言，由于其潜在地做了周期延拓，而$y_0!=y_{n-1}$,因此在末端应该会存在Gibbs现象。这对我们可不是好消息，因为我们分析数据的目标就是找出下一步数据的变化趋势。
但是，我做的实验并没有证明这一点。

fft_y=fft(y)
y_i = ifft(fft_y)
print(y[:10])
print(y_i[:10])
print('BBB')
print(y[-10:])
print(y_i[-10:])
[1805.31 1803.82 1804.   1803.32 1803.22 1803.25 1803.5  1800.5  1800.03
 1801.  ]
[1805.31+1.07190447e-15j 1803.82+5.38895912e-15j 1804.  +1.18559131e-15j
 1803.32+4.49063009e-15j 1803.22-1.38048303e-15j 1803.25+4.60533199e-15j
 1803.5 -8.44530795e-16j 1800.5 +7.61701813e-15j 1800.03+4.88853402e-15j
 1801.  +3.91103023e-15j]
BBB
[1737.88 1737.9  1737.9  1737.88 1737.9  1737.9  1737.88 1737.88 1737.88
 1737.9 ]
[1737.88+3.10808130e-15j 1737.9 -2.96376161e-15j 1737.9 +1.02136553e-14j
 1737.88-1.92948595e-15j 1737.9 +1.39423018e-15j 1737.9 +6.51958768e-15j
 1737.88+4.39617171e-15j 1737.88+2.48603384e-15j 1737.88+6.79271175e-15j
 1737.9 -7.48742370e-16j]

末尾的数据也被很好的复原了。
我猜测原因是现在的fft()算法中已经加入了消除Gibbs现象的机制。
消除Gibbs现象的方法是继续向后补充连续的数据。而现代的算法就是向后补充对称的数据，然后将FFT算法转变成DCT算法（二维余弦变换），可以消除Gibbs现象。
我猜测现在对fft()的实现中什么蝶形算法之类的，就已经采用了这些措施，消除了Gibbs现象。有大神给讲一讲吗？

补充几句（20210403）：这个问题搞清楚了，实际上fft算法的Gibbs现象是很严重的，这是因为sin、cos函数是典型的周期函数，对非周期的数据ifft的结果会把整个数据作为一个周期看待，而且认为应该是连续的，因此会要求开头和结尾的数据趋向于一致。
但是在python的fft包里显然实现了一个功能，就是对于未经修改的频谱可以完美的反向变换出原函数来。这使得前面的实验看不到。但只要一滤波，就会原型毕露了。

好了，啰啰嗦嗦说了很多。其实这篇文章是写给自己看的，在写文章的过程中，原本没有想透的许多东西就想明白了。所以，写博文是一种有效的学习方法。有不正确的地方请指正。