基-2FFT算法的软件实现
- 格式:doc
- 大小:4.62 MB
- 文档页数:13
基2FFT算法一、什么是基2FFT算法基2FFT算法(Fast Fourier Transform)是一种用于将一个离散信号变换成其频域表示的算法。
它通过将离散信号分解为多个较小的离散信号,并对这些信号进行递归处理,最后合并得到频域表示。
基2FFT算法是最常用的FFT算法之一,其时间复杂度为O(nlogn)。
二、基2FFT算法的原理基2FFT算法的核心思想是分治法。
根据离散信号的长度,将信号分解为两个较小的子信号,然后对这两个子信号分别进行FFT变换。
再将得到的两个子信号的频域表示合并,得到整个信号的频域表示。
具体过程如下:1.判断信号长度是否为1,如果是,则直接返回该信号作为其频域表示。
2.将信号分为偶数索引和奇数索引两个子信号。
3.对这两个子信号分别进行FFT变换,并得到它们的频域表示。
4.合并两个子信号的频域表示,得到整个信号的频域表示。
基于分治法,基2FFT算法可以将信号的长度减半,并通过递归的方式处理子信号,直到信号长度为1,再进行合并得到整个频域表示。
三、基2FFT算法的实现下面是基于伪代码的基2FFT算法的实现:function FFT(signal):n = length(signal)if n == 1:return signalelse:even = FFT(signal[0::2])odd = FFT(signal[1::2])for k = 0 to n/2-1:t = odd[k] * exp(-2*pi*i*k/n)signal[k] = even[k] + tsignal[k + n/2] = even[k] - treturn signal四、基2FFT算法的应用基2FFT算法在信号处理、图像处理、数字滤波器等领域具有广泛的应用。
以下是一些基于基2FFT算法的应用示例:1.频谱分析:将时间域信号转换为频域信号,可以分析信号的频率分量及其强弱关系,从而实现信号的频谱分析。
离散时间信号的基2快速傅里叶变换FFT (时间抽取)蝶形算法实现一、一维连续信号的傅里叶变换连续函数f(x)满足Dirichlet (狄利克雷)条件下,存在积分变换:正变换:2()()()()j ux F u f x e dx R u jI u π+∞--∞==+⎰ 反变换:2()()j ux f x F u e du π+∞-∞=⎰其中()()cos(2)R u f t ut dt π+∞-∞=⎰,()()sin(2)I u f t ut dt π+∞-∞=-⎰定义幅值、相位和能量如下:幅度:1222()()()F u R u I u ⎡⎤⎡⎤=+⎣⎦⎣⎦ 相位:()arctan(()/())u I u R u ϕ= 能量:22()()(E u R u I u =+)二、一维离散信号的傅里叶变换将连续信号对自变量进行抽样得到离散信号(理想冲击抽样脉冲),利用连续信号的傅里叶变换公式得到离散时间傅里叶变换DTFT ;再利用周期性进行频域抽样,得离散傅里叶变换DFT (详情参考任何一本《数字信号处理》教材)。
DFT 变换如下:正变换:12/0()(),0,1,2,1N j ux Nx F u f x eu N π--===-∑。
反变换:12/01()(),0,1,2,1N j ux Nu f x F u ex N Nπ-===-∑。
DFT 是信号分析与处理中的一种重要变换,因为计算机等数字设备只能存储和处理离散数据(时域、频域)。
因直接计算DFT 的计算量大(与变换区间长度N 的平方成正比,当N 较大时,计算量太大),所以在快速傅里叶变换(简称FFT)出现以前,直接用DFT 算法进行谱分析和信号的实时处理是不切实际的。
直到1965年发现了DFT 的一种快速算法(快速傅里叶变换,即FFT )以后,情况才发生了根本的变化。
FFT 有时间抽取和频率抽取两种,下面介绍时间抽取FFT 。
三、时间抽取的基2快速傅里叶变换FFT令2j NN W eπ-=,则2jkm km NNWeπ-=称为旋转因子,把DFT 正变换改写为:1[][],0,1,1N km N k F m f k W m N -===-∑将函数记作x ,变换后为X ,则为:10[][],0,1,1N kmN k X m x k W m N -===-∑时间抽取的FFT 算法利用了旋转因子的三个性质:周期性、对称性和可约性。
华北电力大学实验报告||实验名称FFT的软件实现实验(Matlab)IIR数字滤波器的设计课程名称信号分析与处理||专业班级:电气化1308 学生姓名:袁拉麻加学号: 2 成绩:指导教师:杨光实验日期: 2015-12-17快速傅里叶变换实验一、实验目的及要求通过编写程序,深入理解快速傅里叶变换算法(FFT)的含义,完成FFT和IFFT算法的软件实现。
二、实验内容利用时间抽取算法,编写基2点的快速傅立叶变换(FFT)程序;并在FFT程序基础上编写快速傅里叶反变换(IFFT)的程序。
三:实验要求1、FFT和IFFT子程序相对独立、具有一般性,并加详细注释;2、验证例6-4,并能得到正确结果。
3、理解应用离散傅里叶变换(DFT)分析连续时间信号频谱的数学物理基础。
四、实验原理:a.算法原理1、程序输入序列的元素数目必须为2的整数次幂,即N=2M,整个运算需要M 级蝶形运算;2、输入序列应该按二进制的码位倒置排列,输出序列按自然序列排列;3、每个蝶形运算的输出数据军官占用其他输入数据的存储单元,实现“即位运算”;4、每一级包括N/2个基本蝶形运算,共有M*N/2个基本蝶形运算;5、第L级中有N/2L个群,群与群的间隔为2L。
6、处于同一级的各个群的系数W分布相同,第L级的群中有2L-1个系数;7、处于第L级的群的系数是(p=1,2,3,…….,2L-1)而对于第L级的蝶形运算,两个输入数据的间隔为2L-1。
b.码位倒置程序流程图开始检测A序列长度nk=0j=1x1(j)=bitget(k,j);j=j+1Yj<m?Nx1=num2str(x1);y(k+1)=bin2dec(x1);clear x1k=k+1c.蝶形运算程序流程图五、程序代码与实验结果a.FFT程序:%%clear all;close all;clc;%输入数据%A=input('输入x(n)序列','s');A=str2num(A);% A=[1,2,-1,4]; %测试数据%%%%校验序列,%n=length(A);m=log2(n);if (fix(m)~=m)disp('输入序列长度错误,请重新输入!');A=input('输入x(n)序列','s');A=str2num(A);elsedisp('输入正确,请运行下一步')end%%%码位倒置%for k=0:n-1for j=1:m %取M位的二进制数%x1(j)=bitget(k,j); %倒取出二进制数%endx1=num2str(x1); %将数字序列转化为字符串%y(k+1)=bin2dec(x1); %二进制序列转化为十进制数%clear x1endfor k=1:nB(k)=A(y(k)+1); %时间抽取序列%endclear A%%%计算%for L=1:m %分解为M级进行运算%LE=2^L; %第L级群间隔为2^L%LE1=2^(L-1); %第L级中共有2^(L-1)个Wn乘数,进行运算蝶运算的两数序号相隔LE1%W=1;W1=exp(-1i*pi/LE1);for R=1:LE1 %针对第R个Wn系数进行一轮蝶运算,共进行LE1次%for P=R:LE:n %每个蝶的大小为LE% Q=P+LE1;T=B(Q)*W;B(Q)=B(P)-T;B(P)=B(P)+T;endW=W*W1;endendB %输出X(k)%%%验证结果:例6-4b.IFFT程序:%%clear all;close all;clc;%输入数据%A=input('输入X(k)序列','s');A=str2num(A);% A=[6,2+2i,-6,2-2i]; %测试数据%%%%校验序列,%n=length(A);m=log2(n);if (fix(m)~=m)disp('输入序列长度错误,请重新输入!');A=input('输入x(n)序列','s');A=str2num(A);elsedisp('输入正确,请运行下一步')end%%%码位倒置%for k=0:n-1for j=1:m %取M位的二进制数%x1(j)=bitget(k,j); %倒取出二进制数%endx1=num2str(x1); %将数字序列转化为字符串%y(k+1)=bin2dec(x1); %二进制序列转化为十进制数%clear x1endfor k=1:nB(k)=A(y(k)+1); %时间抽取序列%endclear A%%%计算%for L=1:m %分解为M级进行运算%LE=2^L; %第L级群间隔为2^L%LE1=2^(L-1); %第L级中共有2^(L-1)个Wn乘数,进行运算蝶运算的两数序号相隔LE1%W=1;W1=exp(-1i*pi/LE1);for R=1:LE1 %针对第R个Wn系数进行一轮蝶运算,共进行LE1次%for P=R:LE:n %每个蝶的大小为LE%Q=P+LE1;T=B(Q)*W;B(Q)=B(P)-T;B(P)=B(P)+T;endW=W*W1;endendB=conj(B); %取共轭%B=B/n %输出x(n)%验证结果:六、实验心得与结论本次实验借助于Matlab软件,我避开了用C平台进行复杂的复数运算,在一定程度上简化了程序,并添加了简单的检错代码,码位倒置我通过查阅资料,使用了一些函数,涉及到十-二进制转换,数字-文本转换,二-文本转换,相对较复杂,蝶运算我参考了书上了流程图,做些许改动就能直接实现。
一、概述Fast Fourier Transform(快速傅里叶变换,FFT)是一种非常重要的数学算法,它在信号处理、图像处理、通信等领域有着广泛的应用。
其中,radix-2 FFT算法是FFT算法中最为简单、高效的一种算法。
本文将对radix-2 FFT算法的原理进行介绍,并通过具体的例子来演示其工作原理。
二、FFT算法简介FFT算法是一种离散傅里叶变换的快速计算方法,它能够将时域的信号转换为频域的信号,从而方便地进行频域分析和处理。
在数字信号处理中,FFT算法广泛应用于滤波、频谱分析、信号合成等方面。
FFT 算法的应用十分广泛,因此对其进行深入的学习和理解具有重要意义。
三、radix-2 FFT算法原理1. 蝶形运算:radix-2 FFT算法采用了蝶形运算(butterfly operation)的思想。
蝶形运算是一种迭代的运算方式,它将输入信号分为偶数点和奇数点,然后通过相位的旋转实现频域变换。
具体而言,蝶形运算可以表示为:\[ X(k) = X_e(k) + W_N^kX_o(k), \quad 其中,W_N^k为旋转因子,k为频率索引,X_e(k)和X_o(k)分别为偶数点和奇数点的频谱。
蝶形运算是FFT算法中的核心操作,它通过不断的迭代运算,将时域信号转换为频域信号。
2. 蝶形运算的分解:在radix-2 FFT算法中,蝶形运算会被分解为多个小规模的蝶形运算,这样可以大大简化运算的复杂度。
具体而言,对于长度为N的信号,radix-2 FFT算法将其分解为长度为2的小块,然后通过不断的迭代运算,完成整个频域变换过程。
这种分解的方式大大降低了计算的时间复杂度,使得FFT算法能够高效地完成频域变换。
3. 基-2的选择:radix-2 FFT算法之所以被称为radix-2算法,是因为它将信号的长度N限定为2的幂次。
这样一来,在进行频域变换时,可以将信号不断拆分为长度为2的小块,然后通过不断的迭代运算完成整个频域变换。
按时间抽选的基2FFT算法基2FFT算法(Fast Fourier Transform,快速傅里叶变换)是一种高效的算法,用于在计算机上计算离散傅里叶变换(Discrete Fourier Transform,DFT)。
它的核心思想是利用分治策略和递归操作,在O(nlogn)的时间复杂度下完成离散傅里叶变换。
基2FFT算法的关键步骤如下:1. 将输入的序列划分为两个子序列:偶数位置和奇数位置上的元素分别组成两个子序列。
2. 对这两个子序列分别进行离散傅里叶变换,得到两个新的子序列。
3. 将两个新子序列的元素按照原始顺序交替排列,得到最终的结果。
基于以上步骤,可以利用递归操作来实现基2FFT算法。
具体的实现过程如下:1. 如果输入序列的长度为1,则不需要进行任何操作,直接返回该序列作为结果。
2. 如果输入序列的长度大于1,则按照上述步骤进行分割和计算。
3. 首先将输入序列分为偶数位置和奇数位置上的元素组成的两个子序列。
4. 对这两个子序列分别递归调用基2FFT算法,得到两个新的子序列。
5. 将两个新子序列的元素按照原始顺序交替排列,得到最终的结果。
基2FFT算法的时间复杂度分析如下:假设输入序列的长度为n,则每一层递归的时间复杂度为O(n),总共有logn层递归。
因此,基2FFT算法的总时间复杂度为O(nlogn)。
基2FFT算法在信号处理、图像处理等领域具有广泛的应用。
它可以高效地计算离散傅里叶变换,将时域信号转换为频域信号,从而实现信号分析、频谱分析和频域滤波等操作。
同时,基于基2FFT算法的快速傅里叶变换还能够应用于多项式乘法、高效计算卷积等问题的求解。
总之,基2FFT算法是一种高效的离散傅里叶变换算法,通过利用分治策略和递归操作,能够在O(nlogn)的时间复杂度下完成计算。
它在信号处理和图像处理等领域有着重要的应用价值。
基2FFT算法(Fast Fourier Transform,快速傅里叶变换)是离散傅里叶变换(Discrete Fourier Transform,DFT)的一种高效计算方法。
按时间抽取的基2FFT算法分析及MATLAB实现基2FFT算法是一种快速傅里叶变换(Fast Fourier Transform,FFT)的算法,在信号处理、图像处理等领域有着广泛的应用。
该算法通过将N个输入值分解成两个长度为N/2的DFT(离散傅里叶变换)来实现快速的计算。
本文将对基2FFT算法进行分析,并给出MATLAB实现。
基2FFT算法的主要思路是将输入序列分解成奇偶两个子序列,然后分别对这两个子序列进行计算。
具体步骤如下:1.将输入序列拆分成奇数位和偶数位两个子序列。
比如序列x[0],x[1],x[2],x[3]可以拆分成x[0],x[2]和x[1],x[3]两个子序列。
2. 对两个子序列分别进行DFT计算。
DFT的定义为:X[k] = Σ(x[n] * exp(-i * 2π * k * n / N)),其中k为频率的索引,N为序列长度。
3.对得到的两个DFT结果分别进行合并。
将奇数位子序列的DFT结果和偶数位子序列的DFT结果合并成一个长度为N的DFT结果。
4.重复以上步骤,直到计算结束。
基2FFT算法的时间复杂度为O(NlogN),远远小于直接计算DFT的时间复杂度O(N^2)。
这是因为基2FFT算法将问题的规模逐步减半,从而实现了快速的计算。
下面是MATLAB中基2FFT算法的简单实现:```matlabfunction X = myFFT(x)N = length(x);if N == 1X=x;%递归结束条件return;endeven = myFFT(x(1:2:N)); % 偶数位子序列的FFT计算odd = myFFT(x(2:2:N)); % 奇数位子序列的FFT计算W = exp(-1i * 2 * pi / N * (0:N/2-1)); % 蝶形因子temp = W .* odd; % 奇数位子序列的DFT结果乘以蝶形因子X = [even + temp, even - temp]; % 合并得到一个长度为N的DFT结果end```上述代码中,函数myFFT为基2FFT算法的MATLAB实现。
FFT算法研究及基2-FFT算法的C语⾔实现毕业设计 [论⽂]题⽬:FFT算法研究及基2-FFT算法的C语⾔实现学院:电⽓与信息⼯程学院专业:电⽓⼯程及其⾃动化姓名:XXX学号:XXXXXX指导⽼师:XXX完成时间:2015年06⽉01⽇摘要离散傅⽴叶变换(DFT)常常⽤于计算信号处理。
DFT算法可以得到信号的频域特性,因为该算法在计算上是密集的,长时间的使⽤时,计算机不能实时进⾏信号处理。
所以DFT被发现之后的相当长时间内是没被应⽤到实际的项⽬。
到了⼆⼗世纪六⼗年代中期⼀种新的计算⽅法被研究者发现出来,它就是FFT。
FFT 并不是⼀种新的获取频域特征的⽅式,⽽是离散傅⾥叶变换的⼀种快速实现算法。
数字信号处理在当今科技发展中发展很迅速,不但是在传统的通信领域,其他⽅⾯也经常⽤到。
利⽤快速傅⾥叶变换,实现了信号频域的变换处理。
对于信号的处理,往往会和数学中的算法联系到⼀起。
如果处理得当,将会对⽓象,地理信息等的发展,起到举⾜轻重的作⽤,同时对世界其他领域的发展有很⼤的促进作⽤。
关键词: FFT算法,C语⾔,编译实现AbstractDiscrete Fourier Transform (DFT) is often used to calculate the signal processing to obtain frequency domain signals. DFT algorithm can get the frequency domain characteristics of the signal, because the algorithm is computationally intensive, long-time use, the computer is not conducive to real-time signal processing. So DFT since it was discovered in a relatively long period of time is not to be applied to the actual projects until a new fast discrete Fourier calculation method --FFT is found in discrete Fourier transform was able to actually project has been widely used. FFT is not a new way to get the frequency domain, but the discrete Fourier transform of a fast algorithm.Fast Fourier Transform (FFT) is a digital signal processing important tool that the time domain signal into a frequency-domain signal processing. matched filtering has important applications. FFT is a discrete Fourier transform (DFT) is a fast algorithm, it can be a signal from the time domain to the frequency domain. Some signals in the time domain is not easy to observe the characteristics of what is, but then if you change the signal from the time domain to the frequency domain, it is easy to see out of. This design is required to be familiar with the basic principles of FFT algorithm, based on the preparation of C language program to achieve FFT algorithm real number sequence. Keywords: FFT algorithm, C language compiler to achieve⽬录摘要 ................................................................................................................................. I Abstract ............................................................................................................................. I I ⽬录............................................................................................................................. I II 1 引⾔ . (4)1.1 课题背景 (4)1.2 FFT算法的现状 (4)1.3 研究内容 (2)1.4 论⽂的研究成果 (2)2 数字信号处理综述 (3)2.1 数字信号理论 (3)2.2 数字信号处理的实现 (4)2.3 数字信号处理的应⽤及特点 (4)3 基本理论 (6)3.1 FFT算法的基本概念 (6)3.1.1离散傅⾥叶变换(DFT) (6)3.1.2 快速傅⾥叶变换(FFT) (7)3.2 FFT算法的分类 (8)3.2.1按时间抽取(DIT)的FTT (8)3.2.2 按频率抽取(DIF)的FTT (12)3.2.3 快速傅⾥叶分裂基FFT算法 (15)3.2.4 N为组合数的FFT——混合基算法 (18)3.3 傅⾥叶变换的应⽤ (20)4 基-2FFT算法的C语⾔实现及仿真 ........................................ 错误!未定义书签。
2维FFT算法实现——基于GPU的基2快速⼆维傅⾥叶变换上篇讲述了⼀维FFT的GPU实现(),后来我⼜由于需要做了⼀下⼆维FFT,⼤概思路如下。
⾸先看的肯定是公式:如上⾯公式所描述的,2维FFT只需要拆分成⾏FFT,和列FFT就⾏了,其中我在下⾯的实现是假设原点在F(0,0),由于我的代码需要原点在中⼼,所以在最后我将原点移动到了中⼼。
下⾯是原点F(0,0)的2维FFT的伪代码://C2DFFT//被执⾏2DFFT的是⼀个N*N的矩阵,在source_2d中按⾏顺序储存//⽔平⽅向FFTfor (int i=0;i<N;i++){fft1(&source_2d[i*N],&source_2d_1[i*N],N);}//转置列成⾏for (int i=0;i<N*N;i++){int x = i%N;int y = i/N;int index = x*N+y;source_2d[index] = source_2d_1[i];}//垂直FFTfor(int i=0;i<N;i++){fft1(&source_2d[i*N],&source_2d_1[i*N],N);}//转置回来for (int i=0;i<N*N;i++){int x = i%N;int y = i/N;int index = x*N+y;source_2d[index] = source_2d_1[i];}GPU实现⽆⾮把这些东西转换到GPU上。
我基于OpenGL的fragment shader来计算fft;数据都存放在纹理或者FBO⾥⾯。
和1维fft不同的是,NXN的数据⾥⾯,只是对当前列或者当前排做⼀维FFT,所以bit反转表只需要⼀个1*N的buffer就可以了。
对应的蝴蝶图数据也只需要1*N即可。
所以我们有如下的分配:static ofFbo _fbo_bitrev_table;static ofFbo _origin_butterfly_2d;_fbo_bitrev_table.allocate(N,1,GL_RGBA32F);_origin_butterfly_2d.allocate(N,1,GL_RGBA32F);⾸先要做的是把长度为N的bit反转表求出来,这个只需要求⼀次,所以在最开始的时候就⽤CPU求出来:for(int i=0;i<N;i++){_bitrev_index_2d.setColor(i,0,ofFloatColor(bit_rev(i,N-1),0,0,0));}_bitrev_index_2d.update();//翻转后的索引_fbo_bitrev_table.begin();_bitrev_index_2d.draw(0,0,N,1);_fbo_bitrev_table.end();然后初始化最初的蝴蝶图,这个和1维FFT是⼀样的,只是长度不同⽽已:for(int i=0;i<N;i++){//初始化⼆维蝴蝶图if(i%2==0){_data_2d.setColor(i,0,ofFloatColor(0.f,2.f,0,i+1));}else{_data_2d.setColor(i,0,ofFloatColor(1.f,2.f,0,i-1));}}_data_2d.update();/////////////////2D初始化///////////////////初始化2D蝴蝶图_weight_index_2d[0].begin();_data_2d.draw(0,0,N,1);_weight_index_2d[0].end();//备份2D初始蝴蝶图,⽤于下⼀次新的计算_origin_butterfly_2d.begin();_data_2d.draw(0,0,N,1);_origin_butterfly_2d.end();辅助函数:static unsigned int bit_rev(unsigned int v, unsigned int maxv){unsigned int t = log(maxv + 1)/log(2);unsigned int ret = 0;unsigned int s = 0x80000000>>(31);for (unsigned int i = 0; i < t; ++i){unsigned int r = v&(s << i);ret |= (r << (t-i-1)) >> (i);}return ret;}static void bit_reverse_copy(RBVector2 src[], RBVector2 des[], int len){for (int i = 0; i < len;i++){des[bit_rev(i, len-1)] = src[i];}}下⾯定义计算2维IFFT的函数:void GPUFFT::ifft_2d(ofFbo& in,ofFbo& out,int size);其中in是输⼊,out是输出,size就是N,由初始化的时候传⼊了⼀次,在这⾥写是为了⽅便调试的时候临时改变尺⼨。
实验二 基-2FFT 算法的软件实现一、实验目的1、 加深对DFT 算法原理和基本性质的理解;2、 熟悉FFT 算法的流程;3、 了解FFT 算法的应用。
二、基本原理1、 DFT 算法原理 (见教材第三章)2、按时间抽取(DIT )的-2FFT 算法(1)算法原理序列x (n )的N (N =2-M )点DFT 为kn N N n N Wn x n x DFT k X ∑-===1)()]([)(点,k =0, 1, …, N -1 (2.1)将式(2.1)按n 的奇偶性分解为)12(12/212/)12()2()()()(+-=-===∑∑∑∑++=+=l k NN n l k NN n kn Nn kn Nn Wl x Wl x Wn x Wn x k X 奇数偶数奇数偶数(2.2)令)2()(1l x l x =, )12()(2+=l x l x ,因为klN lk NW W 2/2=, 所以式(2.2)可写成)12(2/122/12/012)()()(+--=-=∑∑+=l k N M n k NklN N l Wl x WWl x k X 奇数(2.3)式(2.3)说明,按n 的奇偶性将x (n )分解为两个N /2长的序列x 1(l )和x 2(l ),则N 点DFT 可分解为两个N /2点DFT 来计算。
用X 1(k )和X 2(k )分别表示12,...,1,0 )()]([)(12/02/12/11-===∑-=Nk Wl x l x DFT k X N l kl N N 点 (2.4) 12,...,1,0 )()]([)(12/02/22/22-===∑-=Nk W l x l x DFT k X N l kl N N 点 (2.5)将(2.4)式和(2.5)式代入(2.31)式,并利用kN Nk N W W -=+2和X 1(k )、 X 2(k )的隐含周期性可得到:12,...1,0)()()2()()()(2121-=⎪⎭⎪⎬⎫-=++=N k k X W k X N k X k X W k X k X kN kN (2.6) 这样,将N 点DFT 的计算分解为计算两个N/2点离散傅立叶变换X 1(k )和X 2(k ),再计算(2.6)式。
为了将如上分解过程用运算流图表示,以便估计其运算量,观察运算规律,总结编程方法,先介绍一种表示(2.6)式的蝶形图。
图2.1蝶形运算图图2.2 8点DFT 一次时域抽取分解运算流图根据图2.2可以求得第一次分解后的运算量。
图2.2包括两个N/2点DFT 和N/2个蝶形,每个N/2点DFT 需要2)2/(N 次复数乘法和2/)12/(N N -次复数加法运算,每个蝶形只有一次复数乘法运算和两次复数加法运算。
所以,总的复数乘法次数为总的复数加法次数为2|)1(222)2(212N N N N N N ≈+=+⨯>> (2.7)22222)12(2N N N N =⨯+⨯⨯- (2.8) N=8点DIT-FFT 的运算流图如图2.3(a)所示。
根据W k N/m =W km N ,将图2.3(a)转换成如图2.3(b)所示的标准形式的运算流图图2.3 N=8点DIT-FFT 的运算流图(2)算法效率由图2.3可见,N =2M 时,DIT-FFT 运算流图由M 级蝶形构成,每级有N/2个蝶形。
因此,每级需要N/2次复数乘法运算和N 次复数加法运算,M 级形共需复数乘法次数C M (2)和复数加法次数C A (2)分别为lbN NM N C M 22)2(=⋅=(2.9) C A (2) =N ·M =N lb N (2.10)式中,lb N=log 2 N 。
直接计算N 点DFT 的复数乘法次数为N 2,复数加法次数为(N -1)N 。
当N 1时, N 2/C M(2)>>1,所以N 越大,DIT-FFT 运算效率越高。
DIT-FFT 算法与DFT 所需乘法次数与N 的关系曲线如图2.4所示。
例如,N =210=1024时,DIT-FFT 的运算效率为8.20410210241024)2(22=⨯==-M C N FFT DIT DFT 的乘法次数的乘法次数 (2.11)而当N=211=2048时,37.372112048222)2(22≈⨯==⨯=M N M N N C N M (2.12)图2.4 DIT-FFT 与DFT 所需乘法次数比较曲线(3)运算规律r N W 的确定第m 级运算,一个DIT 蝶型运算的的两接点“距离”为12-m ,所以r Nm m m m m rNm m m m Wk X k X k X W k X k X k X )2()()2()2()()(1111111-------+-=+++= (2.13)r 的求解:(1)将式(2.13)中第一个节点的标号k 表示成L (LN 2=)位二进制数;(2)把此二进制数乘上m L -2,即将L 位二进制数左移m L -位(注意m 是第m 级运算),把右边空出的位置补0,此数即为所求r 的二进制数。
序列的倒序DIT-FFT 算法的输入序列的排序看起来似乎很乱,但仔细分析就会发现这种倒序是很有规律的。
由于LN 2=,因此顺序数可用L 位二进制数(0121n n n n L L --)表示。
M 次偶奇时域抽选过程如图2.5所示。
第一次按最低位0n 的0和1将x(n)分解为偶奇两组,第二次又按次低位1n 的0、 1值分别对偶奇组分解; 依次类推,第L 次按1-L n 位分解,最后所得二进制倒序数如图2.5所示。
图2.5 形成例序的树状图(N =23)形成倒序J 后,将原存储器中存放的输入序列重新按倒序排列。
设原输入序列x (n )先按自然顺序存入数组A 中。
例如,对N =8, A (0),A (1),A (2),…,A (7)中依次存放着x (0),x (1), …, x (7)。
对x(n)的重新排序(倒序)规律如图2.6所示。
图2.6倒序规律三、实验仪器计算机四、实验要求及内容用所学过的编程语言,自行设计出一个按时间抽取的、输入倒位序、输出顺序的基-2 FFT 算法程序。
五、实验报告(1)简述实验目的及原理;(2)画出程序流程框图;(3)主要给出实验内容的程序(要求程序模块化并加注释)。
程序流程框图实验的程序#include<iostream.h> #include<math.h>#include<string.h>#define PI 3.1415926 class Plural//复数类{float real;//实部float img;//虚部public:Plural(float r,float i){real=r;img=i;}Plural(float r)//通过构造函数重载使用数与复数乘法{real=r;img=0;}Plural(){real=0;img=0;}friend Plural* fft(float X[],int n); //fft();friend Plural operator *(Plural p1,Plural p2);//重载乘法运算符friend Plural operator -(Plural p1,Plural p2);//重载减法运算符friend Plural operator +(Plural p1,Plural p2);//重载加法运算符friend Plural* daoxu(Plural* x[],int n); //倒序Plural operator =(Plural p) ; //重载赋值符号friend Plural W(int N,int i) ; //生成旋转因子void show() ; //输出real+imgi};Plural operator *(Plural p1,Plural p2)//复数乘法{return Plural(p1.real *p2.real -p1.img *p2.img ,p1.real *p2.img +p1.img *p2.real ); }Plural operator +(Plural p1,Plural p2)//复数加法{return Plural(p1.real+p2.real,p1.img +p2.img );}Plural operator -(Plural p1,Plural p2)//复数加法{return Plural(p1.real-p2.real,p1.img -p2.img );}Plural Plural::operator =(Plural p) //重载赋值符号{real=p.real ;img=p.img ;return *this;}//********************生成旋转因子***************************Plural W(int N,int i){float r;float im;r=(float)cos(2*PI*float(i)/float(N));im=(float)sin((-1)*2*PI*float(i)/float(N));return Plural(r,im);}//*********************输出函数show()************************* void Plural::show() //输出real+imgi{if(real==0){if(img==0)cout<<0<<" ";elsecout<<img<<"i ";}else{if(img<0)cout<<real<<img<<"i ";elseif(img==0)cout<<real<<" ";elsecout<<real<<"+"<<img<<"i ";}}/************2的n次方*******************/int mi2(int n){int m=1;for(int i=0;i<n;i++){m*=2;//m=2^n}return m;}/*************由二进制数转化为十进制数******************/int twoten(int xu[],int i){int m=0;for(int j=0;j<i;j++){m+=xu[j]*mi2(i-j-1);//m=xu[0]*1+xu[1]*2+xu[2]*4+xu[3]*8+xu[4]*16+...}return m;}/*************对二进制数倒序加一*****************/void add(int xu[],int i){xu[0]++;for(int j=0;j<i;j++){if(xu[j]==2){xu[j]=0;xu[j+1]++;//例如xu[]=0000,1000,0100,1100,0010......}}}/*************倒序**********************************/void daoxu(Plural x[],int n){float m=float (n);int i=0;for(i=0;m>1;i++)//得到n是2的多少次方{m/=2;//m=log(n)/log(2)}int M=mi2(i);int* xu=new int[i];//定义一个长度为i的数组,当做2进制数for(int j=0;j<i;j++){xu[j]=0;}Plural* jia=new Plural[M];//定义一个长度为M的数组jia,以保存倒序for(j=0;j<M;j++){int mm;mm=twoten(xu,i);//得到数组xu[]所对应的十进制数jia[j]=x[mm];add(xu,i);//二进制数组最左端加1}for(j=0;j<n;j++)//撤销动态生成的数组{x[j]=jia[j];}if(jia)delete[] jia;}//***************快速傅里叶变换FFT()************************* Plural* fft(Plural X[],int m){int n=1;for(;m>mi2(n);n++);//得到m等于n次方Plural* A=new Plural[m];for(int i=0;i<n;i++)//共有n级{for(int j=0;j<mi2(n-1-i);j++)//第i级有2^(n-i)族{for(int k=0;k<mi2(i);k++)//每族有2^(m-1)个蝶形运算{int m,n;m=k+j*mi2(i+1);n=k+mi2(i)+j*mi2(i+1);Plural a1,a2;if(i==0){a1=X[m]+W(mi2(i+1),k)*X[n]; //蝶形运算a2=X[m]-W(mi2(i+1),k)*X[n];//蝶形运算A[m]=a1;A[n]=a2;}else{a1=A[m]+W(mi2(i+1),k)*A[n]; //蝶形运算a2=A[m]-W(mi2(i+1),k)*A[n];//蝶形运算A[m]=a1;A[n]=a2;}}}}for(i=0;i<m;i++)//撤销动态生成的数组{X[i]=A[i];}delete[] A;return X;}/****************************主函数*******************/ void main(){int N=0;cout<<"输入变换的点数N:"<<endl;cin>>N;Plural *x=new Plural[N] ;cout<<"数据实部:虚部(少于N点时以#结束):"<<endl;for(int i=0;i<N;i++){char s;float x1,x2;cin>>x1;s=(char)x1;if(s=='#'){ break;}else{cin>>x2;x[i]=Plural(x1,x2);}x1=0;x2=0;}cout<< "您输入的数组:"<<endl;for( i=0;i<N;i++){x[i].show();if((i+1)%4==0)cout<<endl;}cout<<endl;daoxu(x,N);//调用倒序daoxu()函数cout<< "倒序后的数组:"<<endl;for( i=0;i<N;i++){x[i].show();if((i+1)%4==0)cout<<endl;}cout<<endl;x=fft(x,N);//调用fft()函数cout<< "FFT变换后数组"<<endl;for( i=0;i<N;i++){x[i].show();if((i+1)%4==0)cout<<endl;}if(x)delete[] x;}通过本次实验:1、加深了对DFT算法原理和基本性质的理解;2、熟悉了FFT算法的流程;3、了解FFT算法的应用。