fft.cpp

FFT 算法的动态链接库 自己做的 算法复杂度一般

#ifndef FFT_CPP
#define FFT_CPP

#include <math.h>
#define PI 3.1415926
#include "Complex.h"
#include "fft.h"
int Pow(int q ,int m)//int 型别数据求幂
{
	int res;
	res = 1;
	while (m != 0)
	{
		res *= q;
		m--;
	}
	return res;
}
//倒序算法
void InverseSort(int *sortin,int M)
{  
	if (M == 1)
	{
		return;
	}
	int n = (int) Pow(2,M);//求出需要排列的个数
	int sum = sortin[0] + sortin[n - 1];//这个数是对称数字相加的和
	int i;
	for (i = 0 ; i < (n / 2);i++)
	{
		sortin[i] = sortin[ 2 * i];
	}
	InverseSort(sortin,M - 1);
	for (i = 0 ; i < (n / 2);i++)
	{
		sortin[i + n / 2] = sum - sortin[n / 2 - 1 -i];
	}
	return ;
}

void Base2TimeFFT(Complex * orign_result,Complex * l_result, int M,int Isign)
{
	//1 mean fft  -1 mean ifft
	if (!(Isign == 1 || Isign == -1 ))//判断输入是否规范
	{
		return ;
	}	
	Isign = -1 * Isign;//为了后来运算思路清晰 	

	int Nsample = (int) Pow(2,M);
	//////////////////////////////////////////////////////////////////////////
	//倒序	
	int i, j;
	int *inverse_sort;
	inverse_sort = new int[Nsample];
	for (i = 0 ; i < Nsample ; i++)
	{
		inverse_sort[i] = i;
	}
	InverseSort(inverse_sort,M);
	for (i = 0; i < Nsample ; i++)
	{
		l_result[i] = orign_result[(inverse_sort[i])];
	}
	delete[] inverse_sort;



	//计算旋转因子
	Complex *spinelement;
	spinelement = new Complex[Nsample / 2];//用到的旋转因子只有一半
	double N2pi = PI * 2 / Nsample;
	for (i = 0; i < Nsample / 2;i++)
	{
		spinelement[i].m_real = cos(N2pi * i);
		spinelement[i].m_img = Isign *  sin(N2pi * i);
	}
	
	int spinstep;//这个参数代表了该级运算的旋转因子的个数
	//同时也代表了 碟性的跨度
	//运算循环开始
	int L ;
	int spin;
	int temp;
	Complex temp_save1,temp_save2;
	for (L = 1; L <= M;L++)//第L级运算
	{
		temp = (int) Pow(2 , L);
		spinstep = (int) Pow(2,L-1);
		for (j = 0; j < spinstep;j++)
		{
			spin = (int) (j * Pow(2,(M - L)));
			for (i = j ; i <= (Nsample - 1); i+= temp)
			{
				temp_save1 = (l_result[i ] +  
					(l_result[i + spinstep ] * spinelement[spin]));
				temp_save2 = ( l_result[i ] - 
					(l_result[i + spinstep ] * spinelement[spin]));
				l_result[i ] = temp_save1;
				l_result[i + spinstep ] = temp_save2;
			}
		}
	}
	delete[]  spinelement;	
	if (Isign == 1)//如果是IFFT还要归一化 isign 已经在前面反号了的
	{
		for (i = 0 ; i < Nsample;i++)
		{
			l_result[i] /= Nsample;
		}
	}
}

























