快速傅里叶变换及其应用

《快速傅里叶变换及其应用》由会员分享，可在线阅读，更多相关《快速傅里叶变换及其应用（9页珍藏版）》请在金锄头文库上搜索。

1、数值与符号计算实验快速傅里叶变换及其应用1实验要求1. 用C+实现复数类，并为其定义必要的运算符.struct Complex double real_; double imag_;Complex(void)；Complex(double constfe real);Complex(double constfe real, double constfe imag);const;const;const;Complex(Complex constfe v);Complex operator+(Complex const& a) Complex operator-(Complex cortstfe

2、a) Complex operator*(Complex constft a) Complex operator/(irtt n) const;；2. void fft(Complex* dst, Complex+ src, int p);快速傅里叶变换. 求复数数组用0才)的傅里叶变换，结果存放在心却乩逻)中.3. void ifft(Complex+ dst, Complex* src, int p);快速傅 里叶逆 变换.求复数数组刖4口挈)的逆傅里叶变换，结果存放在dst(). 2P)中.4. 利用快速傅里叶变换计算长整数乘法.typedef std：:vector Integer;v

3、oid multiply(Integer* rst. Integer constft a, Integer constfe b);2实验原理计算离散傅里叶变换的一种快速算法，简称FFT。快速傅里叶变换是1965 年由J.W.库利和T.W图基提出的。采用这种算法能使计算机计算离散傅里叶变 换所需要的乘法次数大为减少，特别是被变换的抽样点数N越多，FFT算法计 算量的节省就越显著。快速傅氏变换，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅 立叶变换的算法进行改进获得的。它对傅氏变换的理论并没有新的发现，但是对 于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。设x(n)

4、为N项的复数序列，由DFT变换，任一 X(m)的计算都需要N次 复数乘法和N-1次复数加法，而一次复数乘法等于四次实数乘法和两次实数加 法，一次复数加法等于两次实数加法，即使把一次复数乘法和一次复数加法定义 成一次“运算”(四次实数乘法和四次实数加法)，那么求出N项复数序列的X(m)，即N点DFT变换大约就需要N2次运算。当N=1024点甚至更多的时候， 需要N2=1048576次运算.在FFT中，利用WN的周期性和对称性，把一个N项序列(设N=2k, k为 正整数)，分为两个N/2项的子序列，每个N/2点DFT变换需要(N/2) 2次运 算，再用N次运算把两个N/2点的DFT变换组合成一个N

5、点的DFT变换。这 样变换以后，总的运算次数就变成N+2(N/2)2=N+N2/2。继续上面的例子， N=1024时，总的运算次数就变成了 525312次，节省了大约50%的运算量。如果我们将这种“一分为二”的思想不断进行下去，直到分成两两一组的DFT运算单元，那么N点的DFT变换就只需要Nlog2N次的运算，N在1024点 时，运算量仅有10240次，是先前的直接算法的1%，点数越多，运算量的节约 就越大，这就是FFT的优越性。3算法思想及代码实现3.1复数类及常规运算函数定义了复数类，以及复数的“ + ”，“ ”，“*”，“三”四则运算。复数是指能写成如下形式的数a+bi，这里a和b是实数

6、，i是虚数单位(即-1 开根)。在复数a+bi中，a称为复数的实部，b称为复数的虚部，i称为虚数单位。 当虚部等于零时，这个复数就是实数；当虚部不等于零时，这个复数称为虚数， 复数的实部如果等于零，则称为纯虚数。由上可知，复数集包含了实数集，因而 是实数集的扩张。复数是由意大利米兰学者卡当在十六世纪首次引入，经过达 朗贝尔、棣莫弗、欧拉、高斯等人的工作，此概念逐渐为数学家所接受。复数的四则运算规定为：(a+bi) + (c+di) = (a+c) + (b+d) i，(a+bi) (c+di) = (ac) + (b d) i，(a+bi)(c+di) = (acbd) + (bc+ad) i

7、， (c与d不同时为零)。定义复数类a.b.c.d.a.imag-a.imag-Complex:Complex() real_=O;imag_=0; Complex:Complex(double const& real) this-real_=real;this-imag_=O;Complex:Complex(double const& real, double const& imag) this-real_=real;this-imag_=imag; Complex:Complex(Complex const& v) this-real_ = v.real_;this-imag_ = v.i

8、mag_;加法运算Complex Complex:operator+(Complex const& a) const/符号重定义Complex sum;sum.real_ = this-real_ + a.real_; sum.imag_ = this-imag_ + a.imag_; return sum;减法运算Complex Complex:operator-(Complex const & a) const/符号重定义Complex diff;diff.real_ = this-real_ - a.real_;diff.imag_ = this-imag_ - a.imag_;retu

9、rn diff;乘法运算Complex Complex:operator*(Complex const& a) const /符号重定义Complex muti;muti.real_ = this-real_ * a.real_ - this-imag_ *muti.imag_ = this-imag_ * a.real return muti;除法运算+ this-reala.imag_;* a.imag_;Complex Complex:operator/(Complex const& a) const /符号重定义Complex Div;Div.real_ = (this-real_*a

10、.imag_+this-imag_*a.imag_)/(a.real_*a.real_+*a.imag_);Div.imag_=(this-imag_*a.real_-this-real_ *a.imag_);return Div;a.imag_)/(a.real_*a.real_ +3.2复数数组的傅里叶变换傅立叶变换，也称作傅里叶变换，表示能将满足一定条件的某个函数表示成 三角函数(正弦和/或余弦函数)或者它们的积分的线性组合。在不同的研究领 域，傅立叶变换具有多种不同的变体形式，如连续傅立叶变换和离散傅立叶变换。 最初傅立叶分析是作为热过程的解析分析的工具被提出的。f(t)是t的周期函数

11、，如果t满足狄里赫莱条件：在一个周期内具有有限个 间断点，且在这些间断点上，函数是有限值；在一个周期内具有有限个极值点； 绝对可积。则有下图式成立。称为积分运算f(t)的傅里叶变换，式的积分运算叫做F(3 )的傅里叶逆变换。F(3 )叫做f(t)的像函数, f(t)叫做F(3 )的像原函数。F(3 )是坦)的像。f(t)是F(3 )原像。 傅里叶变换 傅里叶逆变换实现复数数组的傅里叶变换函数void fft (Complex* dst, Complex* src, int p) 快速傅里叶变换/求复数数组src0, 2J)的傅里叶变换，结果存放在dst0, 2S)中int pow_p = po

12、w(double)2,p);int flag,m0,ml,m2,m3;Complex *temp,*w,wO;double real,imag;flag = 0;temp = (Complex *)malloc(sizeof(Complex) * pow_p);w = (Complex *)malloc(sizeof(Complex) * (pow_p/2);计算出w;w0.set_Complex(1,0);real = cos(2*PI/pow_p);imag = sin(2*PI/pow_p);wO.set_Complex(real,imag);for (int i = 1; i 0; i

13、/=2)/从二进制的最高位开始处理，依次往低位计算flag+;应该在偶数次得到结果for (int k=0; ki; k+) for (int j=0; j pow_p/(2*i); j+) m0 = ( k + j*i )%(pow_p/2);/ 低位m1 = pow_p/2 + ( k + j*i )%(pow_p ;高位m2 = ( k + j*i*2 )%(pow_p); m3 = ( k + j*i*2 + i)%(pow_p);讦(flag = 1)/第 一次 tempm0 = srcm2 + srcm3*wj*i;tempml = srcm2 - srcm3*wj*i;else

14、f (flag%2 =0)/ 偶数次 dstm0 = tempm2 + tempm3*wj*i;dstm1 = tempm2 - tempm3*wj*i; else/奇数次 tempm0 = dstm2 + dstm3*wj*i;tempm1 = dstm2 - dstm3*wj*i; f (flag%2 = 1)/奇数次时应当把temp中的变换结果放回dest中 for (int i = 0; ipow_p; i+)dsti = tempi; free(temp);free(w);/ Print_Complex_Array(dst, pow_p);3.3复数数组的逆傅里叶变换实现复数数组的逆傅里叶变换函数void ifft(Complex* dst, Complex* src, int p)/ 快速傅里叶逆变换 求复数数组src0, 2Ap)的傅里叶逆变换，结果存放在dst0, 2J)中 int pow_p = pow(double)2,p);/pow_p = 16int flag,m0,m1,m2,m3;Complex *temp,*w,w0;double real,imag;