復(fù)立葉變換你知道是什么嗎？

下面兩道題關(guān)于使用復(fù)利葉變換的， 這應(yīng)該是很常見的嵌入式問題：A)系統(tǒng)用 adc (小于 16-bit) 采樣 50Hz 交流電流電壓， 采樣頻率800hz, 試求出電流電壓幅值以及功率和功率因數(shù)。B)上面的50hz 電壓中， 混入了另一個(gè) 55hz 的電壓， 求出這兩個(gè)電壓的幅值。這兩道題使用 16-bit, 32-bit 的整數(shù)運(yùn)算， 不使用浮點(diǎn)運(yùn)算， 可以在 mcu 上實(shí)現(xiàn)。C)完成一個(gè) wav 聲音文件的變速不變調(diào)的程序。

（1）復(fù)數(shù)的基礎(chǔ)知識(shí)

在講解 fourier transform 前， 大家必須知道一點(diǎn)基本的復(fù)數(shù)知識(shí)。在復(fù)平面上的一個(gè)點(diǎn) P (x, y) 用復(fù)數(shù)表示為: P = x + i y用極坐標(biāo)表示為： P = r * e^(i a)這里,r = sqrt(x*x + y*) 是點(diǎn) (x, y) 到原點(diǎn)的距離， a = arctan2(x, y) 是角度, e 是自然常數(shù)。這里引出了一個(gè)非常重要的表達(dá)式： e^(i a)= cos(a) + i sin(a)這個(gè)表達(dá)式，是利用復(fù)數(shù)完成角度變換和三角函數(shù)變換的利器。例如，把點(diǎn) P 旋轉(zhuǎn) b 角度，那么新點(diǎn)(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)傅里葉變換的基礎(chǔ)知識(shí)傅里葉變換是一個(gè)積分變換， 公式就不提供了， 有興趣的同學(xué)可以直接訪問下面的連接， 以獲得更詳盡的解釋：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 次諧波的復(fù)數(shù)；N 為周期采樣點(diǎn)數(shù)；x(n)為輸入，n從0 到N-1； 用偽代碼更直觀地說明： void CalculateHarmonic(Complex* X, int harmonic) { for (int i=0; iReal= x(i) * cos( 2*PI* i/N * harmonic) / N; X->Image = x(i) * sin(-2*PI* i/N * harmonic) / N; } }可以看到，離散傅里葉變換基本運(yùn)算其實(shí)很簡(jiǎn)單， 沒有那么復(fù)雜。只要有了 N 個(gè)輸入，比如說通過AD 采樣了 N 個(gè)數(shù)據(jù)后，可以輕易的計(jì)算出各個(gè)諧波，雖然計(jì)算量大了些。下面要做的就是減少計(jì)算量，這可以用兩種方法， 一種當(dāng)然就是熟知的 FFT, 還有一種就是遞推。(3)遞推離散傅里葉 (Recursive DFT)傅里葉變換是一個(gè)積分變換，積分當(dāng)然可以使用迭代遞推來減少運(yùn)算，尤其是周期性的函數(shù)。只要把最后一個(gè)數(shù)據(jù)仍出去，保持其他 N-1 個(gè)數(shù)據(jù)不變，加入一個(gè)新的數(shù)據(jù)就可以了。為了理解這一點(diǎn)，先考慮一下移動(dòng)平均濾波算法： Y(k-1) = (x(k-1) + x(k-2) + … + x(k-N)) /N上面的這個(gè)公式可以寫成迭代也就是遞推的形式： Y(k) = Y(k-1) + (x(k) – x(k-N)) /N同理，由于sin, cosin函數(shù)的周期性，dft 可以由多項(xiàng)式乘法和的形式變換成迭代遞推的形式： 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) / NC 代碼： 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 是周期函數(shù)，所以 cos(2*PI* (i * harmonic) / N) 與cos(2*PI* (i * harmonic % N) / N) 是一樣的，(i * harmonic % N) 的取值范圍：0 to N-1.

總結(jié)一下：傅里葉變換可以很深?yuàn)W， 也可以很淺顯。對(duì)于離散的傅里葉變換的公式， 只要認(rèn)真的看看很容易看明白， 更何況還有代碼說明。通過理解 dft 如何計(jì)算出某一個(gè)諧波， 就可以進(jìn)一步計(jì)算出所有諧波， 再想象一下， 某一個(gè)算法， 可以快速的計(jì)算出所有的諧波， 這樣， 就可以很容易的理解 fft.

(5) 問題 A 的解答在上面的代碼CalculateHarmonic(Complex* X, int harmonic)中可以看出dft 的各次諧波計(jì)算是獨(dú)立的， 不依賴其它次諧波。而且，問題 A 不需要計(jì)算2次(100hz)，3次(150hz)等等諧波，這是 dft 的優(yōu)點(diǎn)之一。首先，定義兩個(gè)復(fù)數(shù)的結(jié)構(gòu)：

typedef short int16;

