傅里叶变换 - 2

最新推荐文章于 2024-02-21 19:46:19 发布

yyx_daytoy

最新推荐文章于 2024-02-21 19:46:19 发布

阅读量956

点赞数 1

分类专栏：傅里叶变换文章标签： python fft

本文链接：https://blog.csdn.net/yyx_daytoy/article/details/113656815

版权

傅里叶变换专栏收录该内容

4 篇文章 0 订阅

订阅专栏

N点离散傅里叶变换(DFT)

说明：在实际应用中，N取值一般为 $2^{n}$ ，如果N不满足要求，一般用0将其补充成 $2^{n}$ 。例如：对音频进行傅里叶变换时，音频采样率是16kHz，分帧窗长是25ms，那么每帧音频的长度 N = 16000 * 0.025 = 400，在后面补充112个0，将其补充为 N = 512，做512点的离散傅里叶变换。当然，如果窗长是32，每帧音频的长度 N = 16000 * 0.032 = 512，正好做512点的离散傅里叶变换。

离散傅里叶变换公式：

$\dpi{150} X(k) = \sum_{n=0}^{N-1} x(n) e^{\frac{-i2\pi kn}{N}}$ k = 0,1,2,...,N-1

根据欧拉变换 $e^{ix} = cos(x) + i*sin(x)$

得到

$X(k)= \sum_{n=0}^{N-1} x(n) \left \{ cos(\frac{-2\pi kn}{N}) + i*sin(\frac{-2\pi kn}{N})\right \} = \sum_{n=0}^{N-1} x(n) \left \{ cos(\frac{2\pi kn}{N}) - i*sin(\frac{2\pi kn}{N})\right \}$ k = 0,1,2,...,N-1

用python代码描述为：

import numpy as np 
from math import pi, sin, cos


class complex:
    def __init__(self):
        self.real = 0
        self.image = 0
 

def my_fft(inpt, N):
    if len(inpt) < N:
        need_zeros = np.zeros((N-len(inpt)))
        inpt = np.hstack((inpt, need_zeros))
    result = []
    for k in range(N):
        res = complex()
        for n in range(N):
            tmp = complex()
            tmp.real = inpt[n] * cos(2*pi*k*n/N)
            tmp.image =  - inpt[n] * sin(2*pi*k*n/N)

            res.real = res.real + tmp.real
            res.image = res.image + tmp.image
        result.append(res)
    return result
    

def py_fft(inpt, N):
    fft_result = np.fft.rfft(inpt, N)
    return fft_result


if __name__ == "__main__":
    inpt = np.array([1,2,3,4,5,6,7])
    N = 8
    my_res = my_fft(inpt, N)
    py_res = py_fft(inpt, N)
    
    print("\n------------------my_fft:")
    for item in my_res:
        print("%f %f"%(item.real, item.image))
    
    print("\n------------------py_fft:")
    for item in py_res:
        print(item)

上面代码 N=8，my_fft是根据离散傅里叶变换公式写的函数，py_fft是调用python的函数，运行结果如下：

------------------my_fft:
28.000000 0.000000                   0
-9.656854 4.000000                   1
-4.000000 -4.000000                  2
1.656854 -4.000000                   3
4.000000 0.000000                    4
1.656854 4.000000                    5
-4.000000 4.000000                   6
-9.656854 -4.000000                  7

------------------py_fft:
(28+0j)                              0
(-9.65685424949+4j)                  1
(-4-4j)                              2
(1.65685424949-4j)                   3
(4+0j)                               4

从结果可以看出，my_fft 和 py_fft 的前5个输出值实部和虚部部分是一致的，但为什么 py_fft 做8点的傅里叶变换，只有5个输出呢？别急，比较一下 my_fft 输出 1和7，2和6，3和5，这几对值都是实部相同，虚部相反数。所以对于8点傅里叶变换，只要知道前5个值，后面3个值也就知道了。

将N值改为16，inpt不变，运行代码得到结果是：

------------------my_fft:
28.000000 0.000000
-0.746035 -22.075230
-9.656854 4.000000
6.441553 -0.091993
-4.000000 -4.000000
1.215301 4.050143
1.656854 -4.000000
-2.910819 2.066906
4.000000 0.000000
-2.910819 -2.066906
1.656854 4.000000
1.215301 -4.050143
-4.000000 4.000000
6.441553 0.091993
-9.656854 -4.000000
-0.746035 22.075230

------------------py_fft:
(28+0j)
(-0.746034924454-22.0752300017j)
(-9.65685424949+4j)
(6.4415530545-0.0919925532372j)
(-4-4j)
(1.21530119499+4.05014307049j)
(1.65685424949-4j)
(-2.91081932504+2.06690562202j)
(4+0j)

可以得到和N=8相同的结论

所以，之后我们说傅里叶变换，一般只取前面的有效值。N点离散傅里叶变换，只取前 N/2+1个值。

my_fft 函数中一共用到了两层for循环，时间复杂度是 O(n^2) ，为了加快离散傅里叶变换，就用到快速傅里叶变换(fast Fourier transform，FFT)