本帖最后由 highgear 于 2010-8-30 08:33 编辑
（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/%E5%82%85%E9%87%8C%E5%8F%B6%E5%8F%98%E6%8D%A2
(3) 离散傅里叶变换(DFT)
http://zh.wikipedia.org/zh-cn/%E7%A6%BB%E6%95%A3%E5%82%85%E9%87%8C%E5%8F%B6%E5%8F%98%E6%8D%A2
离散傅里叶变换的公式：
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<N; 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. |
楼主
