フーリエ変換C++モード
9679 ワード
フーリエ変換と逆変換の3組のコード
フーリエ変換、ローパスフィルタを追加
次もローパスフィルタで、速度が少し速いです
#include
#include
#include
#include
#include
#include
#include
#include
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; i0)
{
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= 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("
");
}
/*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= 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("
");
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 フーリエ変換、ローパスフィルタを追加
// fourie2.cpp : 。
//
#include "stdafx.h"
#include
#include
#include
#include
#include
#include
#include
#include
#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 aa;
for (j = 0; jlength/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< 次もローパスフィルタで、速度が少し速いです
#include
#include
#include
#include
#include
#include
#include
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 (kj)
{
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 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