三組傅立葉變換和反變換的程式碼

#include <stdio.h>

#include <math.h>

#include <stdlib.h>

#include <iostream>

#include <fstream>

#include <vector>

#include  <ctime>

#include <iomanip>

using
namespace  std;



#define PI 3.14159265 

#define N  4095       //取樣256次



typedef
struct               //定義結構體

{

	double
	real;/*實部*/

	double
		img;/*虛部*/

}complex;



complex x[N * 2], *W;                       /*輸出序列的值*/
//序列長度 全域性變數



void
add(complex a, complex b, complex *c)      //  複數加運算

{

	c->real = a.real + b.real;

	c->img = a.img + b.img;

}

void
sub(complex a, complex b, complex *c)      //  複數減運算

{

	c->real = a.real - b.real;

	c->img = a.img - b.img;

}

void
mul(complex a, complex b, complex *c)      //複數乘運算

{

	c->real = a.real*b.real - a.img*b.img;

	c->img = a.real*b.img + a.img*b.real;

}

void
divi(complex a, complex b, complex *c)      //  複數除運算

{

	c->real = (a.real*b.real + a.img*b.img) / (b.real*b.real + b.img*b.img);

	c->img = (a.img*b.real - a.real*b.img) / (b.real*b.real + b.img*b.img);

}



void
initW(int
size)                                //尤拉公式運算  

{

	int
		i;

	W = (complex*)malloc(sizeof(complex)* size);    //分配記憶體空間

	for
		(i = 0; i<size; i++)

	{

		W[i].real = cos(2 * PI / size*i);

		W[i].img = -1 * sin(2 * PI / size*i);

	}

	/*for (i = 0; i < size; i++)

	{

	cout << endl;

	cout << endl;

	cout << endl;

	cout << W[i].real << " ";

	cout << W[i].img <<"j" << " ";

	}*/

}



void
changex(int
size)                      //變址運算 

{

	complex temp;

	unsigned int
		i = 0, j = 0, k = 0;

	double
		t;

	for
		(i = 0; i<size; i++)

	{

		k = i; j = 0;

		t = (log(size) / log(2));

		while
			((t--)>0)

		{

			j = j << 1;

			j |= (k & 1);

			k = k >> 1;

		}

		if
			(j>i)

		{

			temp = x[i];

			x[i] = x[j];

			x[j] = temp;

		}

	}

}





void
fftx()                             //快速傅立葉函式

{

	long
		int  i = 0, j = 0, k = 0, l = 0;

	complex up, down, product;

	changex(N);

	for
		(i = 0; i<log(N) / log(2); i++) /*一級蝶形運算*/

	{

		l = 1 << i;

		for
			(j = 0; j<N; j += 2 * l) /*一組蝶形運算*/

		{

			for
				(k = 0; k<l; k++) /*一個蝶形運算*/

			{

				mul(x[j + k + l], W[N*k / 2 / l], &product);

				add(x[j + k], product, &up);

				sub(x[j + k], product, &down);

				x[j + k] = up;

				x[j + k + l] = down;

			}

		}

	}

}





void
output()   /*輸出x結果*/

{

	int
		i;

	printf("\nx傅立葉變換結果\n");

	for
		(i = 0; i<N; i++)

	{

		if
			(i % 4 == 0 && i != 0) printf("\n");

		printf("  %.2f", x[i].real);

		if
			(x[i].img >= 0.0001)

			printf("+%.2fj  ", x[i].img);

		else
		if (fabs(x[i].img)<0.0001)

			printf("+0.0000j  ");

		else

			printf("%.2fj  ", x[i].img);

	}

	printf("\n");

}



/*void save()           //儲存x傅立葉變換結果

{

int i;

ofstream outfile("D:\\result1.txt", ios::out);

if (!outfile)

{

cerr << "open result1.txt erro!" << endl;

exit(1);

}

outfile << "x傅立葉變換結果:" << endl;

for (i = 0; i<N; i++)

{

outfile << x[i].real;

if (x[i].img >= 0.0001)

outfile << "+" << x[i].img << "j" << " ";

else   if (fabs(x[i].img)<0.0001)

outfile << "+" << "0.00" << "j" << " ";

else

outfile << x[i].img << "j" << " ";

}

printf("\n");

outfile.close();

}*/



int
main()