#include <math.h>
#define PI 3.1415926
/*
按时间抽选的基2-FFT算法
void FFT(double data[] , int Nsample,int Isign)
Nsample  整型变量 【in】 采样点数 要求Nsample为2的幂次
data【】 【in / out】 2 * Nsample的数组，

  倒序算法 用递归实现

  旋转因子 说明：
  一个N级的序列 有1-N的级别运算 L级别有2的L-1次方-1个旋转因子 
  W(N,P),P = 2的log（2，N) - L次方 * J    J 从0到2的L-1次方-1 

 
*/
//typedef double DATATYPE;
//class Complex//定义一个复数类
//{
//public:
//	DATATYPE m_real;
//	DATATYPE m_img;
//public:
//	Complex()
//	{
//		m_real = (DATATYPE)0 ;
//		m_img = (DATATYPE) 0;
//	}
//	Complex(DATATYPE real,DATATYPE img)
//	{
//		m_img = img;
//		m_real = real;
//	}
//	Complex& operator=(const Complex &plex)
//	{
//		m_real = plex.m_real;
//		m_img = plex.m_img;
//		return *this;
//	}
//};
////复数类的几个操作符重载函数
//const Complex operator*(const Complex plex1,const Complex plex2)
//{
//	Complex l_result;
//	l_result.m_real = plex1.m_real * plex2.m_real - plex1.m_img * plex2.m_img;
//	l_result.m_img = plex1.m_real * plex2.m_img + plex1.m_img * plex2.m_real;
//	return l_result;
//}
//
//const Complex operator+(const Complex plex1,const Complex plex2)
//{
//	Complex l_result;
//	l_result.m_real = plex1.m_real + plex2.m_real;
//	l_result.m_img = plex1.m_img + plex2.m_img;
//	return l_result;
//}
//
//const Complex operator-(const Complex plex1,const Complex plex2)
//{
//	Complex l_result;
//	l_result.m_real = plex1.m_real - plex2.m_real;
//	l_result.m_img = plex1.m_img - plex2.m_img;
//	return l_result;
//}
////倒序算法
//void InverseSort(int *sortin,int M)
//{  
//	if (M == 1)
//	{
//		return;
//	}
//	int n = (int) pow(2,M);//求出需要排列的个数
//	int sum = sortin[0] + sortin[n - 1];//这个数是对称数字相加的和
//	int i;
//	for (i = 0 ; i < (n / 2);i++)
//	{
//		sortin[i] = sortin[ 2 * i];
//	}
//	InverseSort(sortin,M - 1);
//	for (i = 0 ; i < (n / 2);i++)
//	{
//		sortin[i + n / 2] = sum - sortin[n / 2 - 1 -i];
//	}
//	return ;
//}
//
//void Base2TimeFFT(Complex * orign_result,Complex * l_result, int M,int Isign = 1)
//{
//	//1 mean fft  -1 mean ifft
//	if (!(Isign == 1 || Isign == -1 ))//判断输入是否规范
//	{
//		return ;
//	}
//	Isign = -1 * Isign;//为了后来运算思路清晰 
//
//
//	int Nsample = (int) pow(2,M);
//	//////////////////////////////////////////////////////////////////////////
//	//倒序	
//	int i, j;
//	int *inverse_sort;
//	inverse_sort = new int[Nsample];
//	for (i = 0 ; i < Nsample ; i++)
//	{
//		inverse_sort[i] = i;
//	}
//	InverseSort(inverse_sort,M);
//	for (i = 0; i < Nsample ; i++)
//	{
//		l_result[i] = orign_result[(inverse_sort[i])];
//	}
//	delete[] inverse_sort;
//
//
//
//	//计算旋转因子
//	Complex *spinelement;
//	spinelement = new Complex[Nsample / 2];//用到的旋转因子只有一半
//	double N2pi = PI * 2 / Nsample;
//	for (i = 0; i < Nsample / 2;i++)
//	{
//		spinelement[i].m_real = cos(N2pi * i);
//		spinelement[i].m_img = Isign *  sin(N2pi * i);
//	}
//	
//	int spinstep;//这个参数代表了该级运算的旋转因子的个数
//	//同时也代表了 碟性的跨度
//	//运算循环开始
//	int L ;
//	int spin;
//	int temp;
//	Complex temp_save1,temp_save2;
//	for (L = 1; L <= M;L++)//第L级运算
//	{
//		temp = (int) pow(2 , L);
//		spinstep = (int) pow(2,L-1);
//		for (j = 0; j < spinstep;j++)
//		{
//			spin = (int) (j * pow(2,(M - L)));
//			for (i = j ; i < Nsample; i+= temp)
//			{
//				temp_save1 = (l_result[i] + 
//					(l_result[i + spinstep] * spinelement[spin]));
//				temp_save2 = ( l_result[i] - 
//					(l_result[i + spinstep] * spinelement[spin]));
//				l_result[i] = temp_save1;
//				l_result[i + spinstep] = temp_save2;
//			}
//		}
//	}
//	delete[]  spinelement;	
// }
#endif