复数基2 DIT FFT程序

好象不断有人问关于FFT的问题，我也发一个以前写的FFT程序吧。先发一个复数运算的，过几天再发一个定点的。我自己觉得这个程序的可读性还是比较好的，如果你能理解蝶形运算，你就应该能理解这个程序。我的程序没有用通用的STL复数类，而用了我以前自己写的一个复数运算类，只是做了一些STL化，大学里的老师应该不会骂我了吧:)。我的这个复数类也有其优点，可以支持直角坐标和极坐标的交互运算。比如说，假设第一个数c1是极坐标（1， PI/2），第二个数c2是笛卡尔坐标（3， 4i），你可以不用考虑转换问题直接运算c1（+，-，*，/）c2，运算结果的坐标类型由第一个操作数决定。

好了，看程序吧。需要指出的是，为了方便，我把复数类的定义和实现写在一个文件里了。

// FILE: my_complex.h

#ifndef _MY_COMPLEX_H_
#define _MY_COMPLEX_H_

#include <cmath>
#include <cassert>
#include <iostream>

static const double PI = 3.1415926536;
static const double epsilon=1e-8;

class Complex
{
double real;
double img;
bool is_polar;

#define IsZero(x) ( (bool)((x)<epsilon && (x)>(-epsilon)) )

void check_zero(void);
public:
Complex() : real(0.0), img(0.0), is_polar(false) {};

//no explicit here, otherwise you cannot use codes like Complex X=32; etc.
Complex(double r) : real(r), img(0.0), is_polar(false) { check_zero(); };

Complex(double r, double i, bool polar=false) : real(r), img(i), is_polar(polar)
{ if(is_polar) assert(real>(-epsilon)); check_zero(); };

inline const double get_real() const;
inline const double get_img() const;
inline const double get_rho() const;
inline const double get_theta() const;
inline const bool IsPolar() const { return is_polar; };

void to_polar(void);
void to_ri(void);
void conj(void) { img = -img; };

const Complex& operator = (double a);
const Complex& operator = (const Complex& a);
const bool operator == (const Complex& a) const;

const Complex operator -() const;

const Complex operator +(const Complex& a) const;
const Complex operator -(const Complex& a) const;
const Complex operator *(double a) const;
const Complex operator *(const Complex& a) const;
const Complex operator /(double a) const;
const Complex operator /(const Complex& a) const;

friend std::ostream& operator <<(std::ostream& os, const Complex& a);

friend const Complex conj(const Complex& a);
};

void Complex::check_zero(void)
{
if(IsZero(real))
real = 0.0;
if(IsZero(img))
img = 0.0;
}

void Complex::to_polar(void)
{
if(is_polar==false)
{
  double temp_real = get_rho();
  img = get_theta();
  real = temp_real;
  is_polar = true;
}
}

void Complex::to_ri(void)
{
if(is_polar)
{
  double temp_real = get_real();
  img = get_img();
  real = temp_real;
  is_polar = false;
}
}

const double Complex::get_real() const
{
return (is_polar ? real*::cos(img):real);
}

const double Complex::get_img() const
{
return (is_polar ? real*::sin(img):img);
}

const double Complex::get_rho() const
{
return (is_polar ? real:std::sqrt(real*real+img*img));
}

const double Complex::get_theta() const
{
if(is_polar)
return img;
else
return ::atan2(img, real);
}

inline const Complex& Complex::operator = (double a)
{
is_polar = false;
real = a;
img = 0;

return *this;
}

inline const Complex& Complex::operator = (const Complex& a)
{
is_polar = a.IsPolar();
real = (is_polar ? a.get_rho():a.get_real());
img = (is_polar ? a.get_theta():a.get_img());

return *this;
}

inline const bool Complex::operator == (const Complex& a) const
{
if(is_polar == a.is_polar)
return (IsZero(real-a.real) && IsZero(img-a.img));
else
return (IsZero(get_real()-a.get_real()) && IsZero(get_img()-a.get_img()));
}

inline const Complex Complex::operator -() const
{
return (is_polar ? Complex(real, PI+img, true):Complex(-real, -img));
}

const Complex Complex::operator +(const Complex& a) const
{
Complex temp(get_real()+a.get_real(), get_img()+a.get_img());
if(is_polar)
temp.to_polar();
return temp;
}

const Complex Complex::operator -(const Complex& a) const
{
Complex temp(get_real()-a.get_real(), get_img()-a.get_img());
if(is_polar)
temp.to_polar();
return temp;
}

const Complex Complex::operator *(double a) const
{
assert(!IsZero(a));

if(is_polar)
return Complex(get_rho()*a, get_theta(), true);
else
return Complex(get_real()*a, get_img()*a);
}

const Complex Complex::operator *(const Complex& a) const
{
if(is_polar)
return Complex(get_rho()*a.get_rho(), get_theta()+a.get_theta(), true);
else
return Complex(get_real()*a.get_real()-get_img()*a.get_img(), get_real()*a.get_img()+get_img()*a.get_real());
}

const Complex Complex::operator /(double a) const
{
assert(!IsZero(a));

if(is_polar)
return Complex(get_rho()/a, get_theta(), true);
else
return Complex(get_real()/a, get_img()/a);
}

const Complex Complex::operator /(const Complex& a) const
{
Complex temp(get_rho()/a.get_rho(), get_theta()-a.get_theta(), true);
if(is_polar==false)
temp.to_ri();
return temp;
}

std::ostream& operator <<(std::ostream& os, const Complex& a)
{
if(a.IsPolar()==false)
os<<'('<<a.get_real()<<"r,"<<a.get_img()<<"i)";
else
os<<'<'<<a.get_rho()<<"h,"<<a.get_theta()<<"t)";
return os;
}