{

	clock_t
		start = clock();

	double
		duration;

	/*對x進行傅立葉變換*/

	int
		i, j = 0, t = 0;

	double
		data[N * 2] = { 0 };

	double
		result[N * 2] = { 0 };

	float
		tx = 0;

	//ofstream outfile1("D:\\result.txt", ios::out);

	for
		(t = 0; t < N; t++)

	{

		data[t] = 5 + 7 * cos(30 * PI / 128 * t - 30 * PI / 180) + 3 * cos(80 * PI / 128 * t - 90 * PI / 180);

	}

	for
		(i = 0; i<N; i++)            //求出x的256個取樣點的值 下一步傅立葉變化

	{

		x[i].real = data[i];

		x[i].img = 0;

		//printf("%.4f ", x[i].real);

		//printf("%.4f ", x[i].img);

	}

	initW(N);

	fftx();

	//output();

	cout << "輸出訊號的模";

	cout << endl;

	for
		(i = 0; i < N; i++)

	{

		result[i] = sqrt(x[i].real*x[i].real + x[i].img*x[i].img);

		if
			(i % 4 == 0)

			cout << endl;

		cout << setprecision(2) << result[i] / N * 2 << " ";

		//outfile1 << result[i] / N * 2 << endl;

	}

	//save();

	clock_t
		end = clock();

	cout << "所用的時間： ";

	duration = (double)(end - start) / CLOCKS_PER_SEC;

	cout << duration << "ms";

	cout << endl;

	return
		0;

}
 

 傅立葉變換，加入一個低通濾波器

// fourie2.cpp : 定義控制檯應用程式的入口點。
//

#include "stdafx.h"
#include <math.h>
#include <malloc.h>
#include <iostream>
#include <istream>
#include <fstream>
#include <sstream>
#include <string>
#include <vector>

#define pi (double) 3.14159265359
#define POWER 9
using namespace std;
/*複數的定義*/
typedef struct
{
	double re;
	double im;
}COMPLEX;
/*複數的加運算*/
COMPLEX Add(COMPLEX c1, COMPLEX c2)
{
	COMPLEX c;
	c.re = c1.re + c2.re;
	c.im = c1.im + c2.im;
	return c;
}
/*負數的減運算*/
COMPLEX Sub(COMPLEX c1, COMPLEX c2)
{
	COMPLEX c;
	c.re = c1.re - c2.re;
	c.im = c1.im - c2.im;
	return c;
}
/*複數的乘運算*/
COMPLEX Mul(COMPLEX c1, COMPLEX c2)
{
	COMPLEX c;
	c.re = c1.re*c2.re - c1.im*c2.im;
	c.im = c1.re*c2.im + c1.im*c2.re;
	return c;
}
/*快速傅立葉變換
TD為時域值，FD為頻域值，power為2的冪數*/
void FFT(COMPLEX *TD, COMPLEX *FD, int power)
{
	int count;
	int i, j, k, bfsize, p;
	double angle;
	COMPLEX *W, *X1, *X2, *X;
	/*計算傅立葉變換點數*/
	count = 1 << power;
	/*分配運算器所需儲存器*/
	W = (COMPLEX *)malloc(sizeof(COMPLEX)*count / 2);
	X1 = (COMPLEX *)malloc(sizeof(COMPLEX)*count);
	X2 = (COMPLEX *)malloc(sizeof(COMPLEX)*count);
	/*計算加權係數*/
	for (i = 0; i<count / 2; i++)
	{
		angle = -i*pi * 2 / count;
		W[i].re = cos(angle);
		W[i].im = sin(angle);
	}
	/*將時域點寫入儲存器*/
	memcpy(X1, TD, sizeof(COMPLEX)*count);
	/*蝶形運算*/
	for (k = 0; k<power; k++)
	{
		for (j = 0; j<1 << k; j++)
		{
			bfsize = 1 << power - k;
			for (i = 0; i<bfsize / 2; i++)
			{
				p = j*bfsize;
				X2[i + p] = Add(X1[i + p], X1[i + p + bfsize / 2]);
				X2[i + p + bfsize / 2] = Mul(Sub(X1[i + p], X1[i + p + bfsize / 2]), W[i*(1 << k)]);
			}
		}
		X = X1;
		X1 = X2;
		X2 = X;
	}
	/*重新排序*/
	cout << "頻域" << endl;
	vector<float> aa;
	for (j = 0; j<count; j++)
	{
		p = 0;
		for (i = 0; i<power; i++)
		{
			if (j&(1 << i))
				p += 1 << power - i - 1;
		}
		FD[j] = X1[p];
		cout << FD[j].re << endl;
		aa.push_back(FD[j].re);
	}
	/*釋放儲存器*/
	free(W);
	free(X1);
	free(X2);
}
/*快速傅立葉反變換，利用快速傅立葉變換
FD為頻域值，TD為時域值，power為2的冪數*/
void IFFT(COMPLEX *FD, COMPLEX *TD, int power)
{
	int i, count;
	COMPLEX *x;
	/*計算傅立葉反變換點數*/
	count = 1 << power;
	/*分配運算所需儲存器*/
	x = (COMPLEX *)malloc(sizeof(COMPLEX)*count);
	/*將頻域點寫入儲存器*/
	memcpy(x, FD, sizeof(COMPLEX)*count);
	/*求頻域點的共軛*/
	for (i = 0; i<count; i++)
	{
		x[i].im = -x[i].im;
	}
	/*呼叫快速傅立葉變換*/
	FFT(x, TD, power);
	/*求時域點的共軛*/
	cout << "反變換" << endl;
	for (i = 0; i<count; i++)
	{
		TD[i].re /= count;
		TD[i].im = -TD[i].im / count;
		cout << TD[i].re << endl;
	}
	/*釋放儲存器*/
	free(x);
}

