傅里叶变换 c语言

（1） 复数的基础知识

在讲解 fourier transform 前， 大家必须知道一点基本的复数知识。

在复平面上的一个点 P (x, y) 用复数表示为: P = x + i y 用极坐标表示为： P = r * e^(i a)

这里, r = sqrt(x*x + y*) 是点 (x, y) 到原点的距离， a = arctan2(x, y) 是角度, e 是自然常数。这里引出了一个非常重要的表达式：

e^(i a) = cos(a) + i sin(a)

这个表达式，是利用复数完成角度变换和三角函数变换的利器。例如，把点 P 旋转 b 角度，那么新点(x1, y1) 的角度为 a+b, 距离仍为 r. 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))

(2) 傅里叶变换的基础知识

傅里叶变换是一个积分变换， 公式就不提供了， 有兴趣的同学可以直接访问下面的连接， 以获得更详尽的解释：

http://zh.wikipedia.org/zh-cn/???é??????????￠

(3) 离散傅里叶变换(DFT)

http://zh.wikipedia.org/zh-cn/?|???￡???é??????????￠

离散傅里叶变换的公式：

X(k) = ∑ x(n) * e^(i -2*PI* n/N * k) / N

这里 X(k) 是第 k 次谐波的复数；N 为周期采样点数；x(n)为输入，n从0 到N-1；

用伪代码更直观地说明：

void CalculateHarmonic(Complex* X, int harmonic) {

for (int i=0; i

X->Real = x(i) * cos( 2*PI* i/N * harmonic) / N; X->Image = x(i) * sin(-2*PI* i/N * harmonic) / N; } }

可以看到，离散傅里叶变换基本运算其实很简单， 没有那么复杂。只要有了 N 个输入，比如说通过AD 采样了 N 个数据后，可以轻易的计算出各个谐波，虽然计算量大了些。下面要做的就是减少计算量，这可以用两种方法， 一种当然就是熟知的 FFT, 还有一种就是递推。

(3) 递推离散傅里叶 (Recursive DFT) 傅里叶变换是一个积分变换，积分当然可以使用迭代递推来减少运算，尤其是周期性的函数。只要把最后一个数据仍出去，保持其他 N-1 个数据不变，加入一个新的数据就可以了。为了理解这一点，先考虑一下移动平均滤波算法： Y(k-1) = (x(k-1) + x(k-2) + … + x(k-N)) /N 上面的这个公式可以写成迭代也就是递推的形式： Y(k) = Y(k-1) + (x(k) – x(k-N)) /N

同理，由于sin, cosin函数的周期性，dft 可以由多项式乘法和的形式变换成迭代递推的形式：

Y(k) = Y(k-1) + x(k) * e^(i -2*PI* k /N * harmonic) / N - x(k - N) * e^(i -2*PI* (k–N) /N * harmonic) / N = Y(k-1) + (x(k) - x(k- N)) * e^(i -2*PI* k /N * harmonic) / N

C 代码：

x(i) = GetFromADC();

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*PI* (i * harmonic) / N) 与cos(2*PI* (i * harmonic % N) / N) 是一样的，(i * harmonic % N) 的取值范围：0 to N-1.