三组傅里叶变换和反变换的代码
#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;
}
}