当前位置: 首页 > news >正文

通过矩阵从整体角度搞懂快速傅里叶变换原理

离散傅里叶变换公式

公式

f[k]=∑n=0N−1g[n]e−i(2π/N)kn,其中(0<=n<N)f[k]=\sum_{n=0}^{N-1}g[n]e^{-i(2\pi/N)kn}, 其中(0<=n<N) f[k]=n=0N1g[n]ei(2π/N)kn,其中(0<=n<N)

逆变换公式

g[n]=1N∑k=0N−1f[k]ei(2π/N)kn,其中(0<=k<N)g[n]=\frac{1}{N}\sum_{k=0}^{N-1}f[k]e^{i(2\pi/N)kn}, 其中(0<=k<N) g[n]=N1k=0N1f[k]ei(2π/N)kn,其中(0<=k<N)

快速傅里叶变换

从以上公式看,如果直接按照公式来求离散傅里叶变换,其时间复杂度是O(N^2)

快速傅里叶变换就是一种能在O(n*log(n))时间复杂度内进行傅里叶变换及其逆变换的算法

离散傅里叶变换公式矩阵表示

G=[g[0]g[1]⋮g[n−1]]F=[f[0]f[1]⋮f[n−1]]G= \begin{bmatrix} g[0] \\ g[1] \\ \vdots \\ g[n-1] \end{bmatrix} \\\ \\\ F= \begin{bmatrix} f[0] \\ f[1] \\ \vdots \\ f[n-1] \end{bmatrix} \\\ G=g[0]g[1]g[n1]  F=f[0]f[1]f[n1] 

E[i]=[e−i(2π/N)∗0∗0e−i(2π/N)∗0∗1⋯e−i(2π/N)∗0∗(n−1)e−i(2π/N)∗1∗095⋯e−i(2π/N)∗1∗(n−1)⋮⋮⋱⋱⋮⋮⋱e−i(2π/N)∗(n−1)∗(n−1)]E[i]=\begin{bmatrix} e^{-i(2\pi/N)*0*0} & e^{-i(2\pi/N)*0*1} & \cdots & e^{-i(2\pi/N)*0*(n-1)} \\ e^{-i(2\pi/N)*1*0} & 95 & \cdots & e^{-i(2\pi/N)*1*(n-1)} \\ \vdots & \vdots & \ddots & \ddots \\ \vdots & \vdots & \ddots & e^{-i(2\pi/N)*(n-1)*(n-1)} \\ \end{bmatrix} E[i]=ei(2π/N)00ei(2π/N)10ei(2π/N)0195ei(2π/N)0(n1)ei(2π/N)1(n1)ei(2π/N)(n1)(n1)

则离散傅里叶公式可改写为

F=E[i]∗GF=E[i]*G F=E[i]G

逆变换可改写为

G=1NE[−i]∗FG=\frac{1}{N}E[-i]*F G=N1E[i]F

注意: E[i]里的 i 是一个变量,i 即正虚数单位,代表正变换,-i 代表逆变换,下同。

递归求解

E[i][n]=[e−i(2π/N)∗k∗0e−i(2π/N)∗k∗1…e−i(2π/N)∗k∗(n−1)]E[i][n]0n=[e−i(2π/N)∗k∗0e−i(2π/N)∗k∗2…e−i(2π/N)∗k∗(n−2)]E[i][n]1n=[e−i(2π/N)∗k∗1e−i(2π/N)∗k∗3…e−i(2π/N)∗k∗(n−1)]E[i][n]=\begin{bmatrix} e^{-i(2\pi/N)*k*0} & e^{-i(2\pi/N)*k*1} & \dots & e^{-i(2\pi/N)*k*(n-1)}\\ \end{bmatrix} \\\ \\\ E[i][n]^{0n}=\begin{bmatrix} e^{-i(2\pi/N)*k*0} & e^{-i(2\pi/N)*k*2} & \dots & e^{-i(2\pi/N)*k*(n-2)}\\ \end{bmatrix} \\\ \\\ E[i][n]^{1n}=\begin{bmatrix} e^{-i(2\pi/N)*k*1} & e^{-i(2\pi/N)*k*3} & \dots & e^{-i(2\pi/N)*k*(n-1)}\\ \end{bmatrix} \\ E[i][n]=[ei(2π/N)k0ei(2π/N)k1ei(2π/N)k(n1)]  E[i][n]0n=[ei(2π/N)k0ei(2π/N)k2ei(2π/N)k(n2)]  E[i][n]1n=[ei(2π/N)k1ei(2π/N)k3ei(2π/N)k(n1)]

