离散傅里叶变换 (DFT)、快速傅里叶变换 (FFT)

离散傅里叶变换 (DFT)

离散傅里叶变换的基

  • 对于周期为 2 π 2\pi 2π 的函数而言,它的傅里叶变换
    f ( t ) = ∑ n = − ∞ ∞ C n e i n t s . t . C n = 1 T ( ∫ 0 T f ( t ) e − i n t d t ) f ( t ) = f ( t + 2 π ) \begin{aligned} f(t)&=\sum_{n=-\infty}^\infty C_ne^{int}\\ s.t.\quad C_n&=\frac{1}{T}\left( \int_{0}^{T} f(t) e^{-int}d t\right)\\ f(t)&=f(t+2\pi) \end{aligned} f(t)s.t.Cnf(t)=n=Cneint=T1(0Tf(t)eintdt)=f(t+2π)但在实际应用中,我们可能并不知道 f ( t ) f(t) f(t),而是只能得到 f ( t ) f(t) f(t) N N N 个采样点,也就是向量 [ f 0 , f 1 , . . . , f N − 1 ] = [ f ( 0 2 π N ) , f ( 1 2 π N ) , . . . , f ( ( N − 1 ) 2 π N ) ] [f_0,f_1,...,f_{N-1}]=[f(0\frac{2\pi}{N}),f(1\frac{2\pi}{N}),...,f((N-1)\frac{2\pi}{N})] [f0,f1,...,fN1]=[f(0N2π),f(1N2π),...,f((N1)N2π)] (可以假定这 N N N 个采样点都是在同一周期内的等距离采样). 离散傅里叶变换就是通过上述向量找到傅里叶系数 C n C_n Cn
    在这里插入图片描述
  • 下面用傅里叶展开写出 f ( t ) f(t) f(t) 在采样点 2 π N \frac{2 \pi}{N} N2π 的函数值
    f ( 2 π N ) = ⋯ + C − 2 e − 2 2 π i N + C − 1 e − 1 2 π i N + C 0 e 0 2 π i N + C 1 e 1 2 π i N + C 2 e 2 2 π i N + ⋯ + C N − 1 e ( N − 1 ) 2 π i N + C N e N 2 π i N + C N + 1 e ( N + 1 ) 2 π i N + ⋯ \left.\begin{aligned} f\left(\frac{2 \pi}{N}\right)&=\cdots+ C_{-2} e^{-2 \frac{2 \pi i}{N}}+C_{-1} e^{-1 \frac{2 \pi i}{N}}+C_{0} e^{0 \frac{2 \pi i}{N}}+C_{1} e^{1 \frac{2 \pi i}{N}}+C_{2} e^{2 \frac{2 \pi i}{N}}+\cdots+C_{N-1} e^{(N-1) \frac{2 \pi i}{N}}+C_{N} e^{N \frac{2 \pi i}{N}}+C_{N+1} e^{(N+1) \frac{2 \pi i}{N}}+\cdots \\ \end{aligned}\right. f(N2π)=+C2e2N2πi+C1e1N2πi+C0e0N2πi+C1e1N2πi+C2e2N2πi++CN1e(N1)N2πi+CNeNN2πi+CN+1e(N+1)N2πi+注意到, e k 2 π i N e^{k\frac{2\pi i}{N}} ekN2πi 实际上是 k k k 的一个周期为 N N N 的函数
    在这里插入图片描述因此可以将相同的项进行合并,得到
    f ( 2 π N ) = ⋯ + C − 2 e − 2 2 π i N + C − 1 e − 1 2 π i N + C 0 e 0 2 π i N + C 1 e 1 2 π i N + C 2 e 2 2 π i N + ⋯ + C N − 1 e ( N − 1 ) 2 π i N + C N e N 2 π i N + C N + 1 e ( N + 1 ) 2 π i N + ⋯ = ( C 0 + C N + C − N + C 2 N + C − 2 N + ⋯   ) e 0 2 π i N   + ( C 1 + C N + 1 + C − N + 1 + C 2 N + 1 + C − 2 N + 1 ⋯   ) e 1 2 π i N   + ( C 2 + C N + 2 + C − N + 2 + C 2 N + 2 + C − 2 N + 2 ⋯   ) e 2 2 π i N   ⋯   + ( C N − 1 + C 2 N − 1 + C − 1 + C 3 N − 1 + C − N − 1 ⋯   ) e ( N − 1 ) 2 π i N \left.\begin{aligned} f\left(\frac{2 \pi}{N}\right)&=\cdots+ C_{-2} e^{-2 \frac{2 \pi i}{N}}+C_{-1} e^{-1 \frac{2 \pi i}{N}}+C_{0} e^{0 \frac{2 \pi i}{N}}+C_{1} e^{1 \frac{2 \pi i}{N}}+C_{2} e^{2 \frac{2 \pi i}{N}}+\cdots+C_{N-1} e^{(N-1) \frac{2 \pi i}{N}}+C_{N} e^{N \frac{2 \pi i}{N}}+C_{N+1} e^{(N+1) \frac{2 \pi i}{N}}+\cdots \\ &= \left(C_{0}+C_{N}+C_{-N}+C_{2 N}+C_{-2 N}+\cdots\right) e^{0 \frac{2 \pi i}{N}} \\ &\ +\left(C_{1}+C_{N+1}+C_{-N+1}+C_{2 N+1}+C_{-2 N+1} \cdots\right) e^{1 \frac{2 \pi i}{N}} \\ &\ +\left(C_{2}+C_{N+2}+C_{-N+2}+C_{2 N+2}+C_{-2 N+2} \cdots\right) e^{2 \frac{2 \pi i}{N}} \\ &\ \cdots \\ &\ +\left(C_{N-1}+C_{2 N-1}+C_{-1}+C_{3 N-1}+C_{-N-1} \cdots\right) e^{(N-1) \frac{2 \pi i}{N}} \end{aligned}\right. f(N2π)=+C2e2N2πi+C1e1N2πi+C0e0N2πi+C1e1N2πi+C2e2N2πi++CN1e(N1)N2πi+CNeNN2πi+CN+1e(N+1)N2πi+=(C0+CN+CN+C2N+C2N+)e0N2πi +(C1+CN+1+CN+1+C2N+1+C2N+1)e1N2πi +(C2+CN+2+CN+2+C2N+2+C2N+2)e2N2πi  +(CN1+C2N1+C1+C3N1+CN1)e(N1)N2πi w = e 2 π i N w=e^{\frac{2\pi i}{N}} w=eN2πi ( w 0 = w N = 1 w^0=w^N=1 w0=wN=1),则有
    f ( 2 π N ) = ( C 0 + C N + C − N + C 2 N + C − 2 N + ⋯   ) w 0   + ( C 1 + C N + 1 + C − N + 1 + C 2 N + 1 + C − 2 N + 1 ⋯   ) w 1   + ( C 2 + C N + 2 + C − N + 2 + C 2 N + 2 + C − 2 N + 2 ⋯   ) w 2   ⋯   + ( C N − 1 + C 2 N − 1 + C − 1 + C 3 N − 1 + C − N − 1 ⋯   ) w ( N − 1 ) \left.\begin{aligned} f\left(\frac{2 \pi}{N}\right) &= \left(C_{0}+C_{N}+C_{-N}+C_{2 N}+C_{-2 N}+\cdots\right) w^{0} \\ &\ +\left(C_{1}+C_{N+1}+C_{-N+1}+C_{2 N+1}+C_{-2 N+1} \cdots\right) w^{1} \\ &\ +\left(C_{2}+C_{N+2}+C_{-N+2}+C_{2 N+2}+C_{-2 N+2} \cdots\right) w^{2} \\ &\ \cdots \\ &\ +\left(C_{N-1}+C_{2 N-1}+C_{-1}+C_{3 N-1}+C_{-N-1} \cdots\right) w^{(N-1)} \end{aligned}\right. f(N2π)=(C0+CN+CN+C2N+C2N+)w0 +(C1+CN+1+CN+1+C2N+1+C2N+1)w1 +(C2+CN+2+CN+2+C2N+2+C2N+2)w2  +(CN1+C2N1+C1+C3N1+CN1)w(N1)也就是说, f ( 2 π N ) f\left(\frac{2 \pi}{N}\right) f(N2π) 的函数值,只需要 N N N 个基就能得到,不需要无穷多个基,只要得到这 N N N 个基的 N N N 个系数即可
  • 按照同样的方法可以写出 f ( t ) f(t) f(t) 在采样点 k 2 π N k\frac{2\pi}{N} kN2π 的函数值
    f ( k 2 π N ) = ∑ n = − ∞ ∞ C n w n k = ( C 0 + C N + C − N + C 2 N + C − 2 N + ⋯   ) w 0 k   + ( C 1 + C N + 1 + C − N + 1 + C 2 N + 1 + C − 2 N + 1 ⋯   ) w 1 k   + ( C 2 + C N + 2 + C − N + 2 + C 2 N + 2 + C − 2 N + 2 ⋯   ) w 2 k   ⋯   + ( C N − 1 + C 2 N − 1 + C − 1 + C 3 N − 1 + C − N − 1 ⋯   ) w ( N − 1 ) k \left.\begin{aligned} f\left(k\frac{2 \pi}{N}\right)&=\sum_{n=-\infty}^\infty C_nw^{nk}\\ &= \left(C_{0}+C_{N}+C_{-N}+C_{2 N}+C_{-2 N}+\cdots\right) w^{0k} \\ &\ +\left(C_{1}+C_{N+1}+C_{-N+1}+C_{2 N+1}+C_{-2 N+1} \cdots\right) w^{1k} \\ &\ +\left(C_{2}+C_{N+2}+C_{-N+2}+C_{2 N+2}+C_{-2 N+2} \cdots\right) w^{2k} \\ &\ \cdots \\ &\ +\left(C_{N-1}+C_{2 N-1}+C_{-1}+C_{3 N-1}+C_{-N-1} \cdots\right) w^{(N-1)k} \end{aligned}\right. f(kN2π)=n=Cnwnk=(C0+CN+CN+C2N+C2N+)w0k +(C1+CN+1+CN+1+C2N+1+C2N+1)w1k +(C2+CN+2+CN+2+C2N+2+C2N+2)w2k  +(CN1+C2N1+C1+C3N1+CN1)w(N1)k由于 w k w^k wk 的周期性, w 0 k , . . . , w ( N − 1 ) k w^{0k},...,w^{(N-1)k} w0k,...,w(N1)k 仍然对应着 w 0 , . . . , w N − 1 w^0,...,w^{N-1} w0,...,wN1 中的一项,因此 f ( k 2 π N ) f\left(k\frac{2 \pi}{N}\right) f(kN2π) 的函数值都只需要 N N N 个基就能得到,基为
    w k = e k 2 π i N , k = 0 , 1 , . . . , N − 1 w^k=e^{k\frac{2\pi i}{N}},\quad\quad k=0,1,...,N-1 wk=ekN2πi,k=0,1,...,N1