float min(float a, float b)
{
	return a < b ? a : b;
}

void convolution(double *input1, double *input2, double *output, int mm, int nn)
{
	double *xx = new double[mm + nn - 1];
	double *tempinput2 = new double[mm + nn - 1];
	for (int i = 0; i < nn; i++)
	{
		tempinput2[i] = input2[i];
	}
	for (int i = nn; i < mm + nn - 1; i++)
	{
		tempinput2[i] = 0.0;
	}
	// do convolution   
	for (int i = 0; i < mm + nn - 1; i++)
	{
		xx[i] = 0.0;
		int tem = (min(i, mm)) == mm ? mm - 1 : min(i, mm);
		for (int j = 0; j <= tem; j++)
		{
			double a = (input1[j] * tempinput2[i - j]);
			xx[i] += (input1[j] * tempinput2[i - j]);
		}
	}
	// set value to the output array   
	for (int i = 0; i < mm + nn - 1; i++)
		output[i] = xx[i];
		delete[] xx;
}

void gaussian_filter(COMPLEX *FD,int length,int frequence)
{
	//double *gauss = new double[frequence];
	//double *current_re = new double[frequence];
	//double *data_out = new double[frequence];
	for (int i = 0; i < length;i++)
	{
		if (i>length/2&&i<=length/2+frequence)
		{
			//gauss[i] = exp(-pow(length / 2 - i, 2) / (2 * frequence*frequence));
			double H = exp(-pow(frequence - i, 2) / (2 * (length / 2)*(length / 2)));
			double M = (0.5 + 0.5*cos(2 * pi*(i - frequence) / (length / 2)));
			FD[i].re = FD[i].re*M;
		}

	}
	//convolution(current_re, gauss, data_out, length, length);
	//for (int i = 0; i < length;i++)
	//{
	//	FD[i].re = data_out[i];
	//}
	//delete[] data_out;
	//delete[] current_re;
	//delete[] gauss;
}

void low_pass_filter(string path,int cut_frequence)
{
	ifstream file_in;
	file_in.open(path);//輸入資料路徑
	string lineStr;
	string str;
	int i = 0;
	double dv6[1 << POWER];
	while (getline(file_in, lineStr))
	{
		stringstream ss(lineStr);
		string str;
		int j = 0;
		double b[2];
		while (getline(ss, str, ','))
		{
			double a = atof(str.c_str());
			b[j] = a;
			j++;
		}
		if (i == 1<<POWER) break;

		if (b[1] < -99)
		{
			b[1] = 0;
		}
		dv6[i] = b[1];
		i++;
	}
	file_in.clear();

	const int N = 1 << POWER;
	double	data[N * 2] = { 0 };
	COMPLEX *x = (COMPLEX*)malloc(sizeof(COMPLEX)* N);
	COMPLEX *y = (COMPLEX*)malloc(sizeof(COMPLEX)* N);
	cout << "原始值" << endl;
	for (int t = 0; t < 1 << POWER; t++)
	{
		x[t].re = dv6[t];
		x[t].im = 0;
		cout << x[t].re << endl;
	}
	FFT(x, y, POWER);
	gaussian_filter(y, N, cut_frequence);
	IFFT(y, x, POWER);

	ofstream file_out;
	path = "./T_out.csv";
	file_out.open(path, ios::out);
	float yy[N];
	for (int i = 0; i < 1 << POWER; i++)
	{
		yy[i] = y[i].im;
		file_out << x[i].re << endl;
	}
	file_out.close();
}

int _tmain(int argc, _TCHAR* argv[])
{
	int cut_frequence = 10000;
	string path = "./T.csv";
	low_pass_filter(path, cut_frequence);
	return 0;
}

 

 下面也是一個低通濾波器，速度快一點

#include  <stdio.h>
#include  <math.h>
#include <iostream>
#include <istream>
#include <fstream>
#include <sstream>
#include <string>

using namespace std;
#define LENGTH 512
//w0=40時達到預期效果