inline const Complex conj(const Complex& a)
{
return Complex(a.real, -a.img, a.is_polar);
}

#endif

// FILE: fft_dit_base2.cpp

#include "my_complex.h"

using namespace std;

static const int N=32;

static Complex fft_coeff[N];
static Complex ifft_coeff[N];

void setup_table(void)
{
for(int i=0; i<N; ++i)
{
fft_coeff[i] = Complex(1.0, -PI*i*2/N, true);
ifft_coeff[i] = Complex(1.0, PI*i*2/N, true);
}
}

void decimation(Complex data[])
{
int r=0;
int i=1;
#define swp(x, y) {Complex temp=(x); (x)=(y); (y)=temp;}

while(i<N/2)
{
  r += N/2;
  swp(data[i], data[r]);
  ++i;

for(int m=N>>1; !((r^=m)&m); m>>=1);
  if(r>i)
  {
   swp(data[i], data[r]);
   swp(data[N-1-i], data[N-1-r]);
  }
  ++i;
}
#undef swp
}

void fft_dit(Complex data[], bool ifft=false)
{
decimation(data);

Complex *coeff;
if(ifft)
coeff = ifft_coeff;
else
coeff = fft_coeff;

int groups=N/2;
int step = 2;
while(groups>=1)
{
  for(int i=0; i<groups; ++i)
  {
   int first = i*step;
   int second = first+step/2;
   int coeff_idx = 0;
   for(int j=0; j<step/2; ++j)
   {
    Complex temp = data[second]*coeff[coeff_idx];

data[second] = data[first] - temp;
data[first] = data[first] + temp;

if(ifft==true)
    {
     data[first] = data[first]/2;
     data[second] = data[second]/2;
    }

++first;
    ++second;
    coeff_idx += groups;
   }
  }
  groups /= 2;
  step *= 2;
}
}

int main()
{
setup_table();

Complex raw_data[N];
for(int i=0; i<N; ++i)
raw_data[i] = i+1;

cout<<"The data to be fft-ed are:"<<endl;
for(int i=0; i<N; ++i)
cout<<raw_data[i]<<' ';
cout<<endl;

fft_dit(raw_data);
cout<<"The fft result are:"<<endl;
for(int i=0; i<N; ++i)
cout<<raw_data[i]<<' ';
cout<<endl;

fft_dit(raw_data, true);
cout<<"The ifft-ed data are:"<<endl;
for(int i=0; i<N; ++i)
cout<<raw_data[i]<<' ';
cout<<endl;

system("PAUSE");
return 0;
}

参考资料:

1. http://www.jjj.de/fxt/, 极具参考价值的网站

2. 任何一本讲信号与系统,或数字信号处理的书

复数基2 DIT FFT程序相关推荐

【 FPGA 】16点并行DIT FFT的实现
目录整体架构介绍旋转因子介绍代码文件结构重点难点易错点整体架构介绍 16点并行FFT分为4级蝶形运算,每一级蝶形运算有一个基本的蝶形单元: 如下是16点DIT FFT的数据流图: 可见,第0 ...
对比DFT程序与FFT程序的效率
徐士良老师编写的c语言算法程序下载链接:https://pan.baidu.com/s/1zDV6iLeYeXmZaoZlP4yRAA 提取码:8opo 一:徐士良老师编写的FFT程序比较 1.FFT ...
请用python写个fft程序
下面是用Python写的一个快速傅里叶变换(FFT)程序:import numpy as np def fft(x): """A recursive implementa ...
基2FFT算法matlab程序编写,频率抽取（DIF）基2FFT算法的MATLAB实现
频率抽取(DIF)基2FFT算法和时间抽取(DIT)基2FFT算法是两种等价的FFT算法,其相同之处: (1)DIF与DIT两种算法均为原位运算. (2)DIF与DIT运算量相同. 不同之处: (1) ...
基2FFT算法matlab程序编写,按时间抽取的基2FFT算法分析及MATLAB实现
按时间抽取的基2FFT 算法分析及MATLAB 实现 1 DIT-FFT 算法的基本原理有限长序列x (n )的N 点DFT 定义为:∑-==10 )()(N n n k N W n x k X , ...
c语言实现1024点fft程序,数字信号处理的步骤与注意事项,并编写1024个采样点的FFT C语言程序...
数字信号处理的步骤与注意事项,并编写1024个采样点的FFT C语言程序 1. 数字信号处理 1.1 数字信号处理概述数字信号处理是研究如何用数字或符号序列来表示信号以及如何对这些序列进行处理的一门 ...
复数矩阵求逆的 C 语言程序
关于复数矩阵求逆,如果使用 MATLAB,就非常简单.我们先用一个 MATLAB 的例子来说明,等会要将 C 语言的程序和 MATLAB 的程序进行对比. close all; clear all; ...
c语言实现1024点fft程序,C语言1024点快速傅里叶变换（FFT）程序，最好经过优化，执行速度快...
voidfft(){intnn,n1,n2,i,j,k,l,m,s,l1;floatar[1024],ai[1024];//实部虚部floata[2050];floatt1,t2,x,y;floatw ...
通信原理及系统系列34——基2-N点FFT蝶形运算结构推演分析
gsl科学计算库文档，翻译了索引，凑合看看。
1.介绍 2.库的使用本章描述如何编译使用GSL的程序,介绍GSL的一般用法. 2.1例子程序 2.2编译和链接 2.3共享库 2.4与ANSI C兼容性 2.5inline函数 2.6长双精度lo ...

复数基2 DIT FFT程序

复数基2 DIT FFT程序相关推荐

最新文章

热门文章