概述

数学基础

微积分

复数

复平面

欧拉公式

单位根及性质

定义

傅里叶变换(FT)

作用

公式

频域的相关解释

傅里叶逆变换(IFT)

公式

离散傅里叶变换(DFT)

公式

DFT的实现

离散傅里叶逆变换(IDFT)实现

DFT的结果含义

快速傅里叶变换(FFT)

通用性

原理推导

递归版FFT (IFFT)

迭代版FFT

递归版FFT分析

数据排序

实现

成品

概述

网上看了很多关于FFT的解释,但他们大都是用于算法加速(快速计算多项式乘法),而本文将从另一个角度来说明FFT的其他作用-----信号处理领域的天才工具

首先是名词解释

缩写	英文	中文
FT	Fourier Transform	傅里叶变换
DFT	Discrete Fourier Transform	离散傅里叶变换
FFT	Fast Fourier Transform	快速傅里叶变换
IFT	Inverse Fourier Transform	傅里叶逆变换
IDFT	Inverse Discrete Fourier Transform	离散傅里叶逆变换
IFFT	Inverse Fast Fourier Transform	快速傅里叶逆变换

之后是他们的关系

FT是最基础的一种变换,IFT则是FT的逆运算

DFT是离散化的FT(计算机只能处理离散数据)

FFT则是DFT的一种快速运算方法

IDFT和IFFT则分别是DFT和FFT的逆运算

数学基础

微积分

微积分基础(高数)还是得有的,请自行收看B站视频

复数

定义 $i\times i=-1$

$z_{1}^{}=a+bi$ 就是一个复数

$z_{2}=c+di$

共轭复数: $\bar{z}=a-bi$

加减法: $z_{1}\pm z_{2}=(a+c)\pm (b+d)i$

乘法: $z_{1}\times z_{2}=(a+bi)(c+di)=ac+adi+bci+bd\cdot i\cdot i=ac+adi+bci-bd\cdot=(ac-bd)+(ad+bc)i$

复平面

我们以实轴为横轴,虚轴为纵轴,建立平面直角坐标系,这个坐标系所在的平面被称为复平面

复数形式 $z=a+bi$ 对应着复平面中唯一确定的一点

欧拉公式

$e^{(ix)}=cos(x)+isin(x)$

通过公式可以看出, $e^{ix}$ 是复平面上的一个单位圆

随着x的增加,代表红色的点沿着逆时针方向进行旋转

关于复数和复平面的详情请看 复变函数

单位根及性质

定义

$z^{n}=1$ 的根

在复数范围内的,n次单位根有n个

它们将复平面单位圆等分为n份

记作 $w_{n}^{k}$

根据定义 $w_{n}^{k}=e^{\frac{2\pi ki}{n}}=cos(\frac{2\pi k}{n})+isin(\frac{2\pi k}{n})$

单位根性质1: $w_{n}^{k}=w_{2n}^{2k}$

证明: $w_{n}^{k}=e^{\frac{2\pi k i }{n}}=e^{\frac{w\pi (wk) i}{(2n)}}=w_{2n}^{2k}$

单位根性质2: $w_{n}^{k+\frac{n}{2}}=-w_{n}^{k}$

根据欧拉公式可以知道 $e^{\pi i}=-1$

证明: $w_{n}^{k+\frac{n}{2}}=e^{\frac{2\pi (k+\frac{n}{2})i}{n}}=e^{\frac{2\pi ki}{n}}\cdot e^{\frac{\pi ni}{n}}=-e^{\frac{2\pi ki}{n}}=-w_{n}^{k}$

单位根性质3: $w_{n}^{0}=w_{n}^{n}$

带进去就行,证明略

单位根性质4: $(w_{n}^{k})^{2}=w_{n}^{2k}$

证明: $(w_{n}^{k})^{2}=(e^{\frac{2\pi ki}{n}})^{2}=e^{\frac{2\pi (2k)i}{n}}=w_{n}^{2k}$

单位根性质5: $w^{a+b}_{n}=w^{a}_{n}\cdot w^{b}_{n}$

带进去就行,证明略

另外 $w_{n}^{n-k}=e^{-\frac{2\pi ki}{n}}$

根据欧拉公式可以知道 $e^{2\pi i}=1$

证明: $w_{n}^{n-k}=e^{\frac{2\pi (n-k)i}{n}}=e^{\frac{2\pi ni}{n}}\cdot e^{-\frac{2\pi ki}{n}}=e^{-\frac{2\pi ki}{n}}$

