一、DFT计算公式
这里就不对DFT概念进行叙述,直接上计算公式
其中N为DFT点数。
公式如此,但是在程序中并非如此运算,而是利用欧拉公式对DFT的计算公式进行了转化
转换后公式变为
利用转化后的公式进行程序编写就变得思路清晰起来。
在进行频谱分析时有一个概念叫做频谱分辨率,其计算公式如下:
∆f为频谱分辨率,N为采样点数,Ts为采样时间间隔,fs为采样频率
比如频谱分辨率为50Hz,那么计算结果中第一个频点(k=1)时对应为50Hz分量,第二个频点(k=2)时对应为100Hz分量,依此类推。k=0时为直流分量。
二、DFT程序实现
按照惯例,直接上程序
#define Pi 3.1416 // DFT计算时的Π值
/*
*==============================================================================
*函数名称:Dft_Calculate
*函数功能:计算DFT
*输入参数:Input:数据数组 N:DFT点数,Num:获取前Num个频点
*输出:前Num频点的幅值
*返回值:无
*============================================================================*/
void Dft_Calculate(int Input[],int N,int Num)
{
int n = 0; // 计算单个点DFT的循环变量
int k = 0; // 计算全部点DFT的循环变量
double Real[N],Imag[N];
double Real_Sum = 0; // 存储实部的和
double Imag_Sum = 0; // 存储虚部的和
double Dft_Result[N]; // 存储DFT结果
// 测试一下int和double的误差
for (k = 0;k < Num;k ++)
{
for (n = 0;n < N;n ++)
{
Real[n] = cos(((2*Pi)/N)*n*k)*Input[n];
Imag[n] = -sin(((2*Pi)/N)*n*k)*Input[n];
Real_Sum = Real_Sum + Real[n]; // 实部求和
Imag_Sum = Imag_Sum + Imag[n]; // 虚部求和
}
Dft_Result[k] = sqrt(Real_Sum * Real_Sum + Imag_Sum * Imag_Sum);
Real_Sum = 0;
Imag_Sum = 0; // 初始化计算值
}
}
该函数逻辑简单,注释明确,这里就不再赘述了。主要功能是输入需要求DFT的数组Input后,输入DFT点数以及想获取的结果长度,即可计算DFT结果值。
值得注意的是该函数计算结果是一个局部变量,在计算完成后会被释放。如果在计算完成后需要对计算结果进行进一步操作,可以将存储计算结果的数组定义放在函数外面。