DFT C语言实现

谁借莪１个温暖的怀抱￠ 2022-01-27 10:13 240阅读 0赞

**DFT的描述**

DFT即离散傅里叶变换，是傅里叶变换在时域和频域上都呈离散的形式，将信号的时域采样变换为其DTFT的频域采样。在形式上，变换两端（时域和频域上）的序列是有限长的，而实际上这两组序列都应当被认为是离散周期信号的主值序列。即使对有限长的离散信号作DFT，也应当将其看作其周期延拓的变换。

**DFT的公式**

X(k)= ∑ n = 0 N − 1 x ( n ) e − j 2 π N n k \\sum\_\{n=0\}^\{N-1\}x(n)\{\\rm e\}^\{-j\\frac\{2π\}\{N\}nk\} ∑n=0N−1x(n)e−jN2πnk

其中，X(k) 是 DFT后的结果

x(n) 是采样信号（一般x(n)是实信号，虚部为0）

**1.根据欧拉公式：** e j θ = c o s θ + j s i n θ \{\\rm e\}^\{jθ\}=cosθ+jsinθ ejθ=cosθ\+jsinθ ，可代换得到：

X(k)= ∑ n = 0 N − 1 x ( n ) c o s 2 π N n k − j x ( n ) s i n 2 π N n k \\sum\_\{n=0\}^\{N-1\}x(n)\{cos\\frac\{2π\}\{N\}nk\}-jx(n)\{sin\\frac\{2π\}\{N\}nk\} ∑n=0N−1x(n)cosN2πnk−jx(n)sinN2πnk ， （其中前半部分为实部，后半部分为虚部）

**2.正余弦合并：** asin(x)+bcos(x)=  a 2 + b 2 \\sqrt\{a^2+b^2\} a2\+b2sin(x+arctanb/a)

**DFT计算程序**

#include <stdio.h>
    #include <stdlib.h>
    #include<math.h>
    #define pi 3.14
    #define N 8000//8000个点
    
    typedef struct
    {
        double real,imag;
    } complex;//复数
    complex dft_out[8020];//单个点计算k
    complex dft_one[8020];//单个点计算n
    double amp[N];//DFT的结果
    int main()
    {
        int n,k;
        double xn;
        for(k=0; k<N; k++)//k循环
        {
            for(n=0; n<N; n++)//n循环
            {
                xn=0.6*sin(n*pi*100)+0.6*sin(n*pi*1000);//原信号
                dft_one[n].real=xn*cos(2*pi/N*n*k);//实部
                dft_one[n].imag=xn*sin(2*pi/N*n*k);//虚部
                dft_out[k].real+=dft_one[n].real;
                dft_out[k].imag+=dft_one[n].imag;//DFT后的实部'虚部相加
            }
            amp[k]=sqrt(dft_out[k].real*dft_out[k].real+dft_out[k].imag*dft_out[k].imag);//正余弦合并（根号下平方之和）
            printf("%d  %f\n",k,amp[k]);
        }
    }

**gnuplot作图：**

![6][]

局部细节放大：

![7][]

[6]: /images/20220127/1238efcbfb1d4ac0b2d11d4c86066bc4.png
[7]: /images/20220127/7a0f35ff14cc4c5580b7044a27a8c16f.png