傅里叶变换(FT)

咱不是数学家,不谈原理,只谈应用

傅里叶变换(FT)就是一种人为发明的数学工具

作用

作用就一句话:将时域信号转为频域

这也是电学老师最常挂在嘴边的一句话

通俗理解:

我们来设想这样一个场景,现在要求这个函数 $f(t)=sin(t)+sin(3t)+sin(5t)+sin(7t)+1$ 的图像

那很简单,只需要将每个点带入函数求出值,并将其画出即可,就像下图

如果这个信号是现实中的某个信号(如声音)

横轴是以时间,纵轴是现实中的物理量(如电压),则这个图像是信号在时域中的表现

而如果我们知道这个图像想要求出上式中的各项的 $\omega$ 及其系数呢?

这就是傅里叶变换所解决的事情

如下图,横轴是 f(t) 的每个asin(wt)的w,FT有对称性,忽略负半轴

纵轴与函数的a相关,则这个图像是信号在频域中的表现

所以,傅里叶变换就是一种将组合起来的信号拆解开来的工具

这是非常反直觉的,但经过数学上的严格证明

它为我们提供了一种新的信号分析的手段

公式

傅里叶变换公式:

$F(\omega )=\int_{-\infty }^{\infty}f(t)e^{-iwt}dt$

计算时只需要把复数分为实部和虚部分别积分即可

频域的相关解释

可能原意动手的小伙伴已经发现了,好像并不能直接计算出 sin(t) 的傅里叶变换

得出上式是根据傅里叶函数的定义和一个另外的构造函数(单脉冲函数)来的

具体的手算FT的过程请看积分变换

每一个函数均可以写成有限或无限个正弦函数的和,即

$f(t)=\sum_{0}^{\infty } a_{i}cos(w_{i}t+\phi _{i})$

下图中的横坐标自然是原函数 f(t)中每一个 $a_{i}cos(w_{i}t+\phi _{i})$ 的 $\omega _{i}$

这个图仅仅表现频率的性质,而还有相位信息没有表现出来,因此表现出

关于0对称的图像

纵坐标是傅里叶变换的结果的模(傅里叶变换的结果是复变函数),这里常用模值

而与幅值的对应关系是互为相反数的两个赋值的和/2 $\pi$

即w=1的幅值 $a=\frac{\pi +\pi }{2\pi }=1$

注意:本图也是网上最常见到的图(幅度谱),只表现出了频率(w),并没有表现出相位

这个图则是刚才函数的相位谱

以 $\left | \omega \right |=1$ 的两个点为例子来解释下

这两个点在幅度谱和相位谱的值分别代表这项展开式 $a_{i}cos(w_{i}t+\phi _{i})$ 中的 $a_{i}$ 和 $\phi _{i}$

所以这项可以表示为 $\frac{\pi }{2\pi }cos(-t+\frac{\pi }{2})+\frac{\pi }{2\pi }cos(t-\frac{\pi }{2})=\frac{1}{2}cos(t-\frac{\pi }{2})+\frac{1}{2}cos(t-\frac{\pi }{2})=cos(\frac{\pi }{2}-t)=sin(t)$

傅里叶逆变换(IFT)

能将时域转化为频域进行处理,这是信号处理的一大进步,

但往往我们需要将这个频域信号转化会时域信号才能让其发挥作用

如音频处理中,如果要消除某个高频(低频)的噪声,我们要得到的肯定是时域信号,频域的信号不能被人耳所识别

这就需要本节提到的傅里叶逆变换了

公式

$f(t)=\frac{1}{2\pi }\int_{-\infty }^{\infty }F(w)e^{iwt}dw$

和傅里叶变换的公式相比,是改变了系数,被积分的量由时域信号换为了频域信号,附带的算子 $e^{-iwt}$ 换成了 $e^{iwt}$

在计算上(尤其是计算机的计算)没有特别大的区别,仅仅是更改了输入输出和某一个(sin项(欧拉公式展开后))的符号而已

离散傅里叶变换(DFT)

因为计算机只能处理离散化的数据(采样不可能做到实时),因此在处理傅里叶变换时需要将其进行离散化,也就是离散傅里叶变换(DFT)

公式

$X(m)=\sum_{n=0}^{N-1}x(n)e^{-jk\frac{2\pi n}{N}}$

X(m)就是DFT之后的频域信号,因为是离散的,因此频域信号是一个个孤立的点

