A) 系统用 adc (小于 16-bit) 采样 50Hz 交流电流电压， 采样频率800hz, 试求出电流电压幅值以及功率和功率因数。
B) 上面的50hz 电压中， 混入了另一个 55hz 的电压， 求出这两个电压的幅值。
这两道题使用 16-bit, 32-bit 的整数运算， 不使用浮点运算， 可以在 mcu 上实现。
C) 完成一个 wav 声音文件的变速不变调的程序。
（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/%E ... 6%E5%8F%98%E6%8D%A2
(3) 离散傅里叶变换(DFT)http://zh.wikipedia.org/zh-cn/%E ... 6%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.
总结一下：
傅里叶变换可以很深奥， 也可以很浅显。对于离散的傅里叶变换的公式， 只要认真的看看很容易看明白， 更何况还有代码说明。通过理解 dft 如何计算出某一个谐波， 就可以进一步计算出所有谐波， 再想象一下， 某一个算法， 可以快速的计算出所有的谐波， 这样， 就可以很容易的理解 fft.
(5) 问题 A 的解答
在上面的代码 CalculateHarmonic(Complex* X, int harmonic) 中可以看出dft 的各次谐波计算是独立的， 不依赖其它次谐波。而且，问题 A 不需要计算2次(100hz)，3次(150hz)等等谐波，这是 dft 的优点之一。首先，定义两个复数的结构：
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<N; i++) {
ui.W.R = (int16) (::cos( 2*3.1415927*i/N) * (1<<14) + 0.5); //应离线计算！！！
ui.W.I = (int16) (::sin(-2*3.1415927*i/N) * (1<<14) + 0.5);
ui.Voltage = 0;
ui.Current = 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计算使用的是 16-bit, 32-bit 的定点运算，这里需要把电压和电流单位化。比如系统最大电压幅值为 Vmax = 400V，最大电流幅值 Imax = 20A, 在数字系统中统一归一化: Q15 = 2^15 = 32768. 即 Vmax，Imax在数字系统对应Q15 = 32768. 因此，演示主程序中的：
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(P*P +Q*Q)
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
6） 问题 B 的解答
现在大家已经知道了， DFT 可以单独的计算各个谐波。这道题，同样可以用 DFT 来做， 当然也可以用 FFT 来做。 50Hz与 55hz 相差 5Hz, 所以必须采用 5Hz 的分辨率。采样频率为800Hz,
周期T800 = 1.25ms;
5Hz周期T5 = 200 ms. 因此，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.R = (int16) (::cos( 2*3.1415927*i/N) * (1<<(14 + LOG2_N)) / N + 0.5);
ui.W.I = (int16) (::sin(-2*3.1415927*i/N) * (1<<(14 + LOG2_N)) / N + 0.5);
ui.Voltage = 0;
}
{
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.R) >> (13 + LOG2_N - 16);
ui.U10_Acc.I += (d * ui.W.I) >> (13 + LOG2_N - 16);
ui.U10.R = (int16) (ui.U10_Acc.R >> 16);
ui.U10.I = (int16) (ui.U10_Acc.I >> 16);
ui.U11_Acc.R += (d * ui.W.R) >> (13 + LOG2_N - 16);
ui.U11_Acc.I += (d * ui.W.I) >> (13 + LOG2_N - 16);
ui.U11.R = (int16) (ui.U11_Acc.R >> 16);
ui.U11.I = (int16) (ui.U11_Acc.I >> 16);
}
注意 (13 + LOG2_N - 16) 的作用。
演示主程序：
{
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
#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 次谐波。
FFT 容我这两天设想一个经典的例子， 另外开一个帖子讲解。
总结：
傅里叶变换的实质是把一个信号通过正交分解（e^(jωt) = cos(ωt) + j sin(ωt) ), 分解成无数的正弦信号， 而这些无数的正弦信号还可以重新被合成为原来的信号。就像白光通过三棱镜分解成光谱， 再通过三棱镜可以被还原成白光一样, 傅里叶变换就是那个三棱镜， 或者说三棱镜就是一个傅里叶变换。
e^(jωt) = cos(ωt) + j sin(ωt)
可以看做钟表的指针以的角速度 ω 旋转时， 指针在纵横两个方向上的投影， 在横轴上的投影就是 sin(ωt) . 假设两个不同时间的钟表叠放在一起， 你坐在其中的一个秒指针上， 你会发现另一块表的秒指针是静止的， 并且在你的指针上的投影是固定的。现在设想一下很多块表的秒指针以不同的速率旋转， 而你所乘坐的秒指针可以控制旋转速率， 那么你会发现， 总可以使某一个秒指针看上去是静止的， 即在你的指针上的投影是常数，与速度无关。
傅里叶变换出来的是什么? 以离散的傅里叶变换ＤＦＴ／ＦＦＴ 来说明，对Ｎ点的数据做傅里叶变换，得到了 N/2 个复数， 这每一个复数实际上代表了一个正弦波, 假设 采样频率为 F, 那么基本频率为 ω0 = 2*PI*F/N
这 N/2 个复数：
Y[0] = x0 + j y0 ： ω = 0, 即 DC.
Y[1] = x1 + j y1 = r1* e^(j a1) : ω = ω0, 代表正弦波 r1* sin( ω0 * t + a1)
Y[2] = x2 + j y2 = r2* e^(j a2) : ω = 2* ω0, 代表正弦波 r2* sin(2* ω0 * t + a2)
....
Y[k] = xk + j yk = rk* e^(j ak) : ω = k* ω0, 代表正弦波 rk* sin(k* ω0 * t + ak)
...
Y[N/2 - 1] =
所以， 这些复数的意义在于正弦波的代表， 不是一般意义上的复数。把上面的正弦波叠加在一起， 又可以得到原来的波形。
首先， 我贴出的DFT程序都是我自己写的， 而且有汇编的版本。 大家已经看到了， c 版本完全使用了16-bit 整数乘法， 32-bit 加法以及少量的移位操作， 除法（主要是 %，用于非 2的次方的 N) 可以完全避开。 可以设定在每次定时中断里采样后计算， 由于递推， 计算量很低（2个16bit乘法， 2个32bit 加法）。 唯一的问题是， 必须使用定点 scale 转换以避免浮点运算， 这不如直接使用浮点直观， 对没有处理经验的程序员可能是一个挑战。虽然演示程序为了通用起见用了 PC, 但是dft 算法程序用在 avr, 51等 8-bit mcu 是完完全全没有问题的。不要再说出单片机不能实现的胡话出来。
FFT程序我也有汇编的版本， C/C++ 版本也是采用了16-bit 整数乘法， 32-bit 加法以及少量的移位操作， 效率很高， 不过在 8-bit mcu 上可能用不上， 因为， 数据窗口点数少了， 用 dft 更好， 数据窗口点数多了， 8-bit mcu 上太慢， 不实用。 因此我就不介绍了。
最后， 我贴出我用 matlab 做的变速不变调的算法验证程序， 作为结束。
简单的讲一讲原理：
下面的程序使用了短时博立叶变换（short time fourier transform), 窗口函数为 hamming。
1） 短时博立叶变换, 这里的片断 segment = N/4, 数据被分割为 0 到 N-1, N/4 to N+N/4-1, N/2 to N+N/2-1， 依次类推。
2） 做 fft, 计算出幅度和相位。
3）计算新的幅度和相位。幅度通过插植， 相位得把 wt 计入： da(2: (1 + NX/2)) = (2*pi*segment) ./ (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);
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;
[Xrows, Xcols] = size(X);
NX = 2 * (Xrows - 1);
Y = zeros(Xrows, round((Xcols - 1) / SPEED));
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;
wavplay(out, fs);
/3 
