cpp实现fft-dit算法，源码。

因为博主懒不想给自己的pc装大头软件matlab，虽然我的hp二代顶配光影精灵性能还ok，但毕竟是带着双系统，存储空间分给Linux一半。所以大软件我一般只装在实验室电脑上，恰巧随机信号的老师给我们讲课让我们练手写dft。我想了一下一个二层循环也太low了也没啥可写的，就开启了不做死就不会死的递归之路。之前也写过一些递归函数，简洁那是相当的简洁。反正我自信dsp学的还说的过去，那就写写fft呗，说不定还能得到老师的加分😀。

我确实得承认我太作了，这个函数调试花了我一个星期时间，一度都进入了呕吐式一句一断执行。好在vs作为宇宙第一ide能让我实时看到每个变量的值，最终我终于是调好了。fft的算法原理很多地方都有介绍，这里我就不废话。就实现起来来说我觉得最难的就是跳跃式数据检索（如果写dif算法那就是跳跃式数据
存储），调试期间我一度搞不清存储位置offset。写完了以后才恍然顿悟，理解了书上所写的层级的概念。不得不说像这种二叉树递归，深度是一个极其重要的概念。

下面上代码：

首先是对<complex.h>中的一些运算符重载，这里建议使用管理员身份运行ide，否则没有更改库函数的权限。

//手动修改<complex.h>中对_C_double_complex结构体的定义。
typedef struct _C_double_complex
{
	double _Val[2];
	

	void multiply(_C_double_complex num)
	{
		double tem = this->_Val[0];
		this->_Val[0] = this->_Val[0] * num._Val[0] - this->_Val[1] * num._Val[1];

		this->_Val[1] = this->_Val[1] * num._Val[0] + tem * num._Val[1];

	}
	void init(double re, double im)
	{
		this->_Val[0] = re;
		this->_Val[1] = im;
	}

	_C_double_complex operator + (_C_double_complex num)
	{
		num._Val[0] = num._Val[0] + this->_Val[0];
		num._Val[1] = num._Val[1] + this->_Val[1];
		return num;
	}
	_C_double_complex operator - (_C_double_complex num)
	{
		num._Val[0] = this->_Val[0] - num._Val[0];
		num._Val[1] = this->_Val[1] - num._Val[1];
		return num;
	}

} _C_double_complex;

接下来是fft.h

#pragma once
#ifndef __FFT
#define __FFT 
#endif // !__FFT
#include <math.h>
#include <complex.h>
#define M_PI 3.14159265358979323846
enum c_cal_state
{
	error, success
};
c_cal_state FFT_DIT(int N, int SIZE_k, _C_double_complex* xn, _C_double_complex* xk,int recall );

下面是fft.c

#include "pch.h"
#include "fft.h"

//#define __RangeXn(name,n,N) n>N?0:*(name+n)


/*
@Function Name:FFT_DIT
@Aurther:QyShen
@Email：shen99855@outlook.com
@Date of version:3/04/2020
@param: N     :xn 长度
		SIZE_k:fft 点数
		xn：  xn基址，只读访问
		xk：  xk基址，读写访问
		recall：递归调用标志，首次调用时，取0，递归调用时，取1；用户应当按recall=0来调用。
@return：error：错误，应考虑停止计算。
		success：成功
@bug说明：目前来讲，应选取N是2的幂，且不要让N！=SIZE_k;否则可能引起xn数组的越界访问，在今后的版本中，会加以修复。
*/
c_cal_state FFT_DIT(int N, int SIZE_k, _C_double_complex* xn, _C_double_complex* xk, int recall)
{

	static int depth = 0;
	static int step = 1;
	static unsigned int bitt;
	static int Commen_N;
	static int base;
	if (SIZE_k < N)
		return error;
	else
	{
		if (recall == 0)
		{

			depth = 0;
			step = 1;
			bitt = 0x80000000;
			Commen_N = N;
			while (!(bitt&N))
				bitt >>= 1;
			if (bitt - N)
				bitt <<= 1;//add 0 to end terminal ,ensure that N=power of 2;
			else
				;
			if (SIZE_k < bitt)
			{
				return error;
			}
			else
				;

			FFT_DIT(bitt, bitt, xn, xk, 1);
			return success;

		}
		else
		{
			depth++;
			step <<= 1;
			if (N > 2)
			{

				_C_double_complex* xe = xn;
				_C_double_complex* xo = xn + step / 2;
				FFT_DIT(N / 2, N / 2, xe, xk, 1);

				depth--; step >>= 1;
				FFT_DIT(N / 2, N / 2, xo, xk + N / 2, 1);
				depth--; step >>= 1;

				for (int i = 0; i < N / 2; i++)
				{
					_C_double_complex teml = xk[i];
					_C_double_complex temr = xk[i + N / 2];
					_C_double_complex temwei;
					temwei.init(0, -2 * M_PI / N * i);
					temr.multiply(cexp(temwei));
					xk[i] = teml + temr;
					xk[i + N / 2] = teml - temr;
				}


			}
			else
			{
				_C_double_complex teml = xn[0];
				//std::cout <<   xn <<".";
				_C_double_complex temr = xn[step / 2];
				//if (step / 2 > Commen_N)
				//	temr.init(0, 0);//判断越界。
				//else
				//	;
				_C_double_complex temwei;
				temwei.init(0, -2 * M_PI / Commen_N * 0);
				temr.multiply(cexp(temwei));
				xk[0] = teml + temr;
				xk[1] = teml - temr;

				return success;

			}
		}

	}

}

简单的例子main.cpp

#include <iostream>
#include "fft.h"

#define __N 256
#define __K 256
using namespace std;
int main()
{
	_C_double_complex xn[__N];
	_C_double_complex xk[__K];
	for (int i = 0; i < __N; i++)
	{
		xn[i].init(sin(2 * M_PI*0.1*i + M_PI / 3) + 10 * sin(2 * M_PI*0.3*i + M_PI / 4), 0);

	}
	FFT_DIT(__N, __K, xn, xk, 0);
	for (int i = 0; i < __K; i++)
	{
		cout << "n=" << i << " xn=" << (xn + i)->_Val[0] << endl;
		cout << "k=" << i << "   re(xk)=" << real(xk[i]) << "  im(xk)=" << imag(xk[i]) << "   amp(xk)=" << abs(xk[i]) << endl;

	}
	system("pause");
	return 0;
}