E[i]=[E[0][n]0nE[0][n]1n⋮⋮E[n−1][n]0nE[n−1][n]1n]\\ E[i]= \begin{bmatrix} E[0][n]^{0n} & E[0][n]^{1n} \\ \vdots & \vdots \\ E[n-1][n]^{0n} & E[n-1][n]^{1n} \\ \end{bmatrix} \\ E[i]=E[0][n]0nE[n1][n]0nE[0][n]1nE[n1][n]1n

将E[i]矩阵竖着再切一刀,平分两半,用数学语言描述就是,

E00(n/2)=[E[0][n]0n⋮E[n/2−1][n]0n]E01(n/2)=[E[0][n]1n⋮E[n/2−1][n]1n]E10(n/2)=[E[n/2][n]0n⋮E[n−1][n]0n]E11(n/2)=[E[n/2][n]1n⋮E[n−1][n]1n]\\ E_{00(n/2)}= \begin{bmatrix} E[0][n]^{0n} \\ \vdots \\ E[n/2-1][n]^{0n}\\ \end{bmatrix} \\\ \\\ E_{01(n/2)}= \begin{bmatrix} E[0][n]^{1n} \\ \vdots \\ E[n/2-1][n]^{1n}\\ \end{bmatrix} \\\ \\\ E_{10(n/2)}= \begin{bmatrix} E[n/2][n]^{0n} \\ \vdots \\ E[n-1][n]^{0n}\\ \end{bmatrix} \\\ \\\ E_{11(n/2)}= \begin{bmatrix} E[n/2][n]^{1n} \\ \vdots \\ E[n-1][n]^{1n}\\ \end{bmatrix} \\\\ E00(n/2)=E[0][n]0nE[n/21][n]0n  E01(n/2)=E[0][n]1nE[n/21][n]1n  E10(n/2)=E[n/2][n]0nE[n1][n]0n  E11(n/2)=E[n/2][n]1nE[n1][n]1n

于是

E[i]=[E00(n/2)E01(n/2)E10(n/2)E11(n/2)]E[i]=\begin{bmatrix} E_{00(n/2)} & E_{01(n/2)} \\ E_{10(n/2)} & E_{11(n/2)} \\ \end{bmatrix} E[i]=[E00(n/2)E10(n/2)E01(n/2)E11(n/2)]

再令

w[k]=e−i(2π/N)∗kw[k]=e^{-i(2\pi/N)*k} w[k]=ei(2π/N)k

W0(n/2)=[w[0]00…00w[1]0…000w[2]…0⋮⋮⋮⋮0000…w[n/2−1]]W1(n/2)=[w[n/2]00…00w[n/2+1]0…000w[n/2+2]…0⋮⋮⋮⋮0000…w[n−1]]W^{0(n/2)}=\begin{bmatrix} w[0] & 0 & 0 & \dots & 0 \\ 0 & w[1] & 0 & \dots & 0 \\ 0 & 0 & w[2] & \dots & 0 \\ \vdots & \vdots & \vdots & \vdots & 0 \\ 0 & 0 & 0 & \dots & w[n/2-1] \\ \end{bmatrix} \\\ \\\ W^{1(n/2)}=\begin{bmatrix} w[n/2] & 0 & 0 & \dots & 0 \\ 0 & w[n/2+1] & 0 & \dots & 0 \\ 0 & 0 & w[n/2+2] & \dots & 0 \\ \vdots & \vdots & \vdots & \vdots & 0 \\ 0 & 0 & 0 & \dots & w[n-1] \\ \end{bmatrix} W0(n/2)=w[0]0000w[1]0000w[2]00000w[n/21]  W1(n/2)=w[n/2]0000w[n/2+1]0000w[n/2+2]00000w[n1]