x(n)是采集到的时域的信号,需要每隔固定的时间进行采集

N是采集的信号和输出信号的数量

欧拉公式 $e^{(ix)}=cos(x)+isin(x)$ 代表着在复平面上一个点随着x的增加沿着单位圆而逆时针旋转

因为在连续傅里叶变换中是 $-\infty\rightarrow \infty$ 的积分,因此 $e^{-iwt}$ 中的t代表将复平面的单位圆分为无数个小份进行求和(积分)

而在DFT中,则需要对整个单位圆间隔一段距离进行采样计算,也就是 $e^{-jm\frac{2\pi n}{N}}$ 中的 $\frac{2\pi n}{N}$

DFT的实现

因为使用到了复数运算,因此我们先建立一个复数结构体

typedef struct __complex
{float real;float imaginary;
}complex;

之后是复数的基本运算

complex complex_multiplication(complex a, complex b)
{complex zj;zj.real = a.real*b.real - a.imaginary*b.imaginary;zj.imaginary = a.real*b.imaginary + a.imaginary*b.real;return zj;
}
complex complex_add(complex a, complex b)
{complex zj;zj.real = a.real + b.real;zj.imaginary = a.imaginary + b.imaginary;return zj;
}
complex complex_subtraction(complex a, complex b)
{complex zj;zj.real = a.real - b.real;zj.imaginary = a.imaginary - b.imaginary;return zj;
}

$X(k)=\sum_{n=0}^{N-1}x(n)e^{-jk\frac{2\pi n}{N}}=\sum_{n=0}^{N-1}x(n)\left [ cos(k\frac{2\pi n}{N})-jsin(k\frac{2\pi n}{N})\right ]$

根据DFT的公式可以看出

每计算一个m的值,需要经历从0到N-1的一个累加

而x(n)是复数, $e^{-jm\frac{2\pi n}{N}}$ 可以使用欧拉公式展开为复数 $e^{-jk\frac{2\pi n}{N}}=cos(k\frac{2\pi n}{N})-jsin(k\frac{2\pi n}{N})$

所以只需要进行一个复数的乘法即可

#define PI (3.1415926)
void DFT(complex *Input, complex *Output,int Len)
{complex zj,zj1;for (int k = 0; k < Len; k++){zj.real = 0;zj.imaginary = 0;for (int n = 0; n < Len; n++){zj1.real= cos(2 * PI*k*n / Len);zj1.imaginary= -sin(2 * PI*k*n / Len);zj1=complex_multiplication(zj1, Input[n]);zj = complex_add(zj, zj1);}Output[k] = zj;}
}

首先输入的是复数数组x(k),输出数组X(k),还有长度N

之后是扫描k 从0到N-1,这步是计算频域中每一个离散点

再之后就是扫描n 也是从0到N-1,是一个累加,再确定k的条件下计算 $\sum_{n=0}^{N-1}x(n)e^{-jk\frac{2\pi n}{N}}$

具体是设置中间变量 zj , zj1,他们都是复数,zj用于存储每个不同k的累加数据

而zj1则用于存储不同k时的算子 $cos(k\frac{2\pi n}{N})-jsin(k\frac{2\pi n}{N})$

之后进行复数乘法和累加到zj中后放入输出的指定位置即可

离散傅里叶逆变换(IDFT)实现

IDFT的公式是

$x(k)=\frac{1}{N}\sum_{n=0}^{N-1}X(k)e^{jk\frac{2\pi n}{N}}=\frac{1}{N}\sum_{n=0}^{N-1}x(n)\left [ cos(k\frac{2\pi n}{N})+jsin(k\frac{2\pi n}{N})\right ]$

从公式可以看出, 我们只需要将DFT的算法进行一点改变即可变为IDFT

首先是在进行复数乘法之前,将-sin项变为sin即可,也就是乘-1

之后是在输出时虚部和实部均 / N(这里是Len)

其余地方无需改动

void DFT(complex *Input, complex *Output,int Len,int inverse)
{complex zj,zj1;for (int k = 0; k < Len; k++){zj.real = 0;zj.imaginary = 0;for (int n = 0; n < Len; n++){zj1.real= cos(2 * PI*k*n / Len);zj1.imaginary = inverse *-sin(2 * PI*k*n / Len);zj1=complex_multiplication(zj1, Input[n]);zj = complex_add(zj, zj1);}if (inverse == -1){zj.real /= Len;zj.imaginary /= Len;}Output[k] = zj;}
}

DFT的结果含义

