随机信号处理-P2-FFT

我通过函数递归实现DIT-FFT并完成对延拓的避免（只取用被递归到的数据，余下数据直接舍弃，或者调整FFT位数并补零），并在处理第一个蝶形运算时同时完成倒位二进制的重新排列。

闲话少叙，放代码：

#include<cstdio>
#include<iostream>
#include<string>
#include<cstdlib>
#include<ctime>
#include<cmath>

using namespace std;
#define N1 512//x序列长度
#define TrueN2 519//真正的抽样点数
#define MaxPeak 100
#define Pai 3.141592653589793238462643383279
#define sold 500
//#define Pai 5//设想：如果pai不是π，会怎样？ 

//定义区 
double x[1025], x_a[2], x_w[2], x_e[2], x_ghals[1025];
double X[3][1025], W[2][1025];
int Eustia,DFThighest_k[MaxPeak],N2;
double sigema_ghals,DFThighest_X[MaxPeak],findf[MaxPeak],mse1[2][21];

//声明区 
int ran(int ran_m, int ran_n, int ran_seek);
double makeghals(int& j);
double rand_back(double i,double j);
void creatx();
void extrax();
void prepareW(int n_W);
void DIT_FFTEustia(int yuan, int wing);//Eustia算法，包括秽羽和展翼两个部分
void DIT_FFT();
void findf_easy(int num_f);
double ran_double(double ran_m, double ran_n, int ran_seek);

int main(void)
{
	printf("Hello,Fred Zhou!\n");
	printf("x采样点数N1=%d\nX抽样点数N1=%d\n",N1,TrueN2);
	double rubbish;
	N2=int(pow(2,ceil(log(TrueN2)/log(2))));//FFT预处理 
	prepareW(N2);
	sigema_ghals=0;
	x_a[0]=2;
	x_a[1]=3;
	x_w[0]=0.4*2*Pai;
	x_w[1]=0.3*2*Pai;
	x_e[0]=0*Pai;
	x_e[1]=0.5*Pai;
	creatx();
	printf("Creat x:\n\n\tx(t)=%.2lfsin(%.2lft+%.2lf)+",x_a[0],x_w[0],x_e[0]);
	printf("%.2lfsin(%.2lft+%.2lf)+N(0,%.2lf^2)\n",x_a[1],x_w[1],x_e[1],sigema_ghals);	
	DIT_FFT();
	findf_easy(2);
	printf("\n解算得到角速度：\n第一条：%.2lf*2Pai  第二条：%.2lf*2Pai\n",findf[0],findf[1]);	
	//展示x与X
	freopen("FFT_X.out","w",stdout);
	for (int i=0;i<TrueN2;i++)
	{
		cout<<X[2][i]<<endl;
	}
	freopen("信号x.out","w",stdout);
	for (int i=0;i<N1;i++)
	{
		cout<<x[i]<<endl;
	}
	freopen("CON","w",stdout);
	printf("请在FFT_X.out提取数据\n");
	printf("//致敬经典:去防疫局门口点一曲《想要拥有羽毛和翅膀》\n");
	return 0;
}

void creatx()
{
	int i=0,j=0;
	for(;i<N1;i++)
	{
		x_ghals[i]=rand_back(0,sigema_ghals*sigema_ghals);
		x[i]=x_a[0]*sin(x_w[0]*i+x_e[0])+x_a[1]*sin(x_w[1]*i+x_e[1])+x_ghals[i];
	}
	return ;
}
//为了补0凑fft，可以N2+1再log2再+1

int ran(int ran_m, int ran_n, int ran_seek) 
{
	srand((unsigned)time(NULL) * ran_seek);
	int ran_int = rand() % (ran_n - ran_m) + ran_m;
	return ran_int;
}

double ran_double(double ran_m, double ran_n, int ran_seek)
{ 
	srand((unsigned)time(NULL) * ran_seek);
	int ran_int = rand()%(10000);
	double ran_dou=ran_int/10000.0;
	ran_dou*=(ran_n-ran_m);
	ran_dou+=ran_m;
	return ran_dou;
}