typedef int int32;

typedef struct SComplex

{

int16 R;

int16 I;

} Complex;

typedef struct SComplex32

{

int32 R;

int32 I;

} Complex32;

接著， 定義兩個(gè)常數(shù)以及電壓電流的結(jié)構(gòu)：

#define N 16 //每周期采樣點(diǎn)數(shù)

#define LOG2_N 4 // log2(N)

struct UI {

ComplexU; //電壓的結(jié)果

ComplexI; //電流的結(jié)果

int16Voltage[N]; //先前的 N 個(gè)電壓

int16Current[N]; //先前的 N 個(gè)電流

Complex32 UAcc;//電壓的累加器

Complex32 IAcc; //電流的累加器

int Index;//迭代索引計(jì)數(shù)器， 8-BIT MCU 可以為 char, 如果 N < 256.

ComplexW[N];//N 點(diǎn)的 cos, sin 系數(shù)

} ui;

初始化，cos, sin 系數(shù)數(shù)組應(yīng)該事先計(jì)算好：

void UI_Init()

{

for (int i=0; i

ui.W[i].R = (int16) (::cos( 2*3.1415927*i/N) * (1<<14) + 0.5); //應(yīng)離線計(jì)算?。?！

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;

}

下面的代碼不斷遞推， 可以求出電壓和電流的復(fù)數(shù)：

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);

}

上面的計(jì)算dft計(jì)算使用的是 16-bit, 32-bit 的定點(diǎn)運(yùn)算，這里需要把電壓和電流單位化。比如系統(tǒng)最大電壓幅值為 Vmax = 400V，最大電流幅值 Imax = 20A, 在數(shù)字系統(tǒng)中統(tǒng)一歸一化:Q15 = 2^15 = 32768. 即 Vmax，Imax在數(shù)字系統(tǒng)對(duì)應(yīng)Q15 = 32768. 因此，演示主程序中的： 8000 ---〉8000/Q15 * Vmax=97.66V 4000 ---〉4000/Q15 * Imax= 2.441A至于功率，很簡(jiǎn)單, 用電壓乘以電流的軛(用 j 來代替復(fù)數(shù)i, 以免混淆)： P + jQ = U*I’ = (ui.U.R + j ui.U.I) * (ui.I.R – j ui.I.I)P是有功功率， Q是無功功率；功率因數(shù)為：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;

}

結(jié)果

Voltage: 8000Cuurent: 3999Power Factor: 500

6）問題B的解答現(xiàn)在大家已經(jīng)知道了，DFT可以單獨(dú)的計(jì)算各個(gè)諧波。這道題，同樣可以用DFT來做，當(dāng)然也可以用FFT來做。50Hz與55hz相差5Hz,所以必須采用5Hz的分辨率。采樣頻率為800Hz,周期T800 = 1.25ms;5Hz周期T5 = 200 ms.因此，5hz數(shù)據(jù)窗口的長度為N = T5 / T800 = 160，這樣50Hz, 55Hz就分別是10，11次諧波。定義常數(shù)：

#define N 160

#define LOG2_N 8

計(jì)算cos, sin系數(shù)。注意(1<<(14 + LOG2_N)) / N　的作用

for (int i=0; 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;

}

計(jì)算：

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;

}

結(jié)果：50Hz: 800055Hz: 4000

上面的例子有一個(gè)問題就是Ｎ＝１６０，這對(duì)于?。颍幔砣萘康模恚悖鮼碚f，不太合適。我們可以做一些改動(dòng)。一是改變采樣頻率，二是保持采樣頻率不變，跳過幾個(gè)點(diǎn)，變相的改變采樣頻率。這里我們可以每采樣５次，計(jì)算一次，這樣Ｎ＝１６０／５＝３２．#define N 32#define LOG2_N 5for (int i=0; i<1000; i++) {??　　if ((i % 5) == 0)?? ?　　　　UI_Calculate((int16) (::sin(Hz50*i)*8000 + ::sin(Hz55*i)* 4000));　}得到了一樣的結(jié)果，而數(shù)據(jù)buffer 為 32， 可以計(jì)算出上到 15 次諧波。

