【转帖】dtf的c代码分析...
时间：2010-08-09 来源：Eason_Hsu
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和我们产生的信号完全一致～～
相关阅读 更多 +