/*-----------------------------------------------------------------//
//                 Fast Fourier Transform program                  //
//-----------------------------------------------------------------//
//  Algrithem: Cooley-Tukey Algrithem ( decimation in frequency )  //
//                                                                 //
//   xr  : Real part of original data as the input;                //
//         Real part of the output.                                //
//   xi  : Imaginary part of original data as the input;           //
//         Imaginary part of the output.                           //
//   nr  : must be an positive integer,                            //
//         N=pow(2,nu) is the length of input array (points).      //
//   T   : =1.0 is forward transform;  =-1.0 is reverse transform  //
//                                                                 //
//-----------------------------------------------------------------*/
void fft(float *xr, float *xi, int nr, float T)
{
	int IBIT(int j, int m);
	int   i, j, k, l, n, n2, nr1, i1, j1, k2, k1n2;
	float c, s, p, tr, ti, trc, tic, ars;

	n = (int)(pow(2.0, nr));
	n2 = n / 2;
	nr1 = nr - 1;

	k = 0;
	for (i = 1; i <= nr; i++)
	{
	loop1:  for (j = 1; j <= n2; j++)
	{
				k2 = k / (int)(pow(2.0, nr1));
				p = (float)(IBIT(k2, nr));
				ars = (float)(6.2831852*p / (float)(n));  //輻角
				c = (float)(cos(ars));
				s = (float)(sin(ars));
				k1n2 = k + n2;
				tr = xr[k1n2] * c + xi[k1n2] * s*T;
				ti = xi[k1n2] * c - xr[k1n2] * s*T;
				xr[k1n2] = xr[k] - tr;
				xi[k1n2] = xi[k] - ti;
				xr[k] = xr[k] + tr;
				xi[k] = xi[k] + ti;
				k++;
	}
			k += n2;
			if (k<n)
			{
				goto loop1;
			}
			else
			{
				k = 0;
				nr1 = nr1 - 1;
				n2 /= 2;
			}
	}

	for (j = 1; j <= n; j++)
	{
		i = IBIT(j - 1, nr) + 1;
		if (i>j)
		{
			j1 = j - 1;
			i1 = i - 1;
			trc = xr[j1];
			tic = xi[j1];
			xr[j1] = xr[i1];
			xi[j1] = xi[i1];
			xr[i1] = trc;
			xi[i1] = tic;
		}
	}

	if (T<0.0)
	{
		for (j = 0; j <n; j++)
		{
			xr[j] /= n;
			xi[j] /= n;
		}
	}
}

int IBIT(int j, int m)
{
	int i, it, j1, j2;
	it = 0;
	j1 = j;
	for (i = 1; i <= m; i++)
	{
		j2 = j1 / 2;
		it = it * 2 + (j1 - 2 * j2);
		j1 = j2;
	}
	return it;
}
/**************************************************************/
void read_csv(float *data,int length)
{
	string path = "E:/Programs/MATLAB/T.csv";
	ifstream file_in;
	file_in.open(path);
	string lineStr;
	string str;
	int i = 0;
	while (getline(file_in, lineStr))
	{
		stringstream ss(lineStr);
		string str;
		int j = 0;
		double b[2];

		while (getline(ss, str, ','))
		{
			double a = atof(str.c_str());
			b[j] = a;
			j++;
		}
		if (b[1] < -99)
		{
			b[1] = 0;
		}
		if (i<length)
		{
			data[i] = b[1];
		}
		else
		{
			int a = 0;
		}

		i++;
	}
	file_in.clear();
}

void main()
{
	float xx[LENGTH], yy[LENGTH], h[LENGTH], w0, w1 = 5;
	int i, j;
	w0 = 10;
	read_csv(xx, LENGTH);
	//====================================//
	ofstream file;
	string path = "E:/Programs/MATLAB/T_out.csv";
	file.open(path,ios::out);

	for (i = 0; i<LENGTH; i++)
		yy[i] = 0.0;

	for (i = 0; i<59; i++)
	{
		fft(xx, yy, 9, 1.0);                    //正向快速傅立葉變換求每道資料頻譜
		for (j = 0; j<LENGTH; j++)              //*窗函式
		{
			if (j>w0&&j <= w0 + w1)
			{
				double aa = (0.54 + 0.46*cos(2 * 3.1415927*(j - w0) / w1));
				xx[j] *= (0.54 + 0.46*cos(2 * 3.1415927*(j - w0) / w1));  
				yy[j] = 0.0;
			}
			else
			if (j>w0 + w1)
				yy[j] = xx[j] = 0.0;
		}

		for (j = 1; j<LENGTH; j++)         //頻譜對稱逆向
		{
			xx[LENGTH - j] = xx[j];
			yy[LENGTH - j] = -yy[j];
		}

		fft(xx, yy, 9, -1.0);         //將濾波結果變換到時間與
		//printf("第%d道資料處理成功\n", i + 1);
	}
	for (int i = 0; i < LENGTH;i++)
	{
		file << xx[i] << endl;
	}
}