首先说明一点,DFT是关于 $x_{\frac{N}{2}}$ 对称的

根据采样定理,想要不失真获得原始信号,则采样频率需要至少为原频率的2倍,我们一般取左半为有效信号

某个信号的幅值谱如下 ,我们要根据这个计算信号的频率幅值和相位(展开成 $a_{i}cos(w_{i}t+\phi _{i})$ )的 $a_{i},w_{i},\phi _{i}$

另外因为截取的信号不是完整周期的信号,因此会出现频谱泄露

所以求相位的方法是先找到频率,在找到频率对应的那个x(k)的虚部和实部求得相位

设采样信号的频率为 $f_{0}$ ,采样点的数量为N

则 $x(k)$ 的对应的频率为 $\frac{f_{0}k}{N}$ ,分辨率为 $\frac{f_{0}}{N}$

幅值在除了直流分量( $x(0)$ )处为 $\frac{x(k)}{\frac{N}{2}}$

在直流分量( $x(0)$ )处为 $\frac{x(0)}{N}$

在上图幅值图中,取样频率为2048Hz,取样点数为2048个

上图尖峰的点对应的分别是下表

因此信号的两个cos分量的频率分别是326Hz和652Hz

他们的幅值分别是1和3

再根据找到的这两个幅值获取实部和虚部,分别为下图

则两个分量的相位均为-90°( $-\frac{\pi }{2}$ ),复平面如下图

原信号的表达式可以表示为 $f(t)=cos(2\pi \cdot 326\cdot t-\frac{\pi }{2})+3cos(2\pi \cdot 652\cdot t-\frac{\pi }{2})=sin(2048.32t)+3sin(4096.64t)$

快速傅里叶变换(FFT)

FFT是DFT的加速算法,所以在意义和结果的用法上和DFT完全相同

通用性

首先说明DFT和IDFT均可以使用FFT进行加速,对应的是FFT和IFFT

对于IDFT

原始公式 $x(k)=\frac{1}{N}\sum_{n=0}^{N-1}X(k)e^{jk\frac{2\pi n}{N}}$

中的算子可以表示为

$e^{\frac{2\pi kj}{N}}=w_{N}^{k}$

对于DFT

原始公式 $X(k)=\sum_{n=0}^{N-1}x(n)e^{-jk\frac{2\pi n}{N}}$

其中的算子

$e^{-\frac{2\pi kj}{N}}=w_{N}^{N-k}$

可以设 t=n-k

则DFT与IDFT的形式得以统一

$w_{N}^{t}$

这两者对于FFT的加速过程来说没有本质区别

原理推导

以DFT为例,IDFT仅仅是在计算时更改了符号和系数而已

$X(k)=\sum_{n=0}^{N-1}x(n)e^{-jk\frac{2\pi n}{N}}=\sum_{n=0}^{N-1}x(n)(e^{-jk\frac{2\pi }{N}})^{n}=\sum_{n=0}^{N-1}x(n)(w_{N}^{t})^{n}$

为了方便起见,我们让 $y=w_{N}^{t}$

注意:下式中的X与X(y)均是DFT输出的某一个值,如X(0),X(1)

则 $X=x_{0}y^{0}+x_{1}y^{1}+x_{2}y^{2}+x_{3}y^{3}+x_{4}y^{5}+\cdots +x_{0n-1}y^{n-1}$

之后我们根据函数的奇偶性将X分为两部分

$X(y)=(x_{0}+x_{2}y^{2}+\cdots +x_{n-2}y^{n-2})+(x_{1}y+x_{3}y^{3}+\cdots +x_{n-1}y^{n-1})=(x_{0}+x_{2}y^{2}+\cdots +x_{n-2}y^{n-2})+y(xx_{1}+x_{3}y^{2}+\cdots +x_{n-1}y^{n-2})$

之后我们构造2个函数

$X_{1}(y)=x_{0}+x_{2}y+\cdots +x_{n-2}y^{\frac{n}{2}-1}$

$X_{2}(y)=x_{1}+x_{3}y+\cdots +x_{n-1}y^{\frac{n}{2}-1}$

则可以知道 $X(y)=X_{1}(y^{2})+yX_{2}(y^{2})$

现在令 $z<\frac{N}{2}$ 构造出一个单位根 $w_{N}^{z}$

因为原来的y也是一个单位根, $y=w_{N}^{t}$ ,且 $t<N$