于是

E[i]=[E00(n/2)E01(n/2)E10(n/2)E11(n/2)]=[E00(n/2)W0(n/2)∗E00(n/2)E10(n/2)W1(n/2)∗E10(n/2)]=[E00(n/2)W0(n/2)∗E00(n/2)−E00(n/2)−W1(n/2)∗E00(n/2)]=[E00(n/2)W0(n/2)∗E00(n/2)−E00(n/2)W0(n/2)∗E00(n/2)]E[i] =\begin{bmatrix} E_{00(n/2)} & E_{01(n/2)} \\ E_{10(n/2)} & E_{11(n/2)} \\ \end{bmatrix} \\\ \\\ =\begin{bmatrix} E_{00(n/2)} & W^{0(n/2)}*E_{00(n/2)} \\ E_{10(n/2)} & W^{1(n/2)}*E_{10(n/2)} \\ \end{bmatrix} \\\ \\\ =\begin{bmatrix} E_{00(n/2)} & W^{0(n/2)}*E_{00(n/2)} \\ -E_{00(n/2)} & -W^{1(n/2)}*E_{00(n/2)} \\ \end{bmatrix} \\\ \\\ =\begin{bmatrix} E_{00(n/2)} & W^{0(n/2)}*E_{00(n/2)} \\ -E_{00(n/2)} & W^{0(n/2)}*E_{00(n/2)} \\ \end{bmatrix} E[i]=[E00(n/2)E10(n/2)E01(n/2)E11(n/2)]  =[E00(n/2)E10(n/2)W0(n/2)E00(n/2)W1(n/2)E10(n/2)]  =[E00(n/2)E00(n/2)W0(n/2)E00(n/2)W1(n/2)E00(n/2)]  =[E00(n/2)E00(n/2)W0(n/2)E00(n/2)W0(n/2)E00(n/2)]
上面这一步就是整个转化的核心了。

再往后,令

G0(n/2)=[g[0]g[2]⋮g[n−2]]G1(n/2)=[g[1]g[3]⋮g[n−1]]G^{0(n/2)}= \begin{bmatrix} g[0] \\ g[2] \\ \vdots \\ g[n-2] \end{bmatrix} \\\ \\\ G^{1(n/2)}= \begin{bmatrix} g[1] \\ g[3] \\ \vdots \\ g[n-1] \end{bmatrix} \\\ G0(n/2)=g[0]g[2]g[n2]  G1(n/2)=g[1]g[3]g[n1] 

F=GE[i]F=[E00(n/2)W0(n/2)∗E00(n/2)−E00(n/2)W0(n/2)∗E00(n/2)]∗[G0(n/2)G1(n/2)]F=GE[i] \\\ F = \begin{bmatrix} E_{00(n/2)} & W^{0(n/2)}*E_{00(n/2)} \\ -E_{00(n/2)} & W^{0(n/2)}*E_{00(n/2)} \\ \end{bmatrix} * \begin{bmatrix} G^{0(n/2)} \\ G^{1(n/2)} \\ \end{bmatrix} \\ \\\\ F=GE[i] F=[E00(n/2)E00(n/2)W0(n/2)E00(n/2)W0(n/2)E00(n/2)][G0(n/2)G1(n/2)]
到这一步,可以看到,我们已经成功将一个 n 阶的离散傅里叶变换降为了 n/2 阶,到这里就实现了O(n*logn)时间复杂度的快速傅里叶算法。
逆变换的和正变换的大同小异,就不在赘述了。

如果想要更简洁的形式,更深刻的理解离散傅里叶变换,可以进一步化简。

