FFT.cpp
#include "fft.h"
#include <math.h>
//精度0.0001弧度
void conjugate_complex(int n,complex_user in[],complex_user out[])
{
int i = 0;
for(i=0;i<n;i++)
{
out[i].imag = -in[i].imag;
out[i].real = in[i].real;
}
}
void c_abs(complex_user f[],float out[],int n)
{
int i = 0;
float t;
for(i=0;i<n;i++)
{
t = f[i].real * f[i].real + f[i].imag * f[i].imag;
out[i] = sqrt(t);
}
}
void c_plus(complex_user a,complex_user b,complex_user *c)
{
c->real = a.real + b.real;
c->imag = a.imag + b.imag;
}
void c_sub(complex_user a,complex_user b,complex_user *c)
{
c->real = a.real - b.real;
c->imag = a.imag - b.imag;
}
void c_mul(complex_user a,complex_user b,complex_user *c)
{
c->real = a.real * b.real - a.imag * b.imag;
c->imag = a.real * b.imag + a.imag * b.real;
}
void c_div(complex_user a,complex_user b,complex_user *c)
{
c->real = (a.real * b.real + a.imag * b.imag)/(b.real * b.real +b.imag * b.imag);
c->imag = (a.imag * b.real - a.real * b.imag)/(b.real * b.real +b.imag * b.imag);
}
#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
void Wn_i(int n,int i,complex_user *Wn,char flag)
{
Wn->real = cos(2*PI*i/n);
if(flag == 1)
Wn->imag = -sin(2*PI*i/n);
else if(flag == 0)
Wn->imag = -sin(2*PI*i/n);
}
//傅里叶变换
void fft_user(int N,complex_user f[])
{
complex_user t,wn;//中间变量
int i,j,k,m,n,l,r,M;
int la,lb,lc;
/*----计算分解的级数M=log2(N)----*/
for(i=N,M=1;(i=i/2)!=1;M++);
/*----按照倒位序重新排列原信号----*/
for(i=1,j=N/2;i<=N-2;i++)
{
if(i<j)
{
t=f[j];
f[j]=f[i];
f[i]=t;
}
k=N/2;
while(k<=j)
{
j=j-k;
k=k/2;
}
j=j+k;
}
/*----FFT算法----*/
for(m=1;m<=M;m++)
{
la=pow(2,m); //la=2^m代表第m级每个分组所含节点数
lb=la/2; //lb代表第m级每个分组所含碟形单元数
//同时它也表示每个碟形单元上下节点之间的距离
/*----碟形运算----*/
for(l=1;l<=lb;l++)
{
r=(l-1)*pow(2,M-m);
for(n=l-1;n<N-1;n=n+la) //遍历每个分组,分组总数为N/la
{
lc=n+lb; //n,lc分别代表一个碟形单元的上、下节点编号
Wn_i(N,r,&wn,1);//wn=Wnr
c_mul(f[lc],wn,&t);//t = f[lc] * wn复数运算
c_sub(f[n],t,&(f[lc]));//f[lc] = f[n] - f[lc] * Wnr
c_plus(f[n],t,&(f[n]));//f[n] = f[n] + f[lc] * Wnr
}
}
}
}
//傅里叶逆变换
void ifft_user(int N,complex_user f[])
{
int i=0;
conjugate_complex(N,f,f);
fft_user(N,f);
conjugate_complex(N,f,f);
for(i=0;i<N;i++)
{
f[i].imag = (f[i].imag)/N;
f[i].real = (f[i].real)/N;
}
}
FFT.h
#ifndef __FFT_H__
#define __FFT_H__
//
#define PI 3.1415926535897932384626433832795028841971
typedef struct _complex_user
{
double real;
double imag;
}complex_user;
///////////////////////////////////////////
void conjugate_complex(int n, complex_user in[], complex_user out[]);
void c_plus(complex_user a, complex_user b, complex_user *c);//复数加
void c_mul(complex_user a, complex_user b, complex_user *c);//复数乘
void c_sub(complex_user a, complex_user b, complex_user *c); //复数减法
void c_div(complex_user a, complex_user b, complex_user *c); //复数除法
void fft_user(int N, complex_user f[]);//傅立叶变换 输出也存在数组f中
void ifft_user(int N, complex_user f[]); // 傅里叶逆变换
void c_abs(complex_user f[], float out[], int n);//复数数组取模
////////////////////////////////////////////
#endif
main.cpp
#include <iostream>
#include "fft.h"
#define N (256)
#define Fs (128.0)
int main(int argc, char **argv)
{
double f1 = 10;
double p1 = 90/180*PI;
double f2 = 12;
double p2 = 0;
double f3 = 18;
double p3 = 30/180*PI;
//
double t[N] = { 0 };
double f[N] = { 0 };
//
complex_user sigIn[N];
complex_user sigOut[N];
float resOut[N];
int i = 0;
//生成时间序列:t=i*Ts
//对应频率序列:f=i/N*Fs
for (i = 0; i < N; i++)
{
t[i] = i / Fs;
}
for (i = 0; i < N; i++)
{
f[i]= i*Fs/N;
}
//生成输入信号
for (i = 0; i < N; i++)
{
sigIn[i].real = 5 + 2 * sin(2 * PI*f1*t[i] + p1) + 8 * sin(2 * PI*f2*t[i] + p2) + 3 * sin(2 * PI*f3*t[i] + p3);
sigIn[i].imag = 0;
//
sigOut[i].real = sigIn[i].real;
sigOut[i].imag = 0;
}
//
//FFT
fft_user(N, sigOut);
//求绝对值
c_abs(sigOut, resOut, N);
for (i = 0; i < N; i++)
{
resOut[i] = resOut[i] / N * 2;
if (i == 0)
{
resOut[i] = resOut[i] / 2;
}
//绝对值小于0.01则过滤
if (resOut[i] < 0.01)
{
resOut[i] = 0;
sigOut[i].real = 0;
sigOut[i].imag = 0;
}
}
//IFFT
std::cin.clear();
std::getchar();
return 0;
}