C#入门 公式法和FFT法求解DTMF信号的DFT

一、涉及数字信号处理公式

离散傅里叶变换定义式
$$X[k]=\sum _{n=0}^{N-1}x[n]e^{-j\frac{2\pi kn}{N}}$$
令
$$W_N=e^{-j\frac{2\pi }{N}}$$
则
$$X[k]=\sum _{n=0}^{N-1}x[n]W_N^{kn}$$
x1[n]由x[n]的奇数项组成
x2[n]由x[n]的偶数项组成
有
$$X[k]=\sum _{n=0}^{\frac{N}{2}-1}{x}_1[n]W_{N/2}^{kn}+W_N^k\sum _{n=0}^{\frac{N}{2}-1}{x}_2[n]W_{N/2}^{kn}$$
即
$$X[k]=X_1[k]+W_N^k{X}_2[k]$$

二、 程序源码

using MathNet.Numerics;
using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Diagnostics;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using System.Windows.Forms;

namespace dft
{
    public partial class Form1 : Form
    {
        public Form1()
        {
            InitializeComponent();
        }


        private void button1_Click(object sender, EventArgs e)
        {
            //使用Math.net类库


            //生成测试信号
            //DTMF信号1
            //697Hz+1209Hz
            //采样率使用5K
            List<Complex32> signal1 = new List<Complex32>();
            for (int i = 0; i < 5000; i++)
            {
                double data = Math.Sin(2 * Math.PI * 697 * i / 5000) + Math.Sin(2 * Math.PI * 1209 * i / 5000);
                Complex32 complexData = new Complex32((float)data,0);
                signal1.Add(complexData);
            }

            //绘制信号
            for (int i = 0; i < 100; i++)
            {
                chart1.Series[0].Points.AddXY(((double)i/5000)+"s", Complex32.Abs(signal1[i]));
            }

            //扩展信号长度到2的整数次幂
            List<Complex32> signal1_expand = expand(signal1);


            Stopwatch stopWatch = new Stopwatch();
            stopWatch.Start();

            TimeSpan span = stopWatch.Elapsed;
            double endTime0 = span.TotalMilliseconds;

            //DFT变换(使用公式)未扩展信号
            List<Complex32> dftResult1 = dft(signal1);

            //stopWatch.Stop();
            span = stopWatch.Elapsed;
            //Console.WriteLine(span.TotalMilliseconds);
            double endTime1 = span.TotalMilliseconds;

            //DFT变换(使用公式)扩展长度信号
            List<Complex32> dftResult2 = dft(signal1_expand);
            span = stopWatch.Elapsed;
            //Console.WriteLine(span.TotalMilliseconds);
            double endTime2 = span.TotalMilliseconds;

            //FFT
            List <Complex32> fftResult = fft(signal1_expand);
            span = stopWatch.Elapsed;
            //Console.WriteLine(span.TotalMilliseconds);
            double endTime3 = span.TotalMilliseconds;

            Console.WriteLine(endTime1-endTime0);
            Console.WriteLine(endTime2-endTime1);
            Console.WriteLine(endTime3-endTime2);

            //绘制dft
            for (int i = 0; i < dftResult1.Count/2; i++)
            {
                chart2.Series[0].Points.AddXY(i * 5000 / dftResult1.Count , Complex32.Abs(dftResult1[i]));
            }
            for (int i = 0; i < dftResult2.Count / 2; i++)
            {
                chart2.Series[1].Points.AddXY(i * 5000 / dftResult2.Count, Complex32.Abs(dftResult2[i]));
            }

            //绘制fft
            for (int i = 0; i < fftResult.Count / 2; i++)
            {
                chart3.Series[0].Points.AddXY(i * 5000 / fftResult.Count, Complex32.Abs(fftResult[i]));
            }
        }


        长度扩充到2的整数次幂
        public List<Complex32> expand(List<Complex32> signal)
        {
            int targetLength = 2;
            while (targetLength < signal.Count)
            {
                targetLength = targetLength * 2;
            }

            //补0
            while (targetLength != signal.Count)
            {
                signal.Add(new Complex32(0, 0));
            }

            return signal;
        }

        //使用公式
        public List<Complex32> dft(List<Complex32> signal)
        {
            List<Complex32> dftResult = new List<Complex32>();


            int N = signal.Count;

            for (int k=0; k<signal.Count;k++)
            {
                Complex32 dftItem = new Complex32(0,0);
                for (int n=0; n<signal.Count;n++)
                {
                    Complex32 temp = new Complex32(0,(float)(-2 * Math.PI * k * n / N));
                    temp = Complex32.Multiply(signal[n], Complex32.Exp(temp));

                    dftItem = Complex32.Add(dftItem, temp);
                }
                dftResult.Add(dftItem);
            }

            return dftResult;
        }


        //fft
        //递归方式实现
        //signal的长度为2的整数次幂
        public List<Complex32> fft(List<Complex32> signal)
        {
            List<Complex32> fftResult = new List<Complex32>();
            List<Complex32> tempList = new List<Complex32>();

            List<Complex32> list1 = new List<Complex32>();
            List<Complex32> list2 = new List<Complex32>();


            if (signal.Count == 2)
            {
                return dft(signal);
            }

            //奇偶分解
            for (int i=0; i<signal.Count;i++)
            {
                if (i % 2 == 0)
                {
                    list1.Add(signal[i]);
                }
                else
                {
                    list2.Add(signal[i]);
                }
            }

            List<Complex32> fftResult1 = fft(list1);
            List<Complex32> fftResult2 = fft(list2);

            int N = signal.Count;
            for(int k=0; k<fftResult1.Count; k++)
            {
                Complex32 temp = new Complex32(0, (float)(-2 * Math.PI * k / N));
                Complex32 wnkfftResult2 = Complex32.Multiply(fftResult2[k], Complex32.Exp(temp));

                //Complex32 tempData = fftResult1[k] + wnkfftResult2;

                fftResult.Add(Complex32.Add(fftResult1[k],wnkfftResult2));
                tempList.Add(Complex32.Add(fftResult1[k],Complex32.Negate(wnkfftResult2)));
            }

            fftResult.AddRange(tempList);
            return fftResult;
        }




    }
}

三、 波形截图

四、 耗时分析

6285.2976
DFT 耗时6274.406
FFT 耗时13.1182000000008
FFT效率远高于公式法直接计算