M=E00(n/2)∗G0(n/2)N=E00(n/2)∗G1(n/2)Z=W0(n/2)M=E_{00(n/2)}*G^{0(n/2)} \\\\ N=E_{00(n/2)}*G^{1(n/2)} \\\\ Z=W^{0(n/2)} \\\\ M=E00(n/2)G0(n/2)N=E00(n/2)G1(n/2)Z=W0(n/2)

F=[M+Z∗N−M+Z∗N]F = \begin{bmatrix} M+Z*N \\ -M+Z*N \\ \end{bmatrix} F=[M+ZNM+ZN]

最后,上代码

public class FFT {/*** 傅里叶变换* @param x* @return*/public static Complex[] fft(Complex[] x){return fftTrans(x, true);}/*** 逆傅里叶变换* @param x* @return*/public static Complex[] ifft(Complex[] x) {int N = x.length;Complex[] y = fftTrans(x, false);for (int n = 0; n < N; n++) {y[n] = y[n].divides(N);}return y;}public static Complex[] fftTrans(Complex[] x, boolean dire) {int N = x.length;Complex[] y = new Complex[N];if (N == 1) {y[0] = x[0];return y;}Complex[] x0 = new Complex[N/2];for (int n = 0; n < N; n += 2) {x0[n/2] = x[n];}Complex[] X0 = fftTrans(x0, dire);Complex[] x1 = new Complex[N/2];for (int n = 1; n < N; n += 2) {x1[n/2] = x[n];}Complex[] X1 = fftTrans(x1, dire);for (int k = 0; k < N / 2; k++) {int ci = -1;if (!dire) {ci=-ci;}Complex wnk = getWnk(N, k, ci);Complex wnkX1 = wnk.multiply(X1[k]);y[k] = X0[k].plus(wnkX1);y[k+N/2] = X0[k].minus(wnkX1);}return y;}private static Complex getWnk(int N, int k, int n) {double T = 2 * Math.PI / N;double kth = T * k * n;Complex wk = new Complex(Math.cos(kth), Math.sin(kth));return wk;}}
http://www.lryc.cn/news/56186.html

相关文章:

  • 【C++从0到1】25、C++中嵌套使用循环
  • FastDFS与Nginx结合搭建文件服务器,并内网穿透实现公网访问
  • 密集场景下的行人跟踪替代算法,头部跟踪算法 | CVPR 2021
  • Matlab与ROS(1/2)---服务端和客户端数据通信(五)
  • 数字化转型的避坑指南:细说数字化转型十二大坑
  • pt05Encapsulationinherit
  • 面向对象编程(基础)9:封装性(encapsulation)
  • fate-serving-server增加取数逻辑并源码编译
  • 循环队列、双端队列 C和C++
  • 正则表达式(语法+例子)
  • Properties和IO流集合的方法
  • python 生成器、迭代器、动态新增属性及方法
  • Java处理JSON
  • 58-Map和Set练习-LeetCode692前k个高频单词
  • 线程生命周期及五种状态
  • OBCP第八章 OB运维、监控与异常处理-灾难恢复
  • 亚马逊云科技Serverless Data:数字经济下的创新动能
  • 【Ruby学习笔记】15.Ruby 异常
  • 聊聊MySQL主从延迟
  • 【C++从0到1】19、C++中多条件的if语句
  • 【多微电网】计及碳排放的基于交替方向乘子法(ADMM)的多微网电能交互分布式运行策略研究(Matlab代码实现)
  • Linux(centos7)安装防火墙firewalld及开放端口相关命令
  • Linux部署.Net Core Web项目
  • 【C++】STL之stack、queue的使用和模拟实现+优先级队列(附仿函数)+容器适配器详解
  • 第⑦讲:Ceph集群RGW对象存储核心概念及部署使用
  • 从异步到promise
  • Linux系统中进行JDK环境的部署
  • Leetcode.1033 移动石子直到连续
  • 【Java】在SpringBoot中使用事务注解(@Transactional)时需要注意的点
  • 找到序列最高位的1和最高位的0并输出位置