介紹完 DFT, 下面輪到FFT. 我現(xiàn)在發(fā)現(xiàn)變速不變調(diào)不適合作為 FFT 的例子， 因?yàn)樽兯俨蛔冋{(diào)涉及了很多其他的概念， 驗(yàn)證程序用 matlab 做的， 雖然不長， 但是用了 matlab 里 FFT, IFFT 的命令， 所以沒有典型性。而DFT/FFT與逆變換 IDFT/IFFT 的定點(diǎn) c/c++ 程序都是我自己做的, 可以更詳細(xì)的講解。FFT 容我這兩天設(shè)想一個(gè)經(jīng)典的例子， 另外開一個(gè)帖子講解。

總結(jié)：傅里葉變換的實(shí)質(zhì)是把一個(gè)信號(hào)通過正交分解（e^(jωt) =cos(ωt)+ j sin(ωt)), 分解成無數(shù)的正弦信號(hào)， 而這些無數(shù)的正弦信號(hào)還可以重新被合成為原來的信號(hào)。就像白光通過三棱鏡分解成光譜， 再通過三棱鏡可以被還原成白光一樣,傅里葉變換就是那個(gè)三棱鏡， 或者說三棱鏡就是一個(gè)傅里葉變換。 e^(jωt) =cos(ωt)+ j sin(ωt)可以看做鐘表的指針以的角速度ω旋轉(zhuǎn)時(shí)， 指針在縱橫兩個(gè)方向上的投影， 在橫軸上的投影就是sin(ωt). 假設(shè)兩個(gè)不同時(shí)間的鐘表疊放在一起， 你坐在其中的一個(gè)秒指針上， 你會(huì)發(fā)現(xiàn)另一塊表的秒指針是靜止的， 并且在你的指針上的投影是固定的。現(xiàn)在設(shè)想一下很多塊表的秒指針以不同的速率旋轉(zhuǎn)， 而你所乘坐的秒指針可以控制旋轉(zhuǎn)速率， 那么你會(huì)發(fā)現(xiàn)， 總可以使某一個(gè)秒指針看上去是靜止的， 即在你的指針上的投影是常數(shù)，與速度無關(guān)。傅里葉變換出來的是什么? 以離散的傅里葉變換ＤＦＴ／ＦＦＴ來說明，對(duì)Ｎ點(diǎn)的數(shù)據(jù)做傅里葉變換，得到了 N/2 個(gè)復(fù)數(shù)， 這每一個(gè)復(fù)數(shù)實(shí)際上代表了一個(gè)正弦波, 假設(shè) 采樣頻率為 F, 那么基本頻率為 ω0 = 2*PI*F/N這 N/2 個(gè)復(fù)數(shù)： 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] =所以， 這些復(fù)數(shù)的意義在于正弦波的代表， 不是一般意義上的復(fù)數(shù)。把上面的正弦波疊加在一起， 又可以得到原來的波形。

首先， 我貼出的DFT程序都是我自己寫的， 而且有匯編的版本。 大家已經(jīng)看到了， c 版本完全使用了16-bit 整數(shù)乘法， 32-bit 加法以及少量的移位操作，除法（主要是 %，用于非 2的次方的 N) 可以完全避開?？梢栽O(shè)定在每次定時(shí)中斷里采樣后計(jì)算， 由于遞推， 計(jì)算量很低（2個(gè)16bit乘法， 2個(gè)32bit 加法）。 唯一的問題是， 必須使用定點(diǎn) scale 轉(zhuǎn)換以避免浮點(diǎn)運(yùn)算， 這不如直接使用浮點(diǎn)直觀， 對(duì)沒有處理經(jīng)驗(yàn)的程序員可能是一個(gè)挑戰(zhàn)。雖然演示程序?yàn)榱送ㄓ闷鹨娪昧?PC, 但是dft 算法程序用在 avr, 51等 8-bit mcu 是完完全全沒有問題的。不要再說出單片機(jī)不能實(shí)現(xiàn)的胡話出來。FFT程序我也有匯編的版本，C/C++ 版本也是采用了16-bit 整數(shù)乘法， 32-bit 加法以及少量的移位操作，效率很高， 不過在8-bitmcu 上可能用不上， 因?yàn)椋?數(shù)據(jù)窗口點(diǎn)數(shù)少了， 用 dft 更好，數(shù)據(jù)窗口點(diǎn)數(shù)多了， 8-bit mcu 上太慢， 不實(shí)用。 因此我就不介紹了。最后， 我貼出我用matlab做的變速不變調(diào)的算法驗(yàn)證程序， 作為結(jié)束。簡(jiǎn)單的講一講原理：下面的程序使用了短時(shí)博立葉變換（short time fourier transform),窗口函數(shù)為 hamming。1） 短時(shí)博立葉變換, 這里的片斷 segment = N/4, 數(shù)據(jù)被分割為 0 到N-1,N/4 to N+N/4-1, N/2 to N+N/2-1， 依次類推。2） 做 fft, 計(jì)算出幅度和相位。3）計(jì)算新的幅度和相位。幅度通過插植， 相位得把wt 計(jì)入：da(2: (1 + NX/2)) = (2*pi*segment) ./ (NX ./ (1: (NX/2)));4）用新的幅度和相位產(chǎn)生新的復(fù)數(shù)， 加窗并作 ifft 生成變速后的音頻數(shù)據(jù)。

SPEED = 2;

[in_rl, fs] = wavread('C:\windows\Media\Windows XP Startup.wav');

in = in_rl(:, 1)';

sizeOfData = length(in);

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;

[Xrows, Xcols] = size(X);

NX = 2 * (Xrows - 1);

Y = zeros(Xrows, round((Xcols - 1) / SPEED));

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;

wavplay(out, fs);