授之以渔：几个例子弄清复立叶变换的应用

只看该作者 · 2010-8-30 21:44

继续听课~~~

只看该作者 · 2010-8-30 22:28

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[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

只看该作者 · 2010-8-30 22:44

上面的例子有一个问题就是　Ｎ＝１６０，这对于小ｒａｍ容量的　ｍｃｕ　来说，　不太合适。我们可以做一些改动。一是改变采样频率，　二是保持采样频率不变，跳过几个点，变相的改变采样频率。这里我们可以每采样５次，计算一次，　这样　Ｎ　＝１６０／５＝３２．

#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 次谐波。

　

只看该作者 · 2010-8-30 22:57

介绍完 DFT, 下面轮到FFT. 我现在发现变速不变调不适合作为 FFT 的例子，因为变速不变调涉及了很多其他的概念，验证程序用 matlab 做的，虽然不长，但是用了 matlab 里 FFT, IFFT 的命令，所以没有典型性。而DFT/FFT与逆变换 IDFT/IFFT 的定点 c/c++ 程序都是我自己做的, 可以更详细的讲解。

FFT 容我这两天设想一个经典的例子，另外开一个帖子讲解。

只看该作者 · 2010-8-30 23:01

再顶！

期待highgear老师讲课。

只看该作者 · 2010-8-31 11:02

公式都不记得了，再学习一遍

只看该作者 · 2010-8-31 11:54

来听课的

只看该作者 · 2010-8-31 14:04

顶楼主

只看该作者 · 2010-8-31 14:34

期待...

只看该作者 · 2010-8-31 15:32

顶一下highgear大师！讲的太经典了

只看该作者 · 2010-8-31 15:58

期待。。。

只看该作者 · 2010-8-31 20:17

非路过

只看该作者 · 2010-8-31 20:59

顶之

只看该作者 · 2010-8-31 22:57

yeah ,get a chair ,and listen to the couch carefully in the chair!

只看该作者 · 2010-8-31 23:35

本帖最后由 highgear 于 2010-8-31 23:42 编辑

总结：

傅里叶变换的实质是把一个信号通过正交分解（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] =

所以，这些复数的意义在于正弦波的代表，不是一般意义上的复数。把上面的正弦波叠加在一起，又可以得到原来的波形。