double rand_back(double i,double j)//高斯创造
{
	double u1=double((rand()%29999)+1.0)/30000,u2=double(rand()%30000)/30000,r;
	static unsigned int seed=0;
	r=i+sqrt(j)*sqrt(-2.0*(log(u1)/log(exp(1.0))))*cos(2*Pai*u2);
//MENTION!这里u1有可能=0导致报错，我找了整整一晚上！
//而且 log(u1)=1也会莫名报错 
	//printf("t=%lf\n",(log(u1)));//调试：显示创造的随机值 
	return r;
}

void findf_easy(int num_f)
{
	int i,j;
	double find_save=X[2][0];
	X[2][0]=0;
	for(i=0;i<MaxPeak;i++)
	{
		DFThighest_X[i]=0;
		DFThighest_k[i]=-1;
	}
	double last;
	for(i=1;i<TrueN2/2;i++)
	{
		if(X[2][i]<=X[2][i-1])
			continue;
		if(X[2][i]<=X[2][i+1])
			continue;
		if(X[2][i]<=DFThighest_X[num_f-1])
			continue;
		for(j=num_f-2;j>-1;j--)
		{
			if(X[2][i]>DFThighest_X[j])
			{
				DFThighest_X[j+1]=DFThighest_X[j];
				DFThighest_k[j+1]=DFThighest_k[j];
			}
			else
			{
				DFThighest_X[j+1]=X[2][i];
				DFThighest_k[j+1]=i;
				break;
			}
		}
		if(j<0)
		{
			DFThighest_X[0]=X[2][i];
			DFThighest_k[0]=i;
		}
	}
	for(j=num_f-1;j>-1;j--)
	{
		findf[j]=(DFThighest_k[j]);
		findf[j]/=TrueN2;
	}
	X[2][0]=find_save;
	return ;
}

void prepareW(int n_W)
{
	for (int i = 0; i < n_W/2; i++)
	{
		W[0][i] = cos((2 * Pai * i) / n_W);
		W[1][i] = -sin((2 * Pai * i) / n_W);
	}
	return;
}

//自主设计的Eustia式递归组
void DIT_FFTEustia(int yuan, int wing)
//翅膀的根和大小
{
	wing /= 2;
	int Tia;
	if (wing != 1)
	{
		DIT_FFTEustia(yuan, wing);
		DIT_FFTEustia(yuan + wing, wing);
	}
	else
	{
//秽羽
		int tag = N2/2, up = 0;
		Tia = Eustia;//Tia是Eustia的小名
		for (; Tia > 0;)
		{
			tag /= 2;
			up += tag * (Tia % 2);
			Tia /= 2;
		}
		tag = up + N2 / 2;
		X[0][yuan] = x[up] + x[tag];
		X[1][yuan] = 0;
		yuan++;
		X[0][yuan] = x[up] - x[tag];
		X[1][yuan] = 0;
		Eustia++;//Eustia吸取伊莲的天使能
		return;
	}
//展翼
	double XX[3][2];//备用
	Tia = N2 / (wing*2);
	for (int i = 0; i < wing; i++)
//C=W(i*(N2/2wing),N2)=W[N2i/2wing]=W[iTia]
	{
		XX[0][0] = X[0][i + yuan + wing] * W[0][i * Tia] - X[1][i + yuan + wing] * W[1][i * Tia];
		XX[0][1] = X[0][i + yuan + wing] * W[1][i * Tia] + X[1][i + yuan + wing] * W[0][i * Tia];
		XX[1][0] = X[0][i + yuan] + XX[0][0];
		XX[1][1] = X[1][i + yuan] + XX[0][1];
		XX[2][0] = X[0][i + yuan] - XX[0][0];
		XX[2][1] = X[1][i + yuan] - XX[0][1];
		X[0][i + yuan] = XX[1][0];
		X[1][i + yuan] = XX[1][1];
		X[0][i + yuan + wing] = XX[2][0];
		X[1][i + yuan + wing] = XX[2][1];
	}
	Tia=sold;//经典...别删
	return;
}

void DIT_FFT()
{
	Eustia = 0;
	DIT_FFTEustia(0, N2);
	for (int i = 0; i < N2; i++)
	{//求模 
		X[2][i] = sqrt(X[0][i] * X[0][i] + X[1][i] * X[1][i]);
	}
	return;
}

运行结果：

同时，展现出了相当强大的抗噪声性能：

哪怕噪声远超过任何一个信号，依旧可以稳定滤波。