离散傅里叶变换

  • 现在我们知道, f ( t ) f(t) f(t) 在采样点处的值 f ( k 2 π N ) f\left(k\frac{2 \pi}{N}\right) f(kN2π) 其实是 N N N 个基的线性组合,并且函数值 f ( k 2 π N ) f\left(k\frac{2 \pi}{N}\right) f(kN2π) N N N 个基是已知的,因此我们就可以列出关于 N N N 个系数 的 N N N 个方程,组成线性方程组,最终就可以解出 N N N 个系数 ( C 0 + C N + C − N + C 2 N + C − 2 N + ⋯   ) , ( C 1 + C N + 1 + C − N + 1 + C 2 N + 1 + C − 2 N + 1 ⋯   ) , . . . , ( C N − 1 + C 2 N − 1 + C − 1 + C 3 N − 1 + C − N − 1 ⋯   ) \left(C_{0}+C_{N}+C_{-N}+C_{2 N}+C_{-2 N}+\cdots\right),\\\left(C_{1}+C_{N+1}+C_{-N+1}+C_{2 N+1}+C_{-2 N+1} \cdots\right),\\...,\\\left(C_{N-1}+C_{2 N-1}+C_{-1}+C_{3 N-1}+C_{-N-1} \cdots\right) (C0+CN+CN+C2N+C2N+),(C1+CN+1+CN+1+C2N+1+C2N+1),...,(CN1+C2N1+C1+C3N1+CN1) { w 0 ( C 0 + C N + C − N + C 2 N + C − 2 N + ⋯   ) + w 0 ( C 1 + C N + 1 + C − N + 1 + C 2 N + 1 + C − 2 N + 1 ⋯   ) + ⋯ + w 0 ( C N − 1 + C 2 N − 1 + C − 1 + C 3 N − 1 + C − N − 1 ⋯   ) = f ( 0 ) w 0 ( C 0 + C N + C − N + C 2 N + C − 2 N + ⋯   ) + w 1 ( C 1 + C N + 1 + C − N + 1 + C 2 N + 1 + C − 2 N + 1 ⋯   ) + ⋯ + w ( N − 1 ) ( C N − 1 + C 2 N − 1 + C − 1 + C 3 N − 1 + C − N − 1 ⋯   ) = f ( 2 π N ) ⋯ w 0 k ( C 0 + C N + C − N + C 2 N + C − 2 N + ⋯   ) + w 1 k ( C 1 + C N + 1 + C − N + 1 + C 2 N + 1 + C − 2 N + 1 ⋯   ) + ⋯ + w ( N − 1 ) k ( C N − 1 + C 2 N − 1 + C − 1 + C 3 N − 1 + C − N − 1 ⋯   ) = f ( k 2 π N ) ⋯ w 0 ( N − 1 ) ( C 0 + C N + C − N + C 2 N + C − 2 N + ⋯   ) + w 1 ( N − 1 ) ( C 1 + C N + 1 + C − N + 1 + C 2 N + 1 + C − 2 N + 1 ⋯   ) + ⋯ + w ( N − 1 ) ( N − 1 ) ( C N − 1 + C 2 N − 1 + C − 1 + C 3 N − 1 + C − N − 1 ⋯   ) = f ( ( N − 1 ) 2 π N ) \left\{\begin{aligned} &w^{0} \left(C_{0}+C_{N}+C_{-N}+C_{2 N}+C_{-2 N}+\cdots\right) +w^{0}\left(C_{1}+C_{N+1}+C_{-N+1}+C_{2 N+1}+C_{-2 N+1} \cdots\right) + \cdots +w^{0}\left(C_{N-1}+C_{2 N-1}+C_{-1}+C_{3 N-1}+C_{-N-1} \cdots\right) =f\left(0\right)\\ &w^{0} \left(C_{0}+C_{N}+C_{-N}+C_{2 N}+C_{-2 N}+\cdots\right) +w^{1}\left(C_{1}+C_{N+1}+C_{-N+1}+C_{2 N+1}+C_{-2 N+1} \cdots\right) + \cdots +w^{(N-1)}\left(C_{N-1}+C_{2 N-1}+C_{-1}+C_{3 N-1}+C_{-N-1} \cdots\right) =f\left(\frac{2 \pi}{N}\right) \\&\cdots \\ &w^{0k} \left(C_{0}+C_{N}+C_{-N}+C_{2 N}+C_{-2 N}+\cdots\right) +w^{1k}\left(C_{1}+C_{N+1}+C_{-N+1}+C_{2 N+1}+C_{-2 N+1} \cdots\right) + \cdots +w^{(N-1)k}\left(C_{N-1}+C_{2 N-1}+C_{-1}+C_{3 N-1}+C_{-N-1} \cdots\right) =f\left(k\frac{2 \pi}{N}\right) \\&\cdots \\ &w^{0(N-1)} \left(C_{0}+C_{N}+C_{-N}+C_{2 N}+C_{-2 N}+\cdots\right) +w^{1(N-1)}\left(C_{1}+C_{N+1}+C_{-N+1}+C_{2 N+1}+C_{-2 N+1} \cdots\right) + \cdots +w^{(N-1)(N-1)}\left(C_{N-1}+C_{2 N-1}+C_{-1}+C_{3 N-1}+C_{-N-1} \cdots\right) =f\left((N-1)\frac{2 \pi}{N}\right) \end{aligned} \right. w0(C0+CN+CN+C2N+C2N+)+w0(C1+CN+1+CN+1+C2N+1+C2N+1)++w0(CN1+C2N1+C1+C3N1+CN1)=f(0)w0(C0+CN+CN+C2N+C2N+)+w1(C1+CN+1+CN+1+C2N+1+C2N+1)++w(N1)(CN1+C2N1+C1+C3N1+CN1)=f(N2π)w0k(C0+CN+CN+C2N+C2N+)+w1k(C1+CN+1+CN+1+C2N+1+C2N+1)++w(N1)k(CN1+C2N1+C1+C3N1+CN1)=f(kN2π)w0(N1)(C0+CN+CN+C2N+C2N+)+w1(N1)(C1+CN+1+CN+1+C2N+1+C2N+1)++w(N1)(N1)(CN1+C2N1+C1+C3N1+CN1)=f((N1)N2π)
  • 但我们想解的其实是傅里叶系数 C k C_k Ck,因此离散傅里叶假设 f ( t ) f(t) f(t) 的傅里叶级数里只包含 C 0 , C 1 , . . , C N − 1 C_0,C_1,..,C_{N-1} C0,C1,..,CN1 (只要 N N N 足够大,就可以非常接近原函数),此时我们就可以写出如下线性方程组并解出 C 0 , C 1 , . . , C N − 1 C_0,C_1,..,C_{N-1} C0,C1,..,CN1
    { w 0 C 0 + w 0 C 1 + ⋯ + w 0 C N − 1 = f ( 0 ) w 0 C 0 + w 1 C 1 + ⋯ + w ( N − 1 ) C N − 1 = f ( 2 π N ) ⋯ w 0 k C 0 + w 1 k C 1 + ⋯ + w ( N − 1 ) k C N − 1 = f ( k 2 π N ) ⋯ w 0 ( N − 1 ) C 0 + w 1 ( N − 1 ) C 1 + ⋯ + w ( N − 1 ) ( N − 1 ) C N − 1 = f ( ( N − 1 ) 2 π N ) \left\{\begin{aligned} &w^{0} C_{0} +w^{0}C_{1} + \cdots +w^{0}C_{N-1} =f\left(0\right)\\ &w^{0} C_{0} +w^{1}C_{1} + \cdots +w^{(N-1)}C_{N-1} =f\left(\frac{2 \pi}{N}\right) \\&\cdots \\ &w^{0k} C_{0} +w^{1k}C_{1} + \cdots +w^{(N-1)k}C_{N-1} =f\left(k\frac{2 \pi}{N}\right) \\&\cdots \\ &w^{0(N-1)} C_{0} +w^{1(N-1)}C_{1} + \cdots +w^{(N-1)(N-1)}C_{N-1} =f\left((N-1)\frac{2 \pi}{N}\right) \end{aligned} \right. w0C0+w0C1++w0CN1=f(0)w0C0+w1C1++w(N1)CN1=f(N2π)w0kC0+w1kC1++w(N1)kCN1=f(kN2π)w0(N1)C0+w1(N1)C1++w(N1)(N1)CN1=f((N1)N2π)写成矩阵形式即为
    [ f 0 f 1 f 2 f 3 ⋮ f N − 1 ] = [ 1 1 1 1 ⋯ 1 1 w w 2 w 3 ⋯ w N − 1 1 w 2 w 4 w 6 ⋯ w 2 ( N − 1 ) 1 w 3 w 6 w 9 ⋯ w 3 ( N − 1 ) ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ 1 w N − 1 w 2 ( N − 1 ) w 3 ( N − 1 ) ⋯ w ( N − 1 ) 2 ] [ C 0 C 1 C 2 C 3 ⋮ C N − 1 ] \left[\begin{array}{c} f_{0} \\ f_{1} \\ f_{2} \\ f_{3} \\ \vdots \\ f_{N-1} \end{array}\right]=\left[\begin{array}{cccccc} 1 & 1 & 1 & 1 & \cdots & 1 \\ 1 & w & w^{2} & w^{3} & \cdots & w^{N-1} \\ 1 & w^{2} & w^{4} & w^{6} & \cdots & w^{2(N-1)} \\ 1 & w^{3} & w^{6} & w^{9} & \cdots & w^{3(N-1)} \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 1 & w^{N-1} & w^{2(N-1)} & w^{3(N-1)} & \cdots & w^{(N-1)^{2}} \end{array}\right]\left[\begin{array}{c} C_{0} \\ C_{1} \\ C_{2} \\ C_{3} \\ \vdots \\ C_{N-1} \end{array}\right] f0f1f2f3fN1 = 111111ww2w3wN11w2w4w6w2(N1)1w3w6w9w3(N1)1wN1w2(N1)w3(N1)w(N1)2 C0C1C2C3CN1 其中,由 w = e 2 π i N w=e^{\frac{2\pi i}{N}} w=eN2πi 组成的 N × N N\times N N×N 矩阵 F N ∗ F_N^{*} FN 即为傅里叶逆变换矩阵,由 f f f 向量到 C C C 向量的变换即为离散傅里叶变换,由 C C C 向量到 f f f 向量的变换即为离散傅里叶逆变换
  • F N ∗ F_N^{*} FN 的共轭转置为 F N F_N FN,则有
    F N ∗ F N = N I N F_N^*F_N=NI_N FNFN=NIN因此,
    F N ∗ − 1 = 1 N F N F_N^{*-1}=\frac{1}{N}F_N FN∗−1=N1FN可以将离散傅里叶变换写为
    1 N F N f = C 1 N [ 1 1 1 1 ⋯ 1 1 w ˉ w ˉ 2 w ˉ 3 ⋯ w ˉ N − 1 1 w ˉ 2 w ˉ 4 w ˉ 6 ⋯ w ˉ 2 ( N − 1 ) 1 w 3 w 6 w 9 ⋯ w 3 ( N − 1 ) ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ 1 w ˉ N − 1 w ˉ 2 ( N − 1 ) w ˉ 3 ( N − 1 ) ⋯ w ˉ ( N − 1 ) 2 ] [ f 0 f 1 f 2 f 3 ⋮ f N − 1 ] = [ C 0 C 1 C 2 C 3 ⋮ C N − 1 ] \frac{1}{N}F_Nf=C \\ \frac{1}{N}\left[\begin{array}{cccccc} 1 & 1 & 1 & 1 & \cdots & 1 \\ 1 & \bar{w} & \bar{w}^2 & \bar{w}^3 & \cdots & \bar{w}^{N-1} \\ 1 & \bar{w}^2 & \bar{w}^4 & \bar{w}^6 & \cdots & \bar{w}^{2(N-1)} \\ 1 & w^3 & w^6 & w^9 & \cdots & w^{3(N-1)} \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 1 & \bar{w}^{N-1} & \bar{w}^{2(N-1)} & \bar{w}^{3(N-1)} & \cdots & \bar{w}^{(N-1)^2} \end{array}\right]\left[\begin{array}{c} f_0 \\ f_1 \\ f_2 \\ f_3 \\ \vdots \\ f_{N-1} \end{array}\right]=\left[\begin{array}{c} C_0 \\ C_1 \\ C_2 \\ C_3 \\ \vdots \\ C_{N-1} \end{array}\right] N1FNf=CN1 111111wˉwˉ2w3wˉN11wˉ2wˉ4w6wˉ2(N1)1wˉ3wˉ6w9wˉ3(N1)1wˉN1wˉ2(N1)w3(N1)wˉ(N1)2 f0f1f2f3fN1 = C0C1C2C3CN1 其中, w ˉ = e − 2 π i N \bar w=e^{-\frac{2\pi i}{N}} wˉ=eN2πi w w w 的共轭。设 X [ k ] = N C k X[k]=NC_k X[k]=NCk x [ k ] = f k x[k]=f_k x[k]=fk,则有
    X [ k ] = ∑ n = 0 N − 1 x [ n ] e − k 2 π n i N X[k]=\sum_{n=0}^{N-1} x[n] e^{-k \frac{2 \pi n i}{N}} X[k]=n=0N1x[n]ekN2πni即为书上的 DFT 公式;而逆离散傅里叶变换 IDFT 则为
    x [ n ] = 1 N ∑ k = 0 N − 1 X [ k ] e k 2 π n i N x[n]=\frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{k\frac{2\pi n i}{N}} x[n]=N1k=0N1X[k]ekN2πni
def DFT_slow(x):
    x = np.asarray(x, dtype=float)  # 输入信号采样点
    N = x.shape[0]  # 采样点个数

    spectro = []    # 信号频谱
    for k in range(N):  # 依次计算 N 个傅里叶系数
        real = 0    # 傅里叶系数的实部
        img = 0     # 傅里叶系数的虚部
        for n in range(N):  # 遍历所有采样点
            t = 2 * math.pi * k * n / N
            real += math.cos(t) * x[n]  # 累加实部
            img -= math.sin(t) * x[n]   # 累加虚部
        
        spectro.append(real + img * 1j)
    return np.asarray(spectro)
def DFT_slow(x):
    # 用矩阵乘积的形式实现 DFT
    x = np.asarray(x, dtype=float)  # 输入信号采样点
    N = x.shape[0]  # 采样点个数
    n = np.arange(N)
    k = n.reshape((N, 1))
    M = np.exp(-2j * np.pi * k * n / N) # 傅里叶变换矩阵
    return np.dot(M, x)

  • 假设周期函数为 f ( t ) = 1 + e i t + e i 2 t + e i 3 t \begin{aligned} f(t)&=1+e^{it}+e^{i2t}+e^{i3t} \end{aligned} f(t)=1+eit+ei2t+ei3t
  • 采样 N = 4 N=4 N=4 个点时,可以计算得到采样点处的函数值
    f ( 0 ) = 4 f ( 1 2 π 4 ) = f ( 2 2 π 4 ) = f ( 3 2 π 4 ) = 0 f(0)=4 \quad f\left(1 \frac{2 \pi}{4}\right)=f\left(2 \frac{2 \pi}{4}\right)=f\left(3 \frac{2 \pi}{4}\right)=0 f(0)=4f(142π)=f(242π)=f(342π)=0 w = e 2 π i N = i w=e^{\frac{2\pi i}{N}}=i w=eN2πi=i,由离散傅里叶变换可以解得 [ C 0 , C 1 , C 2 , C 3 ] = [ 1 , 1 , 1 , 1 ] [C_0,C_1,C_2,C_3]=[1,1,1,1] [C0,C1,C2,C3]=[1,1,1,1]
    [ 1 1 1 1 1 w w 2 w 3 1 w 2 w 4 w 6 1 w 3 w 6 w 9 ] − 1 [ 4 0 0 0 ] = [ 1 1 1 1 ] \left[\begin{array}{cccc} 1 & 1 & 1 & 1 \\ 1 & w & w^2 & w^3 \\ 1 & w^2 & w^4 & w^6 \\ 1 & w^3 & w^6 & w^9 \end{array}\right]^{-1}\left[\begin{array}{l} 4 \\ 0 \\ 0 \\ 0 \end{array}\right]=\left[\begin{array}{l} 1 \\ 1 \\ 1 \\ 1 \end{array}\right] 11111ww2w31w2w4w61w3w6w9 1 4000 = 1111
  • 采样 N = 3 N=3 N=3 个点时,可以计算得到采样点处的函数值
    f ( 0 ) = 4 f ( 2 π 3 ) = 1 f ( 2 2 π 3 ) = 1 f(0)=4 \quad f\left(\frac{2 \pi}{3}\right)=1 \quad f\left(2 \frac{2 \pi}{3}\right)=1 f(0)=4f(32π)=1f(232π)=1 w = e 2 π i N = − 1 2 + 3 2 i w=e^{\frac{2\pi i}{N}}=-\frac{1}{2}+\frac{\sqrt3}{2}i w=eN2πi=21+23 i,由离散傅里叶变换可以解得 [ C 0 + C 3 , C 1 , C 2 ] = [ 2 , 1 , 1 ] [C_0+C_3,C_1,C_2]=[2,1,1] [C0+C3,C1,C2]=[2,1,1]
    [ 1 1 1 1 w w 2 1 w 2 w 4 ] − 1 [ 4 1 1 ] = [ 2 1 1 ] \left[\begin{array}{ccc} 1 & 1 & 1 \\ 1 & w & w^2 \\ 1 & w^2 & w^4 \end{array}\right]^{-1}\left[\begin{array}{l} 4 \\ 1 \\ 1 \end{array}\right]=\left[\begin{array}{l} 2 \\ 1 \\ 1 \end{array}\right] 1111ww21w2w4 1 411 = 211 可以看出,由于采样点个数太少,我们无法解出高频系数 C 3 C_3 C3,只能解出 C 0 + C 3 C_0+C_3 C0+C3,也就是把高频信号当成低频信号

  • 假设周期函数为 f ( t ) = 1 + cos ⁡ t + cos ⁡ 2 t = 1 2 e − i 2 t + 1 2 e − i t + 1 + 1 2 e i t + 1 2 e i 2 t \begin{aligned} f(t)&=1+\cos t+\cos 2t=\frac{1}{2} e^{-i 2 t}+\frac{1}{2} e^{-i t}+1+\frac{1}{2} e^{i t}+\frac{1}{2} e^{i 2 t} \end{aligned} f(t)=1+cost+cos2t=21ei2t+21eit+1+21eit+21ei2t
  • 采样 N = 6 N=6 N=6 个点时,由离散傅里叶变换可以解得 [ C 0 , C 1 , C 2 , C 3 , C − 2 , C − 1 ] = [ 1 , 0.5 , 0.5 , 0 , 0.5 , 0.5 ] [C_0,C_1,C_2,C_3,C_{-2},C_{-1}]=[1,0.5,0.5,0,0.5,0.5] [C0,C1,C2,C3,C2,C1]=[1,0.5,0.5,0,0.5,0.5]
  • 采样 N = 4 N=4 N=4 个点时,由离散傅里叶变换可以解得 [ C 0 , C 1 , C 2 + C − 2 , C − 1 ] = [ 1 , 0.5 , 1 , 0.5 ] [C_0,C_1,C_2+C_{-2},C_{-1}]=[1,0.5,1,0.5] [C0,C1,C2+C2,C1]=[1,0.5,1,0.5]
  • 采样 N = 5 N=5 N=5 个点时,由离散傅里叶变换可以解得 [ C 0 , C 1 , C 2 , C − 2 , C − 1 ] = [ 1 , 0.5 , 0.5 , 0.5 , 0.5 ] [C_0,C_1,C_2,C_{-2},C_{-1}]=[1,0.5,0.5,0.5,0.5] [C0,C1,C2,C2,C1]=[1,0.5,0.5,0.5,0.5]

快速傅里叶变换 (FFT)

  • 快速傅里叶变换就是快速计算如下矩阵乘积的算法:
    F N f = [ 1 1 1 1 ⋯ 1 1 w ˉ w ˉ 2 w ˉ 3 ⋯ w ˉ N − 1 1 w ˉ 2 w ˉ 4 w ˉ 6 ⋯ w ˉ 2 ( N − 1 ) 1 w 3 w 6 w 9 ⋯ w 3 ( N − 1 ) ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ 1 w ˉ N − 1 w ˉ 2 ( N − 1 ) w ˉ 3 ( N − 1 ) ⋯ w ˉ ( N − 1 ) 2 ] [ f 0 f 1 f 2 f 3 ⋮ f N − 1 ] F_Nf=\left[\begin{array}{cccccc} 1 & 1 & 1 & 1 & \cdots & 1 \\ 1 & \bar{w} & \bar{w}^2 & \bar{w}^3 & \cdots & \bar{w}^{N-1} \\ 1 & \bar{w}^2 & \bar{w}^4 & \bar{w}^6 & \cdots & \bar{w}^{2(N-1)} \\ 1 & w^3 & w^6 & w^9 & \cdots & w^{3(N-1)} \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 1 & \bar{w}^{N-1} & \bar{w}^{2(N-1)} & \bar{w}^{3(N-1)} & \cdots & \bar{w}^{(N-1)^2} \end{array}\right]\left[\begin{array}{c} f_0 \\ f_1 \\ f_2 \\ f_3 \\ \vdots \\ f_{N-1} \end{array}\right] FNf= 111111wˉwˉ2w3wˉN11wˉ2wˉ4w6wˉ2(N1)1wˉ3wˉ6w9wˉ3(N1)1wˉN1wˉ2(N1)w3(N1)wˉ(N1)2 f0f1f2f3fN1 其中, F N ∗ { i j } = w ˉ i j F_N^*\{ij\}=\bar w^{ij} FN{ij}=wˉij w ˉ = e − 2 π i N \bar w=e^{-\frac{2\pi i}{N}} wˉ=eN2πi
  • 当计算矩阵乘积 A x Ax Ax 时,改变 A A A 的列的顺序,同时改变 x x x 中对应行的顺序,可以使得计算结果不变。因此对于 F N f F_Nf FNf,可以把 F N F_N FN 的偶数列和 f f f 的偶数行都拿到前面得到 F N f = F ~ N f ~ F_Nf=\tilde F_N\tilde f FNf=F~Nf~。例如,对于 F 4 F_4 F4 ( w ˉ \bar w wˉ),有
    F 4 = [ 1 1 1 1 1 w ˉ w ˉ 2 w ˉ 3 1 w ˉ 2 w ˉ 4 w ˉ 6 1 w ˉ 3 w ˉ 6 w ˉ 9 ] F_4=\left[\begin{array}{cccc} 1 & 1 & 1 & 1 \\ 1 & \bar{w} & \bar{w}^2 & \bar{w}^3 \\ 1 & \bar{w}^2 & \bar{w}^4 & \bar{w}^6 \\ 1 & \bar{w}^3 & \bar{w}^6 & \bar{w}^9 \end{array}\right] F4= 11111wˉwˉ2wˉ31wˉ2wˉ4wˉ61wˉ3wˉ6wˉ9 F ~ 4 = [ 1 1 1 1 1 w ˉ 2 w ˉ w ˉ 3 1 w ˉ 4 w ˉ 2 w ˉ 6 1 w ˉ 6 w ˉ 3 w ˉ 9 ] = [ F 2 D 2 F 2 F 2 − D 2 F 2 ] \tilde F_4=\left[\begin{array}{cccc} 1 & 1 & 1 & 1 \\ 1 &\bar{w}^2& \bar{w} & \bar{w}^3 \\ 1 & \bar{w}^4 & \bar{w}^2& \bar{w}^6 \\ 1 & \bar{w}^6 & \bar{w}^3 & \bar{w}^9 \end{array}\right]=\left[\begin{array}{cc} F_2 & D_2 F_2 \\ F_2 & -D_2 F_2 \end{array}\right] F~4= 11111wˉ2wˉ4wˉ61wˉwˉ2wˉ31wˉ3wˉ6wˉ9 =[F2F2D2F2D2F2]其中, D 2 = [ 1 0 0 w ˉ ] D_2=\left[\begin{array}{cc} 1 & 0 \\ 0 & \bar{w} \end{array}\right] D2=[100wˉ] 为对角矩阵, N = 4 N=4 N=4 时有 w ˉ k + 4 = w ˉ k \bar w^{k+4}=\bar w^k wˉk+4=wˉk. 将上述结论一般化可以得到
    F 2 N ~ = [ F N D N F N F N − D N F N ] \widetilde{F_{2 N}}=\left[\begin{array}{cc} F_N & D_N F_N \\ F_N & -D_N F_N \end{array}\right] F2N =[FNFNDNFNDNFN]其中,
    D N = [ 1 w ˉ ⋱ w ˉ N − 1 ] D_N=\left[\begin{array}{cccc} 1 & & & \\ & \bar{w} & & \\ & & \ddots & \\ & & & \bar{w}^{N-1} \end{array}\right] DN= 1wˉwˉN1
  • 利用上述结论可以得到
    F N f = F ~ N f ~ = [ F N 2 D N 2 F N 2 F N 2 − D N 2 F N 2 ] [ f 0 f 2 ⋮ f N − 2 f 1 f 3 ⋮ f N − 1 ] = [ I D N 2 I − D N 2 ] [ F N 2 F N 2 ] [ f e v e n f o d d ] = [ I D N 2 I − D N 2 ] [ F N 2 f e v e n F N 2 f o d d ] \begin{aligned} F_N f=\tilde{F}_N\tilde f&=\left[\begin{array}{cc} F_{\frac{N}{2}} & D_{\frac{N}{2}} F_{\frac{N}{2}}\\ F_{\frac{N}{2}} & -D_{\frac{N}{2}} F_{\frac{N}{2}} \end{array}\right]\left[\begin{array}{c} f_0 \\ f_2 \\ \vdots \\ f_{N-2} \\ f_1 \\ f_3 \\ \vdots \\ f_{N-1} \end{array}\right] \\&=\left[\begin{array}{cc} I & D_{\frac{N}{2}}\\ I & -D_{\frac{N}{2}} \end{array}\right] \left[\begin{array}{cc} F_{\frac{N}{2}} & \\ & F_{\frac{N}{2}} \end{array}\right]\left[\begin{array}{c} f_{even} \\ f_{odd} \end{array}\right] \\&=\left[\begin{array}{cc} I & D_{\frac{N}{2}}\\ I & -D_{\frac{N}{2}} \end{array}\right] \left[\begin{array}{c} F_{\frac{N}{2}}f_{even} \\ F_{\frac{N}{2}}f_{odd} \end{array}\right] \end{aligned} FNf=F~Nf~=[F2NF2ND2NF2ND2NF2N] f0f2fN2f1f3fN1 =[IID2ND2N][F2NF2N][fevenfodd]=[IID2ND2N][F2NfevenF2Nfodd]可以看到,经过变换后,计算出后面的矩阵 [ F N 2 f e v e n F N 2 f o d d ] \left[\begin{array}{c} F_{\frac{N}{2}}f_{even} \\ F_{\frac{N}{2}}f_{odd} \end{array}\right] [F2NfevenF2Nfodd] 只需要 2 ( N 2 ) 2 2(\frac{N}{2})^2 2(2N)2 次乘法,而由于前面一个矩阵是对角矩阵组成的分块矩阵,因此计算前后矩阵的乘积一共只需要 N N N 次乘法。在使用一次变换后,我们只需要 N + 2 ( N 2 ) 2 N+2(\frac{N}{2})^2 N+2(2N)2 次乘法,远小于原来的 N 2 N^2 N2 次乘法。并且后面的 F N 2 F_{\frac{N}{2}} F2N 还可以被继续分解,一直可以分解到二维矩阵乘二维向量,如下图所示:
    在这里插入图片描述

FFT 实际上只要求 N N N 可以被分解为两个正整数的乘积 (e.g. F 9 F_9 F9 可以被分解为 9 个 3 × 3 3\times3 3×3 的矩阵),但我们一般都选择 N N N 为 2 的整数次幂,如果信号采样次数不够 2 k 2^k 2k,那我们就在信号后面补 0

  • 通过 FFT,我们可以将 DFT 原先 O ( N 2 ) O(N^2) O(N2)时间复杂度降为 O ( N log ⁡ N ) O(N\log N) O(NlogN)
def FFT(x):
    """A recursive implementation of the 1D Cooley-Tukey FFT"""
    x = np.asarray(x, dtype=float)
    N = x.shape[0]
    
    if N % 2 > 0:
        raise ValueError("size of x must be a power of 2")
    elif N <= 32:  # this cutoff should be optimized
        return DFT_slow(x)
    else:
        X_even = FFT(x[::2])
        X_odd = FFT(x[1::2])
        factor = np.exp(-2j * np.pi * np.arange(N) / N)
        return np.concatenate([X_even + factor[:N // 2] * X_odd,
                               X_even + factor[N // 2:] * X_odd])

卷积

  • 卷积本身没有实际意义,它只是一种数学操作,但是卷积是很多场景下都会遇到的一类重复操作,因此单独起名为卷积。下面介绍两个卷积经常出现的场景:线性时不变系统傅里叶级数

线性时不变系统

线性 f ( a x + b y ) = a f ( x ) + b f ( y ) f(ax+by)=af(x)+bf(y) f(ax+by)=af(x)+bf(y); 时不变:系统的属性和时间没有关系

  • 给定线性时不变系统的单位冲激响应 h [ n ] h[n] h[n] (横轴为时间)
    在这里插入图片描述我们想知道输入信号 x [ n ] x[n] x[n] 时,系统的输出信号 y [ n [ y[n[ y[n[
    在这里插入图片描述
  • 由于是线性系统,因此可以分别计算 3 个时刻输入信号的输出 ( h [ n − 1 ] h[n-1] h[n1] 表示 h [ n ] h[n] h[n] 右移 1 个时间单位)
    在这里插入图片描述将 3 个输出相加即可得到 y [ n ] y[n] y[n]
    在这里插入图片描述
  • 可以看出,计算 y [ n ] y[n] y[n] 的公式即为离散卷积公式
    在这里插入图片描述可以简写为
    y = x ∗ h y=x*h y=xh也就是说,在定义了卷积之后,我们就可以说任何一个线性时不变系统的输出就是输入信号和单位冲激响应的卷积

傅里叶级数

  • 考虑两个周期为 2 π 2\pi 2π 的周期函数
    f ( t ) = ∑ n = − ∞ ∞ c n e i n t = ∑ n = − ∞ ∞ c n z n g ( t ) = ∑ n = − ∞ ∞ d n e i n t = ∑ n = − ∞ ∞ d n z n \begin{aligned} f(t)&=\sum_{n=-\infty}^\infty c_ne^{int}=\sum_{n=-\infty}^\infty c_nz^{n}\\ g(t)&=\sum_{n=-\infty}^\infty d_ne^{int}=\sum_{n=-\infty}^\infty d_nz^{n} \end{aligned} f(t)g(t)=n=cneint=n=cnzn=n=dneint=n=dnzn其中, z = e i t z=e^{it} z=eit. 我们想要求得 f ( t ) g ( t ) f(t)g(t) f(t)g(t) 的傅里叶系数
  • 我们将 f ( t ) , g ( t ) , f ( t ) g ( t ) f(t),g(t),f(t)g(t) f(t),g(t),f(t)g(t) 都展开,可得
    f ( t ) = ⋯ + c − 2 z − 2 + c − 1 z − 1 + c 0 z 0 + c 1 z + c 2 z 2 + ⋯ g ( t ) = ⋯ + d − 2 z − 2 + d − 1 z − 1 + d 0 z 0 + d 1 z + d 2 z 2 + ⋯ f ( t ) g ( t ) = ⋯ + m − 2 z − 2 + m − 1 z − 1 + m 0 z 0 + m 1 z + m 2 z 2 + ⋯ \begin{gathered} f(t)=\cdots+c_{-2} z^{-2}+c_{-1} z^{-1}+c_0 z^0+c_1 z+c_2 z^2+\cdots \\ g(t)=\cdots+d_{-2} z^{-2}+d_{-1} z^{-1}+d_0 z^0+d_1 z+d_2 z^2+\cdots \\ f(t) g(t)=\cdots+m_{-2} z^{-2}+m_{-1} z^{-1}+m_0 z^0+m_1 z+m_2 z^2+\cdots \end{gathered} f(t)=+c2z2+c1z1+c0z0+c1z+c2z2+g(t)=+d2z2+d1z1+d0z0+d1z+d2z2+f(t)g(t)=+m2z2+m1z1+m0z0+m1z+m2z2+可知,
    m n = c 0 d n + c 1 d n − 1 + c − 1 d n + 1 + ⋯ = ∑ k = − ∞ ∞ c k d n − k m_n=c_0 d_n+c_1 d_{n-1}+c_{-1} d_{n+1}+\cdots=\sum_{k=-\infty}^{\infty} c_k d_{n-k} mn=c0dn+c1dn1+c1dn+1+=k=ckdnk因此, f ( t ) g ( t ) f(t)g(t) f(t)g(t) 的傅里叶系数即为 f ( t ) f(t) f(t) 傅里叶系数和 g ( t ) g(t) g(t) 傅里叶系数的卷积

卷积这一节后面还有些内容,以后再看

参考文献


more:

评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值