所以我们需要将 $w_{N}^{z}$ 和 $w_{N}^{z+\frac{N}{2}}$ 均带入 $X(y)$ 并让两个结果连接起来才与原式等价

(也就是说将整个原来的式子分为了两部分,一个是 $<\frac{N}{2}$ 的部分,另一个是 $>\frac{N}{2}$ 的部分

搞到这里大多数人(我刚学的时候) 就已经晕了,搞这玩意干啥,和DFT有啥关系?

我们回想以下怎么通过X(y)这个数学上的函数算出DFT来

我们需要将k从0遍历到N-1,也就是t从N-1遍历到0,带入X(y)中,得到从X(0)到X(N-1)的这N个数

$y=w_{N}^{t}=w_{N}^{N-k}=e^{-\frac{2\pi kj}{N}}$

但是这样直接带入,时间复杂度和DFT没有什么区别

下面就是FFT的核心了,注意看好

将构造出来的单位根带入X(y)得

$X(w_{N}^{z}) = X_{1}((w_{N}^{z})^{2})+w_{N}^{z}X_{2}((w_{N}^{z})^{2}) = X_{1}(w_{\frac{N}{2} }^{z})+w_{N}^{z}X_{2}(w_{\frac{N}{2}}^{z})$

$X(w_{N}^{z+\frac{N}{2}}) = X_{1}(w_{N}^{2z+N})+w_{N}^{z+\frac{N}{2}}(w_{N}^{2z+N}) =X_{1}(w_{N}^{2z}\cdot w_{N}^{N})-w_{N}^{z}X_{2}(w_{N}^{2z}\cdot w_{N}^{N}) =X_{1}(w_{N}^{2z})-w_{N}^{z}X_{2}(w_{N}^{2z}) =X_{1}(w_{\frac{N}{2} }^{z})-w_{N}^{z}X_{2}(w_{\frac{N}{2}}^{z})$

为了直观我们只看最后的结果

$X(w_{N}^{z}) =X_{1}(w_{\frac{N}{2} }^{z})+w_{N}^{z}X_{2}(w_{\frac{N}{2}}^{z})$

$X(w_{N}^{z+\frac{N}{2}}) =X_{1}(w_{\frac{N}{2} }^{z})-w_{N}^{z}X_{2}(w_{\frac{N}{2}}^{z})$

看出来什么了没有,这两个(前后半部分)的表达式只有一项的符号不同

也就是说,如果我们知道了 $X_{1}(w_{\frac{N}{2} }^{z}),w_{N}^{z},X_{2}(w_{\frac{N}{2}}^{z})$ 就可以同时知道 $X(w_{N}^{z}) ,X(w_{N}^{z+\frac{N}{2}})$ ,

也就知道了DFT的计算结果了

那 $X_{1}(w_{\frac{N}{2} }^{z}),X_{2}(w_{\frac{N}{2}}^{z})$ 怎么求呢?回看一下我们的定义

$X_{1}(y)=x_{0}+x_{2}y+\cdots +x_{n-2}y^{\frac{n}{2}-1}$

$X_{2}(y)=x_{1}+x_{3}y+\cdots +x_{n-1}y^{\frac{n}{2}-1}$

和X(y)的定义

$X=x_{0}y^{0}+x_{1}y^{1}+x_{2}y^{2}+x_{3}y^{3}+x_{4}y^{5}+\cdots +x_{0n-1}y^{n-1}$

原来如此, $X_{1}(w_{\frac{N}{2} }^{z}),X_{2}(w_{\frac{N}{2}}^{z})$ 原来分别是整个序列的偶数项和奇数项的DFT啊

那就可以分别对他们的奇偶数项进行DFT

以此类推,也就是分治的方法,或者可以叫它递归

而且输入的N必须得是2的整数次幂,不然就不能每次都分为数量相等奇偶数项

等到N==1的时候,只有一项的DFT就是它本身,直接返回即可

递归版FFT (IFFT)

根据上文所述原理,我们可以轻易写出递归版的FFT

操作流程如下:

判断数据数量,如果为1 则无需操作直接返回

将奇数项和偶数项分开

分别对奇数和偶数项进行FFT得出结果(递归调用)

将结果连接起来得到FFT

而IFFT, 是在计算 $w_{N}^{z}$ 时将sin项乘上了个 -1

还需要再输出的时候将每一项均除以 N ,但是在递归里不便实现

需要的话可以在FFT结束后自行计算即可

void FFT(complex *Input, complex *Output,int Len, int inverse)
{if (Len == 1)//长度为1则直接返回即可return;complex *zj_j, *zj_o, zj,zj1;zj_j = (complex *)malloc(sizeof(complex)*Len / 2);//申请动态变量,中间变量zj_o = (complex *)malloc(sizeof(complex)*Len / 2);for (int i = 0; i < Len / 2; i++)//将输入的奇偶项分开{zj_o[i] = Input[i * 2];zj_j[i] = Input[i * 2 + 1];}FFT(zj_j, zj_j, Len / 2, inverse);//对奇偶项分别进行FFTFFT(zj_o, zj_o, Len / 2, inverse);for (int i = 0; i < Len / 2; i++){zj1.real= cos(2 * PI*i / Len);zj1.imaginary= inverse *-sin(2 * PI*i / Len);//X2的系数zj = complex_multiplication(zj1, zj_j[i]);//计算X2的那一项 这样可以减少一次复数乘法Output[i] = complex_add(zj_o[i], zj);//前半部分Output[i + Len / 2]= complex_subtraction(zj_o[i], zj);//后半部分}free(zj_j);//释放中间变量free(zj_o);
}

这段程序看起来非常好,也完美达到了完美的目标

但是占用内存太大了(尤其是对于单片机等设备,内存==钱啊),得优化一下

迭代版FFT

递归版FFT分析

首先我们看一下递归版FFT干了什么

我们将这种操作叫做蝴蝶操作

$X(w_{N}^{z}) =X_{1}(w_{\frac{N}{2} }^{z})+w_{N}^{z}X_{2}(w_{\frac{N}{2}}^{z})$

$X(w_{N}^{z+\frac{N}{2}}) =X_{1}(w_{\frac{N}{2} }^{z})-w_{N}^{z}X_{2}(w_{\frac{N}{2}}^{z})$

在同一层递归中,奇偶数项的FFT的值是确定的,而 $w_{N}^{z}$ 的值是改变的
以数据量 N=8为例介绍

这图中的元素代表一层,也就是最顶层的递归

而不同颜色的中括号代表进行蝴蝶操作中使用到的不同的 $w_{N}^{z}$

他们的奇偶数项的FFT 是相同的

并且计算后将结果放回到x(k)里面

也就是说在图中,x0-x7在经历蝴蝶变换之后就为FFT的输出数据了

而分治之后(下一层的结果)的奇偶数项的FFT已经分别存入了x0-x7

之后再看下两层递归

含义和上层相同

这便是N=8的全部层了

全家福如上

因为要先算出奇偶数项的FFT,才能进行连接(蝴蝶操作)

所以计算的顺序都是从最底层的递归开始算起

那实际上进行的操作便是

这8个输入数据分别经历了3( $\log_{2}{8}$ )次蝴蝶操作

而如果记最底层操作为0层,向上依次为 i 层,则每次蝴蝶操作的变参数为 $2^{i}$ 次单位根

如果我们知道了数据最后的排序 $x_{0},x_{4} ,x_{2} ,x_{6} ,x_{1} ,x_{5}, x_{3} , x_{7}$

那只需要 $x_{0},x_{4}$ 分别经历以 $w_{2}^{0},w_{4}^{0},w_{8}^{0}$ 为变参数的3次蝴蝶操作即可得出最后结果

那现在的问题就是如何得出数据最后的排序

数据排序

咱不是数学家,也不会完备的证明,所以就只说结论了

最后的序号顺序是 位逆序置换 的结果

大白话讲:将编号转化为二进制,再从高位读到低位,写的时候从低位往高位写,就是最后结果的编号

就像下表一样 (N=8)

固然可以直接按照刚才说到倒序摆放数字得出最终编号,但是某些数学家发明了更好的解决方法,可以加快运算

规律:(应该有证明,但是咱不会)

数据的个数 $N=8=2^{3}$ ,最终编号的二进制位数是 3 ( $\log_{2}{N}$ )位

原始编号是偶数,最终编号最高位为0,奇数最高位为1

除最高位外,其他位均为原始编号/2(向下取整) 的最终编号向右移一位的结果

(首先说明:计算按原始编号0-(N-1)进行,也就是说,再计算最终编号第x个时,我们已经知道了x/2个的最终编号是什么了)

例子:想要计算原始编号是3的最终编号,也就是第3个(从0开始数)最终编号,

最终编号最高位是1,因为原始是奇数

又因为 3/2=1(向下取整) 则第3个最终编号的低位是第1个最终编号向右移1位的结果(100>>1=010)

度3个最终编号是110

c语言实现

首先是计算 $\log_{2}{N}$ ,因为是2的整数次幂,所以使用一个循环即可

之后建立一个最终编号的数组 rev ,并按照上文提到规律进行计算

int lim = 1,//用于计数的中间变量len = 0, //最终编号的二进制位数*rev;//最终编号rev = (int *)malloc(sizeof(int)*Len);rev[0] = 0;while (lim < Len)//计算log 获取最终编号的二进制位数{lim <<= 1;len++;}for (int i = 0; i < lim; i++)//获取最终编号{rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (len - 1));}

实现

首先是将输入数组按照最终编号排序

之后需要一个循环来扫描单位根(在例子中是2,4,8次单位根)

在之前这个循环中需要嵌套一个循环,它是用于扫描n次单位根第几个的(如 $w_{8}^{0},w_{8}^{1},w_{8}^{2},w_{8}^{3}$ )

在这个循环中还需嵌套一个循环,用于扫描左半得出全部答案(如 $w_{4}^{0}$ 单位根时,要扫描 $x_{0},x_{2}$ 得出 $x_{0},x_{2},x_{4},x_{6}$ 的答案)

为了减少几次三角运算,我们要用单位根性质5(数学基础中的)

来通过 $w_{8}^{0}$ 进行复数乘法计算出 $w_{8}^{1},w_{8}^{2},w_{8}^{3}$

void FFT(complex *Input, complex *Output, int Len, int inverse)
{complex zj,omega,zj1,zj2;int lim = 1,//用于计数的中间变量len = 0, //最终编号的二进制位数*rev;//最终编号rev = (int *)malloc(sizeof(int)*Len);rev[0] = 0;while (lim < Len)//计算log 获取最终编号的二进制位数{lim <<= 1;len++;}for (int i = 0; i < Len; i++)//获取最终编号并将输入按照最终编号顺序排列{rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (len - 1));Output[i] = Input[rev[i]];}for (int m = 1; m < Len; m *= 2)//单位根扫描{zj.real = cos(PI / m);//本来应该是2PI/2m 进行了约分zj.imaginary = inverse *-sin(PI / m);for (int i = 0; i < Len; i += m * 2)//扫描奇偶项{omega.real =  1;omega.imaginary = 0;for (int j = 0; j < m; j++)//扫描左半{zj1 = Output[i + j];zj2 = complex_multiplication(omega, Output[i + j + m]);Output[i + j] = complex_add(zj1,zj2);Output[i + j + m] = complex_subtraction(zj1, zj2);//蝴蝶操作,同时得出左右半omega = complex_multiplication(omega, zj);//更新单位根}}}if (inverse == -1)//IFFT时需要/N{for (int i = 0; i < Len; i++){Output[i].real /= Len;Output[i].imaginary /= Len;}}
}

成品

FFT的内容就到这了,下面给出完整代码

使用FFT的结果作为信号分析时的方法请看 DFT的结果含义

(使用目录跳转即可)

#define PI (3.1415926)
typedef struct __complex
{float real;float imaginary;
}complex;
complex complex_multiplication(complex a, complex b)
{complex zj;zj.real = a.real*b.real - a.imaginary*b.imaginary;zj.imaginary = a.real*b.imaginary + a.imaginary*b.real;return zj;
}
complex complex_add(complex a, complex b)
{complex zj;zj.real = a.real + b.real;zj.imaginary = a.imaginary + b.imaginary;return zj;
}
complex complex_subtraction(complex a, complex b)
{complex zj;zj.real = a.real - b.real;zj.imaginary = a.imaginary - b.imaginary;return zj;
}
//输入输出的数组指针,长度只能是2的整数次幂,inverse为1是FFT,-1是IFFT
void FFT(complex *Input, complex *Output, int Len, int inverse)
{complex zj,omega,zj1,zj2;int lim = 1,//用于计数的中间变量len = 0, //最终编号的二进制位数*rev;//最终编号rev = (int *)malloc(sizeof(int)*Len);rev[0] = 0;while (lim < Len)//计算log 获取最终编号的二进制位数{lim <<= 1;len++;}for (int i = 0; i < Len; i++)//获取最终编号并将输入按照最终编号顺序排列{rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (len - 1));Output[i] = Input[rev[i]];}for (int m = 1; m < Len; m *= 2)//单位根扫描{zj.real = cos(PI / m);//本来应该是2PI/2m 进行了约分zj.imaginary = inverse *-sin(PI / m);for (int i = 0; i < Len; i += m * 2)//扫描奇偶项{omega.real =  1;omega.imaginary = 0;for (int j = 0; j < m; j++)//扫描左半{zj1 = Output[i + j];zj2 = complex_multiplication(omega, Output[i + j + m]);Output[i + j] = complex_add(zj1,zj2);Output[i + j + m] = complex_subtraction(zj1, zj2);//蝴蝶操作,同时得出左右半omega = complex_multiplication(omega, zj);//更新单位根}}}if (inverse == -1)//IFFT时需要/N{for (int i = 0; i < Len; i++){Output[i].real /= Len;Output[i].imaginary /= Len;}}
}

另外:不知道为什么,网上的FFT的单位根的 sin 使用的都是正值但是我根据公式推出的和计算相位角所得到的都应该是 -sin

请大佬解答

【基础知识】从FT到FFT相关推荐

Xiaojie雷达之路---毫米波雷达基础知识---速度估计
本篇文章主要是介绍毫米波雷达的基础知识中的速度估计,对于本篇文章主要回答4个问题 1. 雷达如何估计一个物体的速度 2. 如果有多个物体,距离雷达的距离一样,但是速度不一样,会发生什么情况 3. 两个 ...
硬件工程师基础知识（http://huarm.taobao.com/ ）
硬件工程师基础知识 1. 请列举您知道的电阻.电容.电感品牌(最好包括国内.国外品牌). 电阻: 美国:AVX.VISHAY威世日本:KOA兴亚.Kyocera京瓷.muRata村田.Pan ...
【学习笔记--FMCW基础知识】
学习笔记--FMCW基础知识前言 mmWave测距原理 mmWave区分多个物体 mmWave的距离分辨率(Range Solution) mmWave的最大测量距离前言由于工作原因需要了解TI ...
计算机基础知识与程序设计二,计算机基础与程序设计.doc
计算机基础与程序设计.doc (17页) 本资源提供全文预览,点击全文预览即可全文预览,如果喜欢文档就下载吧,查找使用更方便哦! 14.9 积分 <计算机基础与稈序设计>是高等教冇H学考 ...
查找元器件_电子元器件-常用电子元器件基础知识大全
电子元器件-常用电子元器件基础知识大全一.电阻器(R) 简称电阻,是指具有一定技术性能的在电路中专起电阻作用的元件,可用来调节电路中的电流和电压,或者作为电路中的负载.1.电阻的参数: a. 阻值: ...
学习笔记：Java 并发编程①_基础知识入门
若文章内容或图片失效,请留言反馈. 部分素材来自网络,若不小心影响到您的利益,请联系博主删除. 视频链接:https://www.bilibili.com/video/av81461839 视频下载: ...
08.GPIO基础知识和工作原理
GPIO基础知识和工作原理参考资料 1.正点原子开发板资料 <STM32FXX开发指南:库函数版本>-第六章跑马灯实验 2.STM32FxXx官方资料: F429:<STM32F4 ...
示波器基础知识100问
www.EET-china.com 点击打开链接编者注: 本文件为电子工程专辑网站编辑部原创,电子工程专辑享有本文章完全著作权, 如需转载该文章,必须经过电子工程专辑网站编辑部同意.联系电子工程专辑 ...
Linux基础第一章：基础知识与基础命令
目录一.虚拟机的三种网卡模式二.基础知识 2.1 shell的内外部命令 2.2可以使用type命令来区分内外部命令 2.3命令行头解释 2.4文件具体信息含义 2.5命令行格式三.基础命令 ...
【RAC】RAC相关基础知识
[RAC]RAC相关基础知识 1.CRS简介从Oracle 10G开始,oracle引进一套完整的集群管理解决方案--Cluster-Ready Services,它包括集群连通性.消息和锁. ...

【基础知识】从FT到FFT

概述

数学基础

微积分

复数

复平面

欧拉公式

单位根及性质

定义

傅里叶变换(FT)

作用

公式

频域的相关解释

傅里叶逆变换(IFT)

公式

离散傅里叶变换(DFT)

公式

DFT的实现

离散傅里叶逆变换(IDFT)实现

DFT的结果含义

快速傅里叶变换(FFT)

通用性

原理推导

递归版FFT (IFFT)

迭代版FFT

递归版FFT分析

数据排序

实现

成品

【基础知识】从FT到FFT相关推荐

最新文章

热门文章