【转帖】DTF的C代码分析

DFT 公式

傅立叶分析在很多领域都得到了很广泛的应用，使大家能够很方便的在信号的频域进行分析、处理，而离散傅立叶变换（DFT）自然就成为数字信号处理最为基本的知识。DFT的原理有很多资料可以查找的到，在此只给出基本的推导过程。

DFT的C代码之函数声明

上一节中给出了DFT的计算公式，很容易转换成C的代码，下面给出学习用的DFT函数声明、实现以及测试的C代码：

 

文件1：“01_MyDFT.h” DFT函数声明头文件

/***********************************************************   File name:    01_MyDFT.h   Author:       sszhou   Version:      v 0.1.1   Date:         2007-08-07   Description:  DFT study function header file   History:     1. Date:    2007-08-07        Author:  sszhou        Modification: create, add my_dft function ***********************************************************/ #ifndef __MyDFT__H__ #define __MyDFT__H__   typedef long double data_t;   void my_dft_01(     data_t *xr,          // input real part    a(n)     data_t *xi,          // input image part: b(n)     data_t *yr,          // output real part:  A(n)     data_t *yi,          // output image part: B(n)     int N              // N );   #endif /* #ifndef __MyDFT__H__ */

DFT的C代码之函数实现

 

文件2：“01_MyDFT.c” DFT函数实现文件

/***********************************************************   File name:    01_MyDFT.c   Author:       sszhou   Version:      v 0.1.1   Date:         2007-08-07   Description:  DFT analyse study   History:     1. Date:    2007-08-07        Author:  sszhou        Modification: create, add my_dft function ***********************************************************/ #include <math.h> #include "01_MyDFT.h"   #define PI      3.1415926535897932384626433832795   /**   * Function:   my_dft_01   * @brief:     DFT basic function, study used No.01   * @param:   *             xr  ->  [in] input data real part   *             xi  ->  [in] input data image part   *             yr  ->  [out]output data real part   *             yi  ->  [out]output data image part   *             N   ->  window samples N   * @return:   *             null   * Others:   *             Calculate expressions:   *             A(k) = sum(n=0, N-1) {[a(n)*cos(Qnk) + b(n)*sin(Qnk)]}   *             B(k) = sum(n=0, N-1) {[b(n)*cos(Qnk) - a(n)*sin(Qnk)]}   *             Q = 2*pi/N   */ void my_dft_01(     data_t *xr,          // input real part    a(n)     data_t *xi,          // input image part:  b(n)     data_t *yr,          // output real part:  A(k)     data_t *yi,          // output image part: B(k)     int N              // N ) {     int n, k;     data_t Q;       Q = 2*PI/N;       for (k=0; k<N; k++)     {         yr[k] = 0.0;         yi[k] = 0.0;           for (n=0; n<N; n++)         {             yr[k] += xr[n]*cos(Q*n*k) + xi[n]*sin(Q*n*k);             yi[k] += xi[n]*cos(Q*n*k) - xr[n]*sin(Q*n*k);         }     } }

DFT的C代码之测试程序

文件3：“01_MyDFT_Demo.c” DFT函数实现文件

// 01_MyDFT_Demo.c #include <stdio.h> #include <math.h> #include "01_MyDFT.h"   #define DFT_N       32   #define PI2         2*3.1415926535897932384626433832795   #define Fs          16000       // sample rate #define Ft          2000        // signal frequency     void main(void) {     int i;     data_t phase, inc, freq, gain;     data_t xr[DFT_N], xi[DFT_N];     data_t yr[DFT_N], yi[DFT_N];       // 1. Test signal generate     phase = 0;     inc = PI2*Ft/Fs;            // sample rate is Fs, signal frequency is Ft     for (i=0; i<DFT_N; i++)     {         xr[i] = sin(phase)*0.5; // real part is sin wave, gain is -6.02(dB)         xi[i] = 0;              // image part is zero           phase += inc;           // update sin wave phase         if (phase > PI2)         {             phase -= PI2;         }     }         // 2. Do DFT     my_dft_01(xr, xi, yr, yi, DFT_N);         // 3. Print DFT result     printf("DFT Result list:/n");     for (i=0; i<DFT_N; i++)    {         printf("A(%d) = %f, /tB(%d) = %f;/n", i, yr[i], i, yi[i]);     }         // 4. Print frequency analyse result (base on DFT result)     printf("/n/nFrequency and gain analyse:/n");     for (i=0; i<DFT_N/2; i++)     {         freq = i*Fs/DFT_N;      // Calc Frequency         gain = 20*log10(2*pow(yr[i]*yr[i] + yi[i]*yi[i], 0.5)/DFT_N);   // Calc Gain (dB)         printf("Freq: %.2f(Hz), /tGain: %.6f(dB);/n", freq, gain);     }         printf("/n__END__/n"); }

 

DFT的C代码之执行结果

程序执行结果：

DFT Result list: A(0) = -0.000000,       B(0) = 0.000000; A(1) = 0.000000,        B(1) = -0.000000; A(2) = 0.000000,        B(2) = -0.000000; A(3) = 0.000000,        B(3) = -0.000000; A(4) = 0.000000,        B(4) = -8.000000; A(5) = 0.000000,        B(5) = 0.000000; A(6) = 0.000000,        B(6) = 0.000000; A(7) = 0.000000,        B(7) = 0.000000; A(8) = 0.000000,        B(8) = 0.000000; A(9) = -0.000000,       B(9) = 0.000000; A(10) = 0.000000,       B(10) = -0.000000; A(11) = 0.000000,       B(11) = -0.000000; A(12) = -0.000000,      B(12) = 0.000000; A(13) = 0.000000,       B(13) = -0.000000; A(14) = -0.000000,      B(14) = 0.000000; A(15) = -0.000000,      B(15) = 0.000000; A(16) = 0.000000,       B(16) = -0.000000; A(17) = 0.000000,       B(17) = -0.000000; A(18) = 0.000000,       B(18) = 0.000000; A(19) = 0.000000,       B(19) = -0.000000; A(20) = 0.000000,       B(20) = 0.000000; A(21) = -0.000000,      B(21) = -0.000000; A(22) = 0.000000,       B(22) = -0.000000; A(23) = 0.000000,       B(23) = 0.000000; A(24) = 0.000000,       B(24) = 0.000000; A(25) = -0.000000,      B(25) = 0.000000; A(26) = -0.000000,      B(26) = 0.000000; A(27) = -0.000000,      B(27) = 0.000000; A(28) = -0.000000,      B(28) = 8.000000; A(29) = 0.000000,       B(29) = 0.000000; A(30) = 0.000000,       B(30) = 0.000000; A(31) = -0.000000,      B(31) = -0.000000;     Frequency and gain analyse: Freq: 0.00(Hz),         Gain: -326.272234(dB); Freq: 500.00(Hz),       Gain: -324.153486(dB); Freq: 1000.00(Hz),      Gain: -331.771472(dB); Freq: 1500.00(Hz),      Gain: -319.116251(dB); Freq: 2000.00(Hz),      Gain: -6.020600(dB); Freq: 2500.00(Hz),      Gain: -338.683270(dB); Freq: 3000.00(Hz),      Gain: -314.219624(dB); Freq: 3500.00(Hz),      Gain: -306.271348(dB); Freq: 4000.00(Hz),      Gain: -309.032553(dB); Freq: 4500.00(Hz),      Gain: -309.435166(dB); Freq: 5000.00(Hz),      Gain: -312.910650(dB); Freq: 5500.00(Hz),      Gain: -310.372396(dB); Freq: 6000.00(Hz),      Gain: -307.054308(dB); Freq: 6500.00(Hz),      Gain: -311.365771(dB); Freq: 7000.00(Hz),      Gain: -305.817757(dB); Freq: 7500.00(Hz),      Gain: -307.467434(dB);   __END__

 

DFT的C代码之结果分析

DFT的测试信号是一个采样率为16KHz的2K正弦波信号，增益为-6.021dB（0.5倍）。从DFT执行结果可以看出，只有在k=4，28的时候A或B的值较大，其他点的值都太小可以忽略不计。我们只调了一个2K的正弦波，为什么会有2个不同的点有值，而且值是相同的呢？这是因为N/2以上的点是信号的镜像能量。

DFT分析的是频域的分析，可是X(k)对应的是什么频率的值呢？这和信号的采样率以及作DFT采用的点数N有关，具体的计算公式为：Freq(k)=k*采样率/N（k<N/2）；以测试程序为例，第4个点（k=4）代表的频率值是：4*16000/32 = 2000Hz，这和我们测试信号频率是一至的。当k>=N/2时计算出来的频率会大于乃奎斯特频率，这是因为k>=N/2以上的点都是镜像的点，因此k>=N/2点的频率的计算方法为：Freq(k) = (N-k)*采样率/N (k>=N/2)（从计算频率的公式可以看出，N值越大所描述的频率就越精细，当“N=采样率”的时候就能够精确到1Hz）。

知道了k点对应的频率值，又如何计算k点对应频率的增益呢？从公式：X(k)=A(k)+jB(k)可以看出|X(k)| = (A(k)^2 + B(k)^2)^0.5；由于有一半的镜像能量的分布，因此在计算增益的时候需要增加一倍。同时是对N个点做的DFT，因此该值中包含了N个点总的增益。得到k点频率对应的增益计算公式为：20×log10（2×|X(k)|/N）；

根据k点的频率以及增益的计算公式得到了DFT分析结果中各个点的频率和相应的增益信息。

由测试结果可以清晰的看到测试信号的分析结果为2KHz处的频率增益是-6.02dB和我们产生的信号完全一致～～