离散傅里叶DFT变换
目录：
**1、**傅里叶变换的问题
2、复数的基础知识
3、傅里叶变换的基础知识
4、离散傅里叶变换
5、递推离散傅里叶
6、问题的解答
7、总结
1、傅里叶变换的问题
下面两道题关于使用傅里叶变换的， 这应该是很常见的嵌入式问题：
该系统通过 ADC(<16bit) 捕捉 50Hz 的交流电流电压信号，并以 800Hz 的频率进行采集。请计算其电流、电压幅值及其相应的功率与功率因数。
B) 上面的 50Hz 电压中， 混入了另一个 50Hz 的电压，求出这两个电压的幅值。
这两道题使用 16bit、32bit 的整数运算，不使用浮点运算，可以在 MCU 上实现。
C) 完成一个 wav 声音文件的变速不变调的程序。
2、复数的基础知识
在讲解 Fourier Transform 前，必须知道一点基本的复数知识。
在复平面上的一个点P(x, y) 用复数表示为：P = x + iy
用极坐标表示为：P = r * e^(ia)
r等于√(x²+y²)，表示点(x,y)与原点之间的距离；而θ等于arctan(x,y)，代表角度；其中e(数学-e的自然之美)表示自然常数。在这一部分展示了关键的关系式：e^(iθ)=cosθ+isinθ。这一发现源于数学的精美中的详细推导过程★
此式通过运用复数来实现角度变换与三角函数变换的强大工具。例如，在将点P转换b角度后，则新点(x₁, y₁)的角度变为a+b，并且其模长保持不变。
P1 = x1 + i y1
= r * e^[i (a+b)]
= r*e^(i a) * e^(i b)
= (x + i y) * [cos(b) + i sin(b)]
= [x * cos(b) - y * sin(b)] + i [ y * cos(b) + x * sin(b)]
3、傅里叶变换的基础知识
傅里叶变换是一种积分运算，通过访问以下链接可以获得更为详细的信息：快速傅里叶FFT变换、[基2 FFT时间抽取和频域抽取算法]( "基2 FFT时间抽取和频域抽取算法")。
4、离散傅里叶变换
1）变换公式
离散傅里叶变换即 DFT 的公式：
其中，X(k) 是第 k次谐波的复数；N为周期采样点数；x(n)为输入，n从0到N-1。
用伪代码更直观地说明：
typedef short int16;
typedef int int32;
typedef struct SComplex
{
int16 Real; //实部
int16 Image; //虚部
} Complex;
void CalculateHarmonic(Complex* X, int harmonic) //计算谐波，harmonic为x次谐波，x=1、2、3……
{
for (int i=0; i<N; i++)
{
X->Real = x(i) * cos( 2*PI* i/N * harmonic)/N; //x(i)为输入，i从0到i-1
X->Image = x(i) * sin(-2*PI* i/N * harmonic)/N; //N=16
}
}看到这里可能会觉得奇怪的是什么呢？那就是说，在信息处理领域中存在一些看似简单却暗含玄机的概念。然而在实际应用中我们常常会遇到一些看似复杂实则简单的技术原理。举个例子大家都知道计算机处理信息时都会使用到二进制数字系统而那些复杂的算法往往都建立在这个基础之上。所以掌握这些基本概念对于深入理解后续内容是非常重要的
2）代码解析
struct ffts
{
double real;
double imag;
}DFTSWN[2][16]={
{/**********一次波（基波）***************/
{1.000000e+00,-0.000000e+00}, //01
{9.238795e-01,-3.826834e-01}, //02
{7.071068e-01,-7.071068e-01}, //03
{3.826834e-01,-9.238795e-01}, //04
{-8.269461e-16,-1.000000e+00}, //05
{-3.826834e-01,-9.238795e-01}, //06
{-7.071068e-01,-7.071068e-01}, //07
{-9.238795e-01,-3.826834e-01}, //08
{-1.000000e+00,1.653892e-15}, //09
{-9.238795e-01,3.826834e-01}, //10
{-7.071068e-01,7.071068e-01}, //11
{-3.826834e-01,9.238795e-01}, //12
{2.480838e-15,1.000000e+00}, //13
{3.826834e-01,9.238795e-01}, //14
{7.071068e-01,7.071068e-01}, //15
{9.238795e-01,3.826834e-01}, //16
},
{ /***********二次谐波***************/
{1.000000e+00,-0.000000e+00}, //01
{7.071068e-01,-7.071068e-01}, //02
{-8.269461e-16,-1.000000e+00}, //03
{-7.071068e-01,-7.071068e-01}, //04
{-1.000000e+00,1.653892e-15}, //05
{-7.071068e-01,7.071068e-01}, //06
{2.480838e-15,1.000000e+00}, //07
{7.071068e-01,7.071068e-01}, //08
{1.000000e+00,-3.307784e-15}, //09
{7.071068e-01,-7.071068e-01}, //10
{-4.134730e-15,-1.000000e+00}, //11
{-7.071068e-01,-7.071068e-01}, //12
{-1.000000e+00,4.961676e-15}, //13
{-7.071068e-01,7.071068e-01}, //14
{4.900444e-15,1.000000e+00}, //15
{7.071068e-01,7.071068e-01}, //16
}
};
void DFT_16(unsigned char id)
{
int i;
double a1,a2;
int N=16;
DFTResults[id].real=0;
DFTResults[id].imag=0;
for(i=0;i<N;i++) //DFTSWN
{
a1=DFTdatas[i]*DFTSWN[id][i].real;
a2=DFTdatas[i]*DFTSWN[id][i].imag;
DFTResults[id].real+=a1;
DFTResults[id].imag+=a2;
}
DFTResults[id].real=DFTResults[id].real*2/N;
DFTResults[id].imag=DFTResults[id].imag*2/N;
}（1）初始化：变量DFTResults[id].real和变量DFTResults[id].imag被初始化为零值，并用于存储所对应的频率分量（其中id表示k值）的实部与虚部。
（2）旋转因子乘法累加 ：循环遍历所有时域样本i（0～15）
旋转因子 ：DFTSWN[id][i] 是预计算的复数，其实部为
，虚部为
（即
的分解）。
乘法累加：将输入数据DFTdatas[i]（即时间域中的第i个样本）与其对应的旋转因子相乘所得的结果依次累加到目标位置DFTResults[id]的实部与虚部中。
（3）归一化处理 ：最终结果乘以系数
（即2/16 = 0.125）。
特定归一化需求 ：为了将结果缩放到特定范围（如幅度谱分析）。
适应实信号特性：对于实信号x(n)，其DFT结果满足共轭对称性特征，在计算时只需关注前N/2+1个频率点（即0到8）。值得注意的是，在直流分量(k=0)以及奈奎斯特频点(k=8)处无需乘以二倍系数；而对于其他各分量的幅度值则应乘以二倍系数以准确反映正负频率成分所携带的能量总量。
对于长度为 N = 16 的DFT，频率分量 k 的旋转因子定义为：
其中：
k 是频率索引（如基波 k = 1，二次谐波 k = 2）
i 是时域样本索引（0～15）
j 是虚数单位。
3）基波（k = 1）计算
倘若 i = 2，则：
因此：
特殊点：i = 8
4）二次谐波（k = 2）计算
倘若：i = 3，则：
因此：
特殊点：i = 4
5）数据生成步骤
（1）确定 N 和 k：N = 16，基波 k = 1，二次谐波 k = 2。
（2）计算每个 i 的旋转因子：对 i = 0 到 15，计算
和
。
（3）数值格式化：使用科学计数法表示，保留 6位有效数字。极小值
近似为 0。
（4）预计算优化：避免实时计算复数指数，提高效率。
5、递推离散傅里叶
傅里叶变换是一种数学工具，在信号处理领域具有重要作用；积分运算可以通过迭代递推方法来优化计算过程；尤其是周期性函数的处理中更为显著；移除最后一个数据点后，在保持其他N-1个数据不变的前提下；然后依次加入新的数据点以完成整个计算流程；为了深入理解这一概念，请让我们先探讨移动平均滤波器的工作原理：
上面的这个公式可以写成迭代也就是递推的形式：
同样地，在基于 sin 和 cos 函数的周期性的基础上，DFT 通过将多项式乘法的总和转换为逐次迭代的过程得以实现。
C 代码：
x(i) = GetFromADC(); //获取AD的值
X->Real += (x(i) – x(i-N)) * cos( 2*PI* i/N * harmonic)/N;
X->Image += (x(i) – x(i-N)) * sin(-2*PI* i/N * harmonic)/N; 基于 cos 和 sin 是周期函数的事实可知, 因此可以推知, 表达式 cos(2π(i * harmonic)/N) 和 cos(2π(i * harmonic % N)/N) 会给出相同的结果. 此外, 其中 (i * Harmonic % N) 的具体取值范围是从 0 到 N-1 的整数.
总结一下：
傅里叶变换既可以非常复杂难以理解（deep），也可以非常简单易懂（shallow）。对于离散傅里叶变换（DFT）的具体公式部分，在深入研究后会发现其本质并不复杂（complicated），甚至可以说是非常直观（intuitive）。此外，在理解DFT的基本原理后（be used to calculate each harmonic component），还可以进一步掌握如何快速地完成这一过程（efficiently complete）。为了更加深入的理解快速傅里叶变换（FFT），想象一下这个算法的设计思路：它通过巧妙地分解问题来显著提高计算效率（significant improvement in computational efficiency）。
6、问题的解答
1）问题A的解答
在该函数CalculateHarmonic(Complex* X, int harmonic)中可以看出DFT的各项谐波计算相互独立不受其他次谐波的影响然而，在问题A中并不需要计算2倍频（100Hz）、3倍频（150Hz）等谐波这正是DFT的一个显著优势接下来我们将详细阐述如何定义两个复数变量的结构
typedef short int16;
typedef int int32;
typedef struct SComplex
{
int16 R; //实部
int16 I; //虚部
}Complex;
typedef struct SComplex32
{
int32 R; //实部
int32 I; //虚部
}Complex32;接着， 定义两个常数以及电压电流的结构：
#define N 16 //每周期采样点数
#define LOG2_N 4 // log2(N)
struct UI
{
Complex U; //电压的结果
Complex I; //电流的结果
int16 Voltage[N];//先前的 N 个电压
int16 Current[N];//先前的 N 个电流
Complex32 UAcc; //电压的累加器
Complex32 IAcc; //电流的累加器
int Index; //迭代索引计数器， 8-BIT MCU 可以为 char, 如果 N < 256
Complex W[N]; //N 点的 cos, sin 系数
}ui;初始化，cos、sin 系数数组应该事先计算好：
void UI_Init()
{
for (int i=0; i
{
ui.W[i].R = (int16) (::cos( 2*3.1415927*i/N) * (1<<14) + 0.5); //应离线计算！！！
ui.W[i].I = (int16) (::sin(-2*3.1415927*i/N) * (1<<14) + 0.5);
ui.Voltage[i] = 0;
ui.Current[i] = 0;
}
ui.UAcc.R = 0; ui.UAcc.I = 0;
ui.IAcc.R = 0; ui.IAcc.I = 0;
ui.Index = 0;
}下面的代码不断递推， 可以求出电压和电流的复数：
void UI_Calculate(int16 voltage, int16 current)
{
int16 d;
d = voltage - ui.Voltage[ui.Index];
ui.Voltage[ui.Index] = voltage;
ui.UAcc.R += (d * ui.W[ui.Index].R) >> (13 + LOG2_N - 16);
ui.UAcc.I += (d * ui.W[ui.Index].I) >> (13 + LOG2_N - 16);
ui.U.R = (int16) (ui.UAcc.R >> 16);
ui.U.I = (int16) (ui.UAcc.I >> 16);
d = current - ui.Current[ui.Index];
ui.Current[ui.Index] = current;
ui.IAcc.R += (d * ui.W[ui.Index].R) >> (13 + LOG2_N - 16);
ui.IAcc.I += (d * ui.W[ui.Index].I) >> (13 + LOG2_N - 16);
ui.I.R = (int16) (ui.IAcc.R >> 16);
ui.I.I = (int16) (ui.IAcc.I >> 16);
ui.Index = (ui.Index + 1) & (N-1);
}上面的计算采用定点运算技术实现DFT计算过程，并支持使用不同位宽的运算方案以适应不同的需求。为此处需将电压与电流进行归一化处理。例如，在系统中最大的电压幅度设定为Vmax=400V（伏特），而最大的电流幅度设定为I_max=20A（安培）。在数字系统中统一采用Q15格式表示数值：Q₁₅=₂¹⁵=₃₂7₆₈（十六进制）。其中V_max和I_max分别对应Q₁₅格式中的数值。因此，在演示主程序中：
8000 ---〉8000/Q15 * Vmax = 97.66V
4000 ---〉4000/Q15 * Imax = 2.441A
至于功率，很简单, 用电压乘以电流的轭(用 j 来代替复数i, 以免混淆)：
P + jQ = U*I’ = (ui.U.R + j ui.U.I) * (ui.I.R – j ui.I.I)
P是有功功率， Q是无功功率；
功率因数为：cos(theta) = P / sqrt(PP +QQ)
///
计算有功功率、无功功率两种方法:
DFT方法可以通过富里叶变换计算电流电压矢量，并从而能够计算出有功功率和无功功率；电流电压的方向性不再是问题。
U*I. 在一个完整的周期内, 即如楼主所述, 电压与电流在同一相位下的乘积即为有效功率；当电压与移相后的电流呈90度相位差时, 即当电流提前或滞后16个点的位置时, 可以计算出无功功率值。
定期采样采集V和I值，并计算它们在每个时间段内的瞬时功率乘以时间间隔之积之和作为单位时间内的电能总量。将电压序列进行相位移90度后得到新的电压序列; 同理将该新的电压序列与电流序列逐点相乘后再累加求和就可以获得无功电量。
///
Visual C++下的演示主程序 ：
#include "stdafx.h"
#include "Math.h"
#include "stdio.h"
#include "UI.h"
#define Magnitude(c) ((int) sqrtf(c.R*c.R + c.I*c.I))
#define PI 3.14159265f
int _tmain(int argc, _TCHAR* argv[])
{
int16 voltage, current;
Complex PQ;
UI_Init();
for (int i=0; i<1000; i++)
{
voltage = (int16) (::sin(2*PI*i/N)*8000); //模拟采样电压
current = (int16) (::sin(2*PI*i/N + PI/3)* 4000); //模拟采样电流
UI_Calculate(voltage, current);
}
PQ.R = (ui.U.R * ui.I.R + ui.U.I * ui.I.I) >> 15;
PQ.I = (ui.U.I * ui.I.R - ui.U.R * ui.I.I) >> 15;
printf("Voltage: %d\n", Magnitude(ui.U));
printf("Cuurent: %d\n", Magnitude(ui.I));
printf("Power Factor: %d\n", PQ.R * 1000 / Magnitude(PQ));
::getchar();
return 0;
}结果
Voltage: 8000
Cuurent: 3999
Power Factor: 500
2）问题B的解答
目前大家已掌握DFT单独实现各谐波计算的能力。对于这道题而言，使用 DFT 是可行的选择之一；而 FFT 则是更为高效的替代方案。需要注意的是，在这种情况下最低分辨率要求为5Hz；具体来说，在本题中采样频率设定在800Hz时，则周期长度精确计算得出为1.25毫秒。
5Hz 周期 T5 = 200ms。因此，5Hz 数据窗口的长度为 N = T5/T800 = 160，这样 50Hz、55Hz 就分别是 10、11次谐波。
定义常数：
#define N 160
#define LOG2_N 8计算 cos、sin 的系数，注意：**(1 <<(14 + LOG2_N))/N **的作用。
for (int i=0; i<N; i++)
{
ui.W[i].R = (int16) (::cos( 2*3.1415927*i/N) * (1<<(14 + LOG2_N)) / N + 0.5);
ui.W[i].I = (int16) (::sin(-2*3.1415927*i/N) * (1<<(14 + LOG2_N)) / N + 0.5);
ui.Voltage[i] = 0;
}计算：
void UI_Calculate(int16 voltage)
{
int16 d;
int i;
d = voltage - ui.Voltage[ui.Index];
ui.Voltage[ui.Index] = voltage;
i = (ui.Index * 10) % N;
ui.U10_Acc.R += (d * ui.W[i].R) >> (13 + LOG2_N - 16);
ui.U10_Acc.I += (d * ui.W[i].I) >> (13 + LOG2_N - 16);
ui.U10.R = (int16) (ui.U10_Acc.R >> 16);
ui.U10.I = (int16) (ui.U10_Acc.I >> 16);
i = (ui.Index * 11) % N;
ui.U11_Acc.R += (d * ui.W[i].R) >> (13 + LOG2_N - 16);
ui.U11_Acc.I += (d * ui.W[i].I) >> (13 + LOG2_N - 16);
ui.U11.R = (int16) (ui.U11_Acc.R >> 16);
ui.U11.I = (int16) (ui.U11_Acc.I >> 16);
ui.Index = (ui.Index + 1) % N;
}注意：(13 + LOG2_N - 16) 的作用。
演示主程序：
int _tmain(int argc, _TCHAR* argv[])
{
UI_Init();
float Hz50 = 2 * PI * 50 / 800;
float Hz55 = 2 * PI * 55 / 800;
for (int i=0; i<1000; i++)
{
UI_Calculate((int16) (::sin(Hz50*i)*8000 + ::sin(Hz55*i)* 4000));
}
printf("50Hz: %d\n", Magnitude(ui.U10));
printf("55Hz: %d\n", Magnitude(ui.U11));
::getchar();
return 0;
}结果：
50Hz: 8000
55Hz: 4000
这个案例的一个问题是参数设置为N=160这对于具有较小内存容量的微控制器来说可能会导致性能问题为此我们需要进行一些调整一种方法是改变采样频率另一种方法则是不更改采样率但跳过若干数据点以模拟地调整采样频率例如每隔五次采样取一个数据点这样总的数据量就会减少到原来的三分之一
#define N 32
#define LOG2_N 5
for (int i=0; i<1000; i++)
{
if ((i % 5) == 0)
UI_Calculate((int16) (::sin(Hz50*i) * 8000 + ::sin(Hz55*i) * 4000));
}得到了一样的结果，而数据 buffer 为 32， 可以计算出上到 15次谐波。
介绍完 DFT 后轮到 FFT。我意识到变速不变调不适合作为 FFT 的例子, 因为变速不变调涉及了很多其他的概念, 而基于 Matlab 编程的验证程序虽然简短, 但采用了 Matlab 中的 FFT 和 IFFT 命令, 因此缺乏代表性。相比之下, 这些 DFT/FFT 与 IDFT/IFFT 的定点 C/C++ 程序均为本人编写, 可以为它们提供更详细的讲解。
7、总结
傅里叶变换的本质是将一个信号通过正交分解（e^(jωt) = cos(ωt) + j sin(ωt)）转化为无限多个正弦波的叠加形式；这些无限多个正弦波经过重新组合又能够精确重构出原始信号。与白光透过三棱镜分离出不同色光并可借由同一装置恢复为白光的过程相似地傅里叶变换相当于一种能够分离与重组信号频率成分的独特工具
e^(jωt) = cos(ωt) + j sin(ωt)
类似于钟表的指针以角速度 ω 转动时，在横轴方向上的投射结果为 sin(ωt) 。设想两个不同时间点叠加在一起的钟表盘面，在你的观察视角下（假设你置身于其中一个秒针时），另一块表的秒针呈现静止状态，并且其在横向轴上的投射结果为固定数值。
设想一下多台运转中的计时器（手表）其秒针运行的速度各不相同 而你所乘用的手表中的秒针则可调节其转速 那么你将发现 无论怎样调节自己的计时器（手表） 都可在某一特定时间点让某一秒钟计时标记看似静止不动 其在你的视角下的显示位置保持恒定不变 这种现象与转速设置无关
傅里叶变换的结果是什么？以离散傅里叶变换（DFT）或快速傅里叶变换（FFT）为例说明：对N点的数字信号进行快速傅里叶变换（FFT），计算后得到N/2个复数结果。每个复数结果实际上对应一个正弦波分量。假设采样频率为F，则基频计算公式为ω₀ = 2πF/N。这些N/2个复数结果：
当k=0时,Y[k]=x₀+jy₀,其频率分量ω= 0 ，这相当于直流分量。
当k=1时,Y[k]=x₁+jy₁=r₁e^{jθ₁} ，其频率分量ω= ω₀ ，对应一个振幅为r₁、相位为θ₁的正弦波。
同样地，
当k=2时,Y[k]=x₂+jy₂=r₂e^{jθ₂} ，其频率分量ω=2ω₀ ，对应一个振幅为r₂、相位为θ₂的正弦波。
以此类推，
对于一般的k值，
有Y[k]=x_k +j y_k = r_k e^{jθ_k} ，其对应的频率分量是ω=kω_{
⁰} ，这也是一幅度为r_k、相位角为θ_k的正弦波。
特别地，
当k=N/2 -1时,
有
Y[N/2 - 1] =
由此可见，在这里所指的一些复数的意义上来说它们是代表正弦波的原因而不是一般意义上的复数也就是说将这些正弦波结合起来就可以恢复原始形态
首先, 我发布的DFT程序均为本人编写, 并且提供了汇编语言版本. C版本则主要采用了16位整数乘法和32位加法运算, 同时配合少量移位操作完成. 对于除法操作（尤其是使用模运算%而非2的幂次方的情况）, 本方案能够有效规避. 在每次定时中断触发时进行采样并计算, 这一过程由于采用了递推算法, 计算量相对较低（仅需执行两个16位乘法指令和两个32位加法指令）. 然而, 为了防止浮点运算带来的潜在精度问题, 必须采用定点缩放技术进行处理, 这一做法相较于直接采用浮点运算具有一定的复杂性. 尽管演示程序为了通用性采用了PC平台作为开发环境, 但DFT算法程序在AVR、51等8位微控制器上的实现完全没有问题. 不要再提出单片机无法实现的观点
对于FFT程序来说，在汇编语言中有相应的实现。其中C/C++版本主要采用16位整数乘法、32位加法以及少量移位操作的方式进行计算运算，并且其运算效率较高。然而，在运行于8位MCU设备中时可能会遇到性能瓶颈。这是因为当数据窗口点数较小时，在这种环境下使用可能会遇到性能瓶颈问题。此时DFT方法更为合适；而当数据窗口点数较大时，则会导致计算速度跟不上需求。因此我就不再深入介绍了。
最后，请展示我通过 Matlab 编程实现的速度恒定且无 distortion 的算法验证程序，并以此收尾。简单阐述其原理：
该系统采用了短时Fourier变换（Short Time Fourier Transform）进行处理，并运用了Hamming window作为窗口函数。
The short-time Fourier transform (STFT) is employed here, where the time segments are defined as segment = N/4. The data is partitioned into consecutive blocks: the first block spans from 0 to N-1, the second from N/4 to N+N/4-1, the third from N/2 to N+N/2-1, and so on.
2）做 FFT，计算出幅度和相位。
3）计算新的幅度和相位。幅度通过插植，相位得把 ωt 计入：da(2: (1 + NX/2)) = (2pisegment) ./ (NX ./ (1: (NX/2)));
4）用新的幅度和相位产生新的复数，加窗并作 IFFT 生成变速后的音频数据。
SPEED = 2;
[in_rl, fs] = wavread('C:\windows\Media\Windows XP Startup.wav');
in = in_rl(:, 1)';
sizeOfData = length(in);
7. N = 2048;
segment = N/4;
window = hamming(N)';
X = zeros((1+ N/2), 1 + fix((sizeOfData - N)/segment));
c = 1;
for i = 0: segment: (sizeOfData - N)
fftx = fft(window .* in((i+1): (i+N)));
X(:, c) = fftx(1: (1+N/2))';
c = c + 1;
end;
18. [Xrows, Xcols] = size(X);
NX = 2 * (Xrows - 1);
Y = zeros(Xrows, round((Xcols - 1) / SPEED));
22. da = zeros(1, NX/2+1);
da(2: (1 + NX/2)) = (2*pi*segment) ./ (NX ./ (1: (NX/2)));
phase = angle(X(:, 1));
c = 1;
for i = 0: SPEED: (Xcols-2)
X1 = X(:, 1 + floor(i));
X2 = X(:, 2 + floor(i));
df = i - floor(i);
magnitude = (1-df) * abs(X1) + df * (abs(X2));
dangle = angle(X2) - angle(X1);
dangle = dangle - da' - 2 * pi * round(dangle/(2*pi));
Y(:, c) = magnitude .* exp(j * phase);
phase = phase + da' + dangle;
c = c + 1;
end
[Yrows, Ycols] = size(Y);
out = zeros(1, N + (Ycols - 1) * segment );
c = 1;
for i = 0: segment: (segment * (Ycols-1))
Yc = Y(:, c)';
Ynew = [Yc, conj(Yc((N/2): -1: 2))];
out((i+1): (i+N)) = out((i+1)i+N)) + real(ifft(Ynew)) .* window;
c = c + 1;
end;
48. wavplay(out, fs);已经完成了FFT相关内容的讲解。关于FFT的内容较为基础且易于理解，在此不再赘述细节。如果感兴趣的话，请参考附带的定点C++实现及其逆变换程序作为参考代码。
class CComplex
{
public:
short R;
short I;
CComplex() {}
CComplex(short r, short i) { R = r; I = i; }
void Set(short r, short i) { R = r; I = i; }
CComplex(CComplex& complex) { R = complex.R; I = complex.I; }
inline void operator=(CComplex& complex) { R = complex.R; I = complex.I; }
void operator+=(CComplex& complex)
{
R += complex.R;
I += complex.I;
}
void operator-=(CComplex& complex)
{
R -= complex.R;
I -= complex.I;
}
inline void Add(CComplex& c1, CComplex& c2)
{
R = c1.R + c2.R;
I = c1.I + c2.I;
}
inline void Sub(CComplex& c1, CComplex& c2)
{
R = c1.R - c2.R;
I = c1.I - c2.I;
}
inline void Add(CComplex& c1, CComplex& c2, int shift)
{
R = (c1.R + c2.R) >> shift;
I = (c1.I + c2.I) >> shift;
}
inline void Sub(CComplex& c1, CComplex& c2, int shift)
{
R = (c1.R - c2.R) >> shift;
I = (c1.I - c2.I) >> shift;
}
inline void Mul(CComplex& c1, CComplex& c2, int shift)
{
R = (c1.R * c2.R + c1.I * c2.I) >> shift;
I = (c1.I * c2.R - c1.R * c2.I) >> shift;
}
inline void MulConj(CComplex& c1, CComplex& c2, int shift)
{
R = (c1.R * c2.R - c1.I * c2.I) >> shift;
I = (c1.I * c2.R + c1.R * c2.I) >> shift;
}
inline void Shift(int n) { R >>= n; I >>= n; }
};在嵌入式，构造函数中的 BitReversal, W 数组应该预先离线计算
class CFFT
{
protected:
static const int N = 32;
CComplex W[N];
int BitReversal[N];
public:
CComplex Y[N];
protected:
void ReverseBit()
{
int rev = 0;
int t = N/2;
int mask;
for (int i = 0; i < N-1; i++) {
BitReversal = rev;
mask = t;
while (rev >= mask) {
rev -= mask;
mask >>= 1;
}
rev += mask;
}
BitReversal[N-1] = N-1;
}
public:
CFFT()
{
for (int i=0; i<N; i++) {
W.R = (short) ( ::cos(PI * 2*i / N) * 16384 + 0.5);
W.I = (short) (-::sin(PI * 2*i / N) * 16384 + 0.5);
}
ReverseBit();
}
void Calculate(short* x)
{
short t;
int i, j, n, index1, index2, indexW, di;
CComplex result;
for(i=0; i<N; i++) {
n = BitReversal;
if(n > i) {
t = x;
x = x[n];
x[n] = t;
}
}
for (i=0; i<N; i+=2) {
t = x[i+1];
Y.Set(x + t, 0);
Y[i+1].Set(x - t, 0);
}
for (di=N>>2, n=2; n<N; n<<=1, di >>= 1) {
for (j=0; j<N; j+=(n<<1)) {
indexW = 0;
for (i=0; i<n; i++) {
index1 = i + j;
index2 = index1 + n;
result.Mul(Y[index2], W[indexW], 14);
Y[index2].Sub(Y[index1], result, 1);
Y[index1].Add(Y[index1], result, 1);
indexW += di;
}
}
}
}
void Inverse()
{
int i, j, n, index1, index2, indexW, di;
CComplex result;
for(i=0; i<N; i++) {
n = BitReversal;
if(n > i) {
result = Y;
Y = Y[n];
Y[n] = result;
}
}
for (i=0; i<N; i+=2) {
result = Y[i+1];
Y[i+1].Sub(Y, result, 1);
Y.Add(Y, result, 1);
}
for (di=N>>2, n=2; n<N; n<<=1, di >>= 1) {
for (j=0; j<N; j+=(n<<1)) {
indexW = 0;
for (i=0; i<n; i++) {
index1 = i + j;
index2 = index1 + n;
result.MulConj(Y[index2], W[indexW], 14);
Y[index2].Sub(Y[index1], result);
Y[index1].Add(Y[index1], result);
indexW += di;
}
}
}
}
};本内容源自21ic的21ic_highgear ， fellowes 老铁们觉得内容很有参考价值吗？点赞支持一下哦！