光谱学与光谱分析
Spectroscopy and Spectral Analysis
基金项目：国家自然科学基金：基于多特征关联的复杂地形高分辨率遥感图像匹配技术研究(41201349) 作者简介：杨雪，1986年生，主要从事：遥感算法和空间数据处理
基于GPU 的巨幅影像分块技术快速傅里叶变换算法研究 杨雪1,2，李学友3,4，李家国2*，马骏1，余涛2，杨健2
1.河南大学计算机信息工程学院，河南开封 475004
2.中国科学院遥感应用研究所, 北京 100101;
3.北京四维空间数码科技有限公司, 北京 100039;
4.中国测绘科学研究院, 北京 100830;
摘 要 快速傅里叶变换(FFT)是遥感影像处理的基础方法，随着高光谱、高空间和高时间分辨率遥感影像获取能力的提升，如何利用快速傅里叶变换技术处理巨幅遥感影像是当前遥感影像处理技术中的重要环节和研究热点。本文在分析传统快速傅里叶变换算法的基础上，提出了一种基于GPU 加速的巨幅遥感影像快速傅里叶变换（Huge Remote Fast Fourier Transform, HRFFT ）算法。通过对CUFFT 函数库中的FFT 算法进行改进，解决了巨幅图像内存或显存溢出的问题。并结合HJ-1A 卫星CCD 遥感影像，进行对比分析研究，证明了该方法的合理性。该算法在遥感图像处理项目进行实际应用，取得了很好的效果。 关键词 GPU ，HRFFT ，快速傅里叶变换，遥感影像 中图分类号: TP751 文献标志码: A
引 言
傅里叶变换是数字图像处理和技术分析的基础，是实现图像滤波、噪声去除、图像压缩和MTF 补偿等过程的基本环节[1]。Cooley 等1965年提出了快速傅里叶变换（FFT ）算法[2]，加快了离散傅里叶变换（DFT ）的计算速度，减少了数据的计算量，使整个图像处理技术的研究有了突破性的进展。此后，人们不断研究一些新的算法和软件来提高海量数据的快速傅里叶变换实现算法。1999年MIT 的MatteoFrigo 和Steven G.Johnson 共同开发的FFTW （Fastest Fourier Transform in the West ）包，是目前世界公认的执行速度最快的FFT 软件，支持一维和多维实数及复数的DFT 并包含对共享和分布式存储系统的并行变换[3]，成功应用于ENVI 和MATLAB 等商业软件。近年来，随着可编程图像处理器（GPU ）的出现，其强大
的并行运算能力、高精度的浮点运算能力以及高速的存储器带宽，将海量数据的快速傅里叶变换推向新的台阶。研究表明，在GPU 环境下运用统一设备架构（CUDA ）的CUFFT 进行FFT 处理的速度大约是CPU 上的20倍左右[4]。
虽然GPU 处理的速度有了质的飞跃，但与CPU 受到物理内存的限制一样，GPU 在进行FFT 运算时也受到显存容量的限制，特别是随着科学技术的不断发展，遥感影像的光谱、空间和时间分辨率不断提高，单景影像的存储量已超过2GB ，合成后的影像数据甚至高达几十GB 。如何在有限的条件下实现GB 级遥感影像的FFT 快速运算，是当前研究的热点问题之一。本文从FFT 基本原理着手，通过将数据分块处理和GPU 加速相结合，突破FFT 运算过程中对巨幅影像大小的限制，并有效地降低了运算的时间。
1 HRFFT 算法
1.1 FFT
算法基础
()x n （01n N ≤≤-）的离散傅里叶变换（DFT ）
可以表示为[5]：
1
0()[()]()N nk N
n X k DFT x n x n W -===∑，
01k N ≤≤- （1）
其中：
2j k n kn N
N
W
e
π-= （2）
将一维傅里叶变换扩展到二维为：
11
2(
)00
(,)(,)M N ux vy
j M N
x y F u v f x y e
π---+===∑∑，
（01x M ≤≤-，01y N ≤≤-） （3） 根据二维傅里叶变换的可分离性，可以将一个二维傅里叶变换转化为两个一维傅里叶变换。快速傅里叶变换（FFT
）算法是运用因子的周期性和对称性，逐步分解为短序列的DFT ，再对其组合成原序列DFT ，进而减少运算量。FFT 算法可以分为时间抽取算法（DFT)和频率抽取算法（DIT ）。不同的N 值选择的FFT 算法基数也不同。对于定点FFT 有许多不同的算法结构[6,7]：基-2、基-4及任意按时间抽取，基-2、基-4及任意按频率抽取等。基-2按频率抽取算法模型蝶形图如图1所示。
N x (0)x (4)x (2)x (6)x (1)x (5)x (3)x (7)
X (0)X (1)X (2)X (3)
X (4)X (5)X (6)X (7)
图1 8点DIT―FFT 运算流图
Fig.1 8 DIT - FFT operation flow chart
从图1可以看出，通过两个点FFT 蝶形算法可以求出下一个点的值，不断递归此过程，直到经过第logN 步并行计算后，才可以算出一个N 点的一维FFT 。在计算过程中，并行节点之间需要进行数据传递，而且每一次循环都需要同步执行，才可以进行下一步运算。由于频域抽取算法的输出结果与正常顺序存在码位倒置的关系，因此需要将其转换后才可以得到最终的输出结果。
1.2 HRFFT 算法
受计算机内存容量的限制，巨幅遥感影像需要经过分块处理才能进行FFT 运算。例如文献[8]提出了一种对影像FFT 处理的中间结果分块存储的方法。该方法虽然能够处理巨幅影像，但由于全部处理过程由CPU 完成，执行效率比较低。实验表明，对于一副大约1.6G 的40000×40000大小的影像，处理时间甚至超过24小时，远远不能满足实际应用的需要。
GPU 的出现为图像处理提供了新的实现方法。Moreland 和Angel 尝试利用GPU 实现快速傅里叶变换[9]，由于受当时GPU 软硬件的限制，虽然最终的结果还不如那时的FFTW 库的效率高，但该文献却展示了运用GPU 在实现FFT 算法上的一个很好的思路和应用前景。2007年NVIDIA 公司开发了利用GPU 实现影像FFT 运算的CUFFT 库，文献[10]对基于GPU 的CUFFT 库与FFTW 算法进行了实验对比，文献[11-16]均运用GPU 对FFT 算法进行改进，实验表明，这些改进的FFT 算法对一般容量的图像而言均具有较高的实现效率，但对于巨幅影像的处理，仍存在一些问题。
本文在上述文献的基础上，提出了基于GPU 加速的巨幅遥感影像快速傅里叶变换的HRFFT 算法。基本处理思路是先将巨幅影像按行分块读入内存，然后利用GPU 对每一行块的多个子块进行一维傅里叶变换和渲染；所有行的块处理完毕后，再对每一列块的图像进行一维傅里叶变换；最终实现二维傅里叶变换。该算法的关键是内存块的大小和显存子块大小的确定，以及在GPU 中处理和渲染子块时如何将中间结果暂存到磁盘的。 1.3
分块处理方法
巨幅影像数据每次读入内存中的块的大小和在GPU 中进行处理的子块的大小，需要从影像数据存储格式、每像素占用字节数、内存可用空间大小以及GPU 可用显存容量等多个方面综合进行考虑[17]。具体的分块方法如图2所示。
图2 巨幅影像FFT 变换分块处理示意图 Fig.2 The huge image FFT transform block diagram
1.3.1 CPU 处理分块大小的确定
FFT 分块算法分为行优先分块和列优先分块两种形式。本文采用固定列数，行优先分块的原则对FFT 进行分块处理，即先对行进行变换，再对行变换后的结果进列变换。
假设每个像素点在内存中所占的字节数为p ，则一个m ×n 的影像分块数可以用公式表示为：
1m n p b a ??=+???? （4） 其中，a 为处理图像时系统的可用内存，该值可
以通过函数自动获取。此时读入内存的每一块影像的字节数为：
m n p
s b =
（5）
1.3.2 GPU 处理子块大小的确定
由于显存容量较小，因此还需要将读入到内存中的块划分为更小的子块，才能利用GPU 进行高效处理。子块的大小会影响到数据在GPU 中运算的速度和数据读写时与I/O 的交互速度。分的子块太小，将中间结果暂存到磁盘上时，文件读写的次数就很多，增加了I/O 之间的频繁交互；分的子块太大，就会导致显存溢出。
子块大小的计算与读入到内存中的块大小计算相似，不同的是还需要考虑显存的纹理数据存储结构以及在GPU 中进行渲染时所占用的显存空间。 1.3.3 分块处理流程
第1步，首先，根据遥感影像和GPU 显存的大小，判断是否对遥感影像进行分块。如果遥感影像能够被GPU 处理，则将遥感影像全部读入GPU 显存中，利用GPU 并行计算进行遥感影像的二维傅立叶变换；如果遥感影像过大，不能被GPU 处理，则转入下步。
第2步，根据显存大小将整个图像分成若干块。先将第一块读入内存中。
第3步，将第一块划分成多个可以在GPU 中处理的子块。
第4步，对第一子块数据进行一维傅里叶行变换，并将变换的中间结果渲染到纹理，接下来对第二子块进行一维傅里叶行变换，直到变换的结果数据与纹理大小相等时，然后将结果进行分块，写入磁盘中，记下地址。
第5步，重复第4步，将第1块的所有子块进行一维傅里叶行变换，并对每一次变换的结果进行分块，写入磁盘，追加到之前的结果块的地址后面，存入到磁盘中。
第6步，重复以上步骤，直到整幅图像的一维傅里叶行变换结束。
第7步，根据影像结果存储的位置，按照之前所分块的结果数据，读取每一列块的影像数据。
第8步，对一列块的影像数据按列进行一维傅里叶列变换，记下每次变换结果的地址，写入磁盘。
第9步，重复步骤7～8，完成所有列块的处理完，得到傅里叶正变换后的结果。
傅里叶逆变换和正变换过程相反。 1.4
渲染加速与中间结果保存
显存中的数据都是以纹理方式进行存储的。由于CPU 中数据是以数据形式进行存储的，实现中间结果的存储很容易，统一的存储器模型决定它可以在程序中任何地方进行存储器的读写操作，GPU 采用的是流式编程模型，运算通常是线性的[18]，运算时需要先把数据存放到纹理中，由于数据运算存在依赖性，实现迭代操作就很困难。因此，为了减少读写操作的次数，本文采用渲染到纹理技术对CUFFT 库进行改进。
首先构造一个高速显存的缓存区Pbuffer ，将每一次CUFFT 算法的处理结果直接输出到Pbuffer ，并绑定到纹理，作为CUFFT 逆变换的输入数据，当计算的影像数据比较大时，可以同时构造两个Pbuffer ，交替存放输出的计算结果和输入的结果。这样就避免了将中间结果不断写入磁盘，减少了与I/O 交互的次数。大大提高了运算的速度。
2 试验与分析
本文以2009年4月12日的HJ-1A 卫星CCD1影像2级产品（序列号为96068）第一波段影像数据为例，影像大小（行×列）16000×14000，左上角经度是117.33，右下角的经度120.62，左上角的纬度是：25.48，右下角的纬度是21.67。
实验环境是：CPU ： Intel （R ）Core （TM ）i7-2600 Duo 3.4GHZ ，硬盘转速7200RPM ，4.00GB 内存，Windows 7，64位操作系统；GPU ：NVIDIA
GEFORCE GT 520，容量1024MB 的显存。测试程序用Visual Studio 2008开发环境编写。通过对比实验验证HRFFT 算法的加速效率。 2.1 不分块的对比试验
验证在不分块情况下FFTW 和HRFFT 算法对图像进行FFT 变换的效率。采用对原始影像进行重采样得到的影像大小（行×列）为64×64，128×128，256×256，512×512，1024×1024，2048×2048，4096×4096，8192×8192，16000×14000，16000×
16000的数据，分别进行试验，并记录其处理时间如表1所示。
表1 FFTW 算法与GPU 算法处理时间对比（单位：s ） Table 1 FFTW algorithm and GPU algorithm processing
time contrast (unit: s)
从表1可以看出，在计算较小的图像时，运用FFTW 的时间小于GPU 算法处理的时间，但是随着影像的增大，运用GPU 就大大缩短了FFT 处理的时间，当处理图像（行×列）小于8192×7000时，FFTW 和GPU 都可以处理，当图像（行×列）大于8192×7000时，影像的大小就超过了可用内存的大小，用FFTW 处理时主要依靠在内存中进行数据传输和计算，而运用GPU 处理时，进行结果显示时采用异步传输的方式将数据从内存中直接拷贝到事先开辟好高速缓冲区的显存中，通过驱动程序构造纹理，减少大量数据在内存的操作。当图像（行×列）大于16000×10000时，用GPU 计算时，数据的处理超过了可用显存的限制，因此，出现显存溢出的问题。
为了实现更大遥感影像的FFT 变换，就需要采用分块存储中间结果的方式进行处理。 2.2
与分块结果比对
本文仍以HJ-1A 卫星CCD1第一波段影像为例，对原始影像进行重采样，得到（行×列）16000×12000、20000×20000、30000×20000、40000×30000、40000×40000、80000×80000的数据，分别对其进行处理，由于实际影像的大小小于表中图像的大小，还需要对图像进行重采样得到被测试的数据。
实验中采用三种方法进行巨幅影像的处理比较。
第一种方法：FFTW 分块算法。主要思想是将图像处理的中间结果以固定缓存来存储，依据缓存的最大容量，每次处理多行数据，并将处理结果组成一块，先在缓存中进行转置，并写入到磁盘中。 第二种方法：CUFFT 分块算法。该算法主要思想是将中间结果分块存储、在读写数据时进行缓存优化并采用多线程技术进行并行处理。
第三种方法：本文提出的HRFFT 算法。该方法在传统方法的基础上，对CUFFT 算法进行了改进，由于CUFFT 算法库本身提供的FFT 运算并没有对算法进行改进，只是用传统的实数FFT 的接口形式对CUFFT 的复数FFT 进行了封装，算法效率也没有提高，因此，本文采用运用了GPU 的并行处理能力，并采用内存分块和GPU 显存分块的形式将处理的中间结果直接渲染到纹理来减少磁盘读写的次数。并可以进行巨幅影像的处理。
实验结果如表2和表3所示。
表2 三种算法正变换运算时间（单位：s ） Table 2 Three algorithms are transform operation time
(unit: s)
图3 遥感影像FFT 算法正变换运算时间对比图 Fig 3. Remote sensing image FFT algorithm transform
computation time comparison chart
表3 三种算法逆变换运算时间（单位：s ） Table 3 The three algorithm inversion computation
time (unit: s)
图4 遥感影像FFT算法逆变换运算时间对比图
Fig.4 Remote sensing image FFT algorithm inversion
computation time comparison chart
从上述表中和图中可以看出，进行FFT处理时，对于（行×列）16000×12000的影像三种方法的时间对比还不是很明显，但是随着影像数据的增大，运用本文算法所用的时间明显比其他两种算法所用的时间短，并且保证了三种方法处理结果的一致性。充分发挥了GPU并行处理的性能。特别是处理巨幅影像时，无论是正变换还是逆变换，本文提出的进一步改进算法在时间上都有很大提高。影像越大，表现的优势就越明显。
3 结论
通过对不同影像快速傅里叶变换算法对比分析，采用FFTW分块算法与CUFFT算法，在处理巨幅影像时都受到限制，而采用本文提出的HRFFT 算法一方面解决了巨幅遥感影像处理时遇到的内存溢出和显存溢出的问题，另一方面采用渲染到纹理的方法来存放CUFFT运算结果，使得遥感影像进行FFT处理的速度明显高于前两种算法，并且影像越大，效果越明显。而FFTW算法又优越于一些图像处理软件如ERDAS、PCI等。因此，采用HRFFT 算法，为以后遥感图像的处理带来了很大便利。解决了以往巨幅影像不能处理或者处理速度慢的问题。应用在实际的遥感图像处理项目中，节省了处理时间，提高了处理效率。Reference
[1]LING Y R,EHLERS M,USERY E L,et al.FFT-enhanced HIS
transform method for fusing high-resolution satellite images[J].ISPRS Journal of Photogrammetry & Remote Sensing,2007(61):381-392.
[2]J. W. Cooley and J. W. Tukey. An Algorithm for the Machine
Calculation of Complex Fourier Series. Math. Comput., Vol.
19:297–301, 1965.
[3]K. Singh, J.P. Walters, J. Hestness, J. Suh,C.M. Rogers,S.P.
Crago, FFTW and Complex Ambiguity Function performance on the Maestro processor. Aerospace Conference, 2011 IEEE.5-12 March 2011, pages 1-8.
[4]D. Pairman. S.E. Belliss. J. Cuff. S.J. McNeill. Detection and
mapping of irrigated farmland in Canterbury, New Zealand.
Geoscience and Remote Sensing Symposium (IGARSS),2011, IEEE International. 24-29 July 2011.pages 696-699.
[5]李仕,王晶,孙辉.基于图形处理器的实数FFT在图形处理中的
应用[J].光学精密工程.2008,12,16(12).2414-2420.
[6]Yue Wu, Weile Jia, Lin-Wang Wang, Weiguo Gao, Long Wang,
Xuebin Chi. GPU Tuning for First-Principle Electronic Structure Simulations.Springer Berlin Heidelberg. 2013. pages 235-246.
[7]A. Nukada. Y. Ogata, T. Endo and S. Matsuoka. Bandwidth
intensive 3-D FFT kernel for GPUs using CUDA. In Proceedings of the 2008 ACM/IEEE conference on Supercomputing, IEEE Press, 2008. pp 1-11.
[8]宋江洪,赵忠明,王刚.提高海量图像数据变换域算法速度的方
法研究[J].计算机辅助设计与图形学学报.
2005.17(7).1517-1521.
[9]Kenneth Moreland and Edward Angel. The FFT on GPU [A].In
Proceedings of Graphics Hardware, San Diego,2003.112-119.
[10]肖江,胡柯良,邓元勇.基于CUDA的矩阵乘法和FFT性能测试
[J].计算机工程.2009.5.35(10):7-10.
[11]Steven M. de Jong and Freek D. van der Meer. Remote
Sensing Image Analysis[M]. P. O. Box 17,3300AA Dordrecht.
The Netherlands. 2006:26-29.
[12]N. K. Govindaraju. B. Lloyd. Y. Dotsenko. B. Smith and J.
Manferdelli. High performance discrete Fourier transforms on graphics processors. In Proceedings of the 2008 ACM/IEEE conference on Supercomputing, pages 1-12. IEEE Press, 2008.
[13]Y. Chen, X. Cui and H. Mei. Large-scale FFT on GPU
clusters. 24th International Conference on Supercomputing (ICS'10), pp 315-324, ACM 2010.
[14]Long Wang, Weile Jia, Xuebin Chi, Yue Wu, Weiguo Gao,
Lin-Wang Wang, Large Scale Plane Wave Pseudopotential Density Functional Theory Calculations on GPU Clusters，Supercomputing Conference 2011. （SCA2011会议论文）[15]杨靖宇.遥感影像GPU并行化处理技术与实现方法[D].郑州：
解放军信息工程大学硕士学位论文.2008.
[16]Y. Ogata, T. Endo, N. Maruyama, and S. Matsuoka. An efficient,
model-based CPU-GPU heterogeneous FFT library. In Parallel and Distributed Processing, 2008. IPDPS 2008. IEEE
International Symposium on. April 2008. Pp 1–10.
[17]Gac, N.S. Mancini, M. Desvignes and D. Houzet. "High speed
3D tomography on CPU, GPU and FPGA", EURASIP Journal on Embedded Systems-Special issue on design and architectures for signal and image processing.Vol.2008, Article
5 (2008).
[18]汤颖,张宏鑫,张美玉.基于图形硬件的纹理图像编码与实时绘
制算法[J].计算机学报,2007.30(2):272-280.
Research on Fast Fourier Transforms Algorithm of Huge Remote Sensing Image Segmentation Technology with GPU
Xue Yang 1,2，Xueyou Li 3,4，Jaguo Li2*，Jun Ma1，Tao Yu2
1. School of Computer and Information Engineering, Henan University, KaiFeng, HeNan, China, 475004;
2. Institute of Remote Sensing Applications Chinese Academy Of Sciences, Beijing, China, 100101;
3.Beijing Siwei Spatial Data Technology CO.,Ltd, Beijing, China, 100039;
4.Chinese Academy of Surveying &Mapping, BeiJing,China, 100830 ;
Abstract:Fast Fourier Transforms (FFT) is base method of image processing. As the productive of the obtaining ability of remote sensing image that has the features of hyperspectrum、high spatial resolution、high time resolution, utilizing FFT technology to handle huge remote sensing image is the critical part and research hot spot of image processing field today. On the basis of analyzing FFT algorithms, this paper put forward a algorithm of remote sensing image FFT based on GPU acceleration HRFFT algorithm by the improvement of the FFT algorithm in the CUFFT function library of CUDA, this algorithm solved the memory and video memory overflow problem faced by huge image. Combined with the CCD image of HJ-1A satellite, this paper proved the feasibility of this method by comparing it with other algorithms through experiments.
Key words: GPU, HRFFT, FFT, Remote Sensing Image
信号与系统的基本思想：把复杂的信号用简单的信号表示，再进行研究。 怎么样来分解信号？任何信号可以用Delta 函数的移位加权和表示。只有系统是线性时不变系统，才可以用单位冲激函数处理，主要讨论各个单位冲激函数移位加权的响应的叠加能得到总的响应。 线性系统（齐次性，叠加定理） 时不变系统 对一个系统输入单位冲激函数，得到的响应为h(t).表征线性时不变系统的非常重要的东西，只要知道了系统对单位冲击函数的响应，就知道了它对任何信号的响应，因为任何信号都可以表示为单位冲激函数的移位加权和。 例如：d(t)__h(t) 那么a*d(t-t0)__a*h(t-t0) -()= ()(t-)d f t f τδττ∝∝? 的响应为-y()=()(-)t f h t d τττ∝ ∝ ? 记为y(t)=f(t)*h(t),称为f(t)和h(t)的卷积 总结为两点：对于现行时不变系统，任何信号可以用单位冲激信号的移位加权和表示，任何信号的响应可以用输入函数和单位冲激函数响应的卷积来表示 连续时间信号和系统的频域分析 时域分析的重点是把信号分解为单位冲激函数的移位加权和，只讨论系统对单位冲激函数的响应。而频域的分析是把信号分解为各种不同频率的正弦函数的加权和，只讨论系统对sinwt 的响应。都是把信号分解为大量单一信号的组合。
周期函数可以展开为傅里叶级数，将矩形脉冲展开成傅里叶级数，得到傅里叶级数的系数 n A sin F = T x x τ 其中0=2 nw x τ。 取样函数sin ()=x S a x 。产生一种震荡，0点的值最大，然后渐渐衰减直至0 第一：对于傅里叶级数的系数，n 是离散的，所以频谱也是离散状的每条谱线都出现在基波频率的整数倍上，其包络是取样函数。 第二：谱线的间距是0w .。零点是0=2nw x τ，02w =T π是谱的基波频率。如果τ不变，T 增大，那么0w 减小，当T 非常大的时候，0w 非常小，谱线近似连续，越来越密，幅度越来越小。 傅里叶变换：非周期函数 正变换：--F jw)= ()iwt f t e dt ∝ ∝?（ 反变换：-1()=()2jnwt f t F jw e dw π ∝∝ ? 常用函数的傅里叶变换（典型非周期信号的频谱）
实验八 利用FFT 实现快速卷积 一、 实验目的 （1） 通过这一实验，加深理解FFT 在实现数字滤波（或快速卷积）中的重要作用，更好的利用FFT 进行数字信号处理。 （2） 进一步掌握循环卷积和线性卷积两者之间的关系。 二、 实验原理与方法 数字滤波器根据系统的单位脉冲响应h(n)是有限长还是无限长可分为有限长单位脉冲响应（Finite Impulse Response ）系统（简记为FIR 系统）和无限长单位脉冲响应（Infinite Impulse Response ）系统（简记为IIR 系统）。 对于FIR 滤波器来说，除了可以通过数字网络来实现外，也可以通过FFT 的变换来实现。 一个信号序列x(n)通过FIR 滤波器时，其输出应该是x(n)与h(n)的卷积： ∑+∞ -∞ =-= =m m n h m x n h n x n y )()()(*)()( 或 ∑+∞ -∞ =-= =m m n x m h n x n h n y ) ()()(*)()( 当h(n)是一个有限长序列，即h(n)是FIR 滤波器，且10-≤≤N n 时 ∑-=-=1 0) ()()(N m m n x m h n y 在数字网络（见图6.1）类的FIR 滤波器中，普遍使用的横截型结构（见下图6.2 图6.1 滤波器的数字网络实现方法 图6.2 FIR 滤波器横截型结构 y(n) y(n) -1-1-1-1
应用FFT 实现数字滤波器实际上就是用FFT 来快速计算有限长度列间的线性卷积。 粗略地说，这种方法就是先将输入信号x(n)通过FFT 变换为它的频谱采样 值X(k)，然后再和FIR 滤波器的频响采样值H(k)相乘，H(k)可事先存放在存储器中，最后再将乘积H(k)X(k)通过快速傅里叶变换（简称IFFT ）还原为时域序列，即得到输出y(n)如图6.3所示。 图6.3 数字滤波器的快速傅里叶变换实现方法 现以FFT 求有限长序列间的卷积及求有限长度列与较长序列间的卷积为例来讨论FFT 的快速卷积方法。 （1） 序列)(n x 和)(n h 的列长差不多。设)(n x 的列长为1N ，)(n h 的列长为2N ，要求 )()(n x n y =N ∑-=-==1 ) ()()(*)()(N r r n h r x n h n x n h 用FFT 完成这一卷积的具体步骤如下： i. 为使两有限长序列的线性卷积可用其循环卷积代替而不发生混叠，必须选择循环卷积长度121-+≥N N N ，若采用基2-FFT 完成卷积运 算，要求m N 2=（m 为整数）。 ii. 用补零方法使)(n x ，)(n h 变成列长为N 的序列。 ?? ?-≤≤-≤≤=10 10)()(11N n N N n n x n x ?? ?-≤≤-≤≤=10 1 0)()(22N n N N n n h n h iii. 用FFT 计算)(),(n h n x 的N 点离散傅里叶变换 )()(k X n x FFT ??→? )()(k H n h FFT ??→? iv. 做)(k X 和)(k H 乘积，)()()(k H k X k Y ?= v. 用FFT 计算)(k Y 的离散傅里叶反变换得 y(n)
第五章 离散傅里叶变换及其快速算法 1 离散傅里叶变换(DFT)的推导 (1) 时域抽样： 目的：解决信号的离散化问题。 效果：连续信号离散化使得信号的频谱被周期延拓。 (2) 时域截断： 原因：工程上无法处理时间无限信号。 方法：通过窗函数(一般用矩形窗)对信号进行逐段截取。 结果：时域乘以矩形脉冲信号，频域相当于和抽样函数卷积。 (3) 时域周期延拓： 目的：要使频率离散，就要使时域变成周期信号。 ! 方法：周期延拓中的搬移通过与)(s nT t -δ的卷积来实现。 表示：延拓后的波形在数学上可表示为原始波形与冲激串序列的卷积。 结果：周期延拓后的周期函数具有离散谱。 (4) 1。 图1 DFT 推导过程示意图 (5) 处理后信号的连续时间傅里叶变换：∑∑ ∞ -∞=-=π--δ???? ? ????= k N n N kn j s kf f e nT h f H )()()(~ 010/2
(i) )(~ f H 是离散函数，仅在离散频率点S NT k T k kf f == =00处存在冲激，强度为k a ，其余各点为0。 (ii) )(~ f H 是周期函数，周期为s s T NT N T N Nf 100===，每个周期内有N 个不同的幅值。 (iii) 时域的离散时间间隔(或周期)与频域的周期(或离散间隔)互为倒数。 2 DFT 及IDFT 的定义 (1) ， (2) DFT 定义：设()s nT h 是连续函数)(t h 的N 个抽样值1,,1,0-=N n ，这N 个点的宽度为N 的DFT 为：[])1,...,1,0(,)()(1 /2-=???? ??==? -=π-∑N k NT k H e nT h nT h DFT s N n N nk j s s N (3) IDFT 定义：设??? ? ??s NT k H 是连续频率函数)(f H 的N 个抽样值1,,1,0-=N k ， 这N 个点 的宽度为N 的IDFT 为： ())1,...,1,0(,11 0/21 -==??? ? ? ?=??? ???? ?? ??? ???-=π--∑ N k nT h e NT k H N NT k H DFT s N k N nk j s s N (4) N nk j e /2π-称为N 点DFT 的变换核函数，N nk j e /2π称为N 点IDFT 的变换核函数。它们互 为共轭。 (5) 同样的信号，宽度不同的DFT 会有不同的结果。DFT 正逆变换的对应关系是唯一的， 或者说它们是互逆的。 (6) 引入N j N e W /2π-= (i) 用途： (a) 正逆变换的核函数分别可以表示为nk N W 和nk N W -。 (b) 核函数的正交性可以表示为：() )(* 1 0r n N W W kr N N k kn N -δ=∑-= (c) DFT 可以表示为：)1,,1,0(,)(10 -==? ??? ??∑ -=N k W nT h NT k H N n nk N s s (d) IDFT 可以表示为：)1,,1,0(,1 )(10 -=??? ? ??=∑-=-N n W NT k H N nT h N k nk N s s (ii) ） (iii) 性质：周期性和对称性： (a) 12==π-j N N e W (b) 12 /-==π-j N N e W (c) r N r N N N r N N W W W W ==+ (d) r N r N N N r N N W W W W -=-=+2/2/ (e) )(1Z m W m N ∈?= (f) ),(/2/2Z n m W e e W n N N n j mN mn j mn mN ∈?===π-π- 3 离散谱的性质 (1) 离散谱定义：称)(Z k NT k H H S k ∈??? ? ??=? 为离散序列)0)((N n nTs h <≤的DFT 离散谱，简称离散谱。 (2) 性质：
从头到尾彻底理解傅里叶变换算法、上 从头到尾彻底理解傅里叶变换算法、上 前言 第一部分、DFT 第一章、傅立叶变换的由来 第二章、实数形式离散傅立叶变换（Real DFT） 从头到尾彻底理解傅里叶变换算法、下 第三章、复数 第四章、复数形式离散傅立叶变换 前言： “关于傅立叶变换，无论是书本还是在网上可以很容易找到关于傅立叶变换的描述，但是大都是些故弄玄虚的文章，太过抽象，尽是一些让人看了就望而生畏的公式的罗列，让人很难能够从感性上得到理解”---dznlong， 那么，到底什么是傅里叶变换算法列?傅里叶变换所涉及到的公式具体有多复杂列? 傅里叶变换（Fourier transform）是一种线性的积分变换。因其基本思想首先由法国学者傅里叶系统地提出，所以以其名字来命名以示纪念。 哦，傅里叶变换原来就是一种变换而已，只是这种变换是从时间转换为频率的变化。这下，你就知道了，傅里叶就是一种变换，一种什么变换列?就是一种从时间到频率的变化或其相互转化。 ok，咱们再来总体了解下傅里叶变换，让各位对其有个总体大概的印象，也顺便看看傅里叶变换所涉及到的公式，究竟有多复杂： 以下就是傅里叶变换的4种变体（摘自，维基百科） 连续傅里叶变换 一般情况下，若“傅里叶变换”一词不加任何限定语，则指的是“连续傅里叶变换”。连续傅里叶变换将平方可积的函数f（t）表示成复指数函数的积分或级数形式。
这是将频率域的函数F(ω)表示为时间域的函数f（t）的积分形式。 连续傅里叶变换的逆变换(inverse Fourier transform)为： 即将时间域的函数f（t）表示为频率域的函数F(ω)的积分。 一般可称函数f（t）为原函数，而称函数F(ω)为傅里叶变换的像函数，原函数和像函数构成一个傅里叶变换对（transform pair）。 除此之外，还有其它型式的变换对，以下两种型式亦常被使用。在通信或是信号处理方面，常以来代换，而形成新的变换对: 或者是因系数重分配而得到新的变换对： 一种对连续傅里叶变换的推广称为分数傅里叶变换（Fractional Fourier Transform）。分数傅里叶变换(fractional Fourier transform,FRFT)指的就是傅里叶变换(Fourier transform,FT)的广义化。 分数傅里叶变换的物理意义即做傅里叶变换a 次，其中a 不一定要为整数；而做了分数傅里叶变换之后，信号或输入函数便会出现在介于时域(time domain)与频域(frequency domain)之间的分数域(fractional domain)。
C语言实现FFT(快速傅里叶变换) 函数原型:空快速傅立叶变换(Struct Compx *xin，Intn) 函数函数:对输入复数组执行快速傅立叶变换(FFT)输入参数:*xin复结构组的第一个地址指针。结构输出参数:no * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *结构compx u，w，t。 nv2 =快速傅立叶变换_ N/2；nm1 =快速傅立叶变换_ N-1；(I = 0；i
MATLAB实验傅里叶分析
实验七 傅里叶变换 一、实验目的 傅里叶变换是通信系统、图像处理、数字信号处理以及物理学等领域内的一种重要的数学分析工具。通过傅里叶变换技术可以将时域上的波形分 布变换为频域上的分布，从而获得信号的频谱特性。MATLAB 提供了专门的函数fft 、ifft 、fft2(即2维快速傅里叶变换)、ifft2以及fftshift 用于实现对信号的傅里叶变换。本次实验的目的就是练习使用fft 、ifft 以及fftshift 函数，对一些简单的信号处理问题能够获取其频谱特性（包括幅频和相频特性）。 二、实验预备知识 1. 离散傅里叶变换(DFT)以及快速傅里叶变换(FFT)简介 设x (t )是给定的时域上的一个波形，则其傅里叶变换为 2()() (1)j ft X f x t e dt π∞--∞=? 显然X ( f )代表频域上的一种分布(波形)，一般来说X ( f )是复数。而傅里叶逆变换定义为： 2()() (2)j ft x t X f e df π∞-∞ =?
因此傅里叶变换将时域上的波形变换为频域上的波形，反之，傅里叶逆变换则将频域上的波形变换为时域上的波形。 由于傅里叶变换的广泛应用，人们自然希望能够使用计算机实现傅里叶变换，这就需要对傅里叶变换(即(1)式)做离散化处理，使 之符合电脑计算的特征。另外，当 把傅里叶变换应用于实验数据的分 析和处理时，由于处理的对象具有 离散性，因此也需要对傅里叶变换 进行离散化处理。而要想将傅里叶 变换离散化，首先要对时域上的波 形x (t )进行离散化处理。采用一个 时域上的采样脉冲序列: δ (t －nT ), n = 0, 1, 2, …, N －1； 可以实现上述目的，如图所示。其中N 为采样点数，T 为采样周期；f s = 1/T 是采样频率。注意采样时，采样频率f s 必须大于两倍的信号频率（实际是截止频率），才能避免混迭效应。 接下来对离散后的时域波形()()()(x t x t t n T x n T δ= -=的傅里叶变换()X f 进行离散处理。与上述做法类 似，采用频域上的δ脉冲序列： x (t δ x (t )δ t t t
第一章 信号与系统的基本概念 1．信号、信息与消息的差别？ 信号：随时间变化的物理量； 消息：待传送的一种以收发双方事先约定的方式组成的符号，如语言、文字、图像、数据等 信息：所接收到的未知内容的消息，即传输的信号是带有信息的。 2．什么是奇异信号？ 函数本身有不连续点或其导数或积分有不连续点的这类函数统称为奇异信号或奇异函数。例如： 单边指数信号 （在t =0点时，不连续）， 单边正弦信号 （在t =0时的一阶导函数不连续）。 较为重要的两种奇异信号是单位冲激信号δ(t )和单位阶跃信号u(t )。 3．单位冲激信号的物理意义及其取样性质？ 冲激信号：它是一种奇异函数，可以由一些常规函数的广义极限而得到。 它表达的是一类幅度很强，但作用时间很短的物理现象。其重要特性是筛选性，即： ()()()(0)(0)t x t dt t x dt x δδ∞ ∞ -∞ -∞ ==? ? 4．什么是单位阶跃信号？ 单位阶跃信号也是一类奇异信号，定义为： 10()00t u t t >?=?
12()()()x t ax t bx t =+，其中a 和b 是任意常数时， 输出信号()y t 是1()y t 和2()y t 的线性叠加，即：12()()()y t ay t by t =+； 且当输入信号()x t 出现延时，即输入信号是0()x t t -时， 输出信号也产生同样的延时，即输出信号是0()y t t -。 其中，如果当12()()()x t x t x t =+时，12()()()y t y t y t =+，则称系统具有叠加性； 如果当1()()x t ax t =时，1()()y t ay t =则称系统具有均匀性。 线性时不变系统是最基本的一类系统，是研究复杂系统，如非线性、时变系统的基础。 6．线性时不变系统的意义与应用？ 线性时不变系统是我们本课程分析和研究的主要对象，对线性时不变性进行推广，可以得到线性时不变系统具有微分与积分性质，假设系统的输入与输出信号分别为()x t 和()y t ，则 当输入信号为 ()dx t dt 时，输出信号则为() dy t dt ； 或者当输入信号为()t x d ττ-∞ ?时，输出信号则为()t y d ττ-∞ ?。 另外，线性时不变系统对信号的处理作用可以用冲激响应（或单位脉冲响应）、系统函数或频率响应进行描述。而且多个系统可以以不同的方式进行连接，基本的连接方式为：级联和并联。 假设两个线性时不变系统的冲激响应分别为：1()h t 和2()h t ， 当两个系统级联后，整个系统的冲激响应为：12()()*()h t h t h t =； 当两个系统并联后，整个系统的冲激响应为：12()()()h t h t h t =+； 当0t <时，若()0h t =， 则此系统为因果系统； 若|()|h t dt ∞ -∞<∞?， 则此系统为稳定系统。 第二章 连续时间系统的时域分析 1．如何获得系统的数学模型？ 数学模型是实际系统分析的一种重要手段，广泛应用于各种类型系统的分析和控制之中。 不同的系统，其数学模型可能具有不同的形式和特点。对于线性时不变系统，其数学模型
快速傅里叶变换FFT的FPGA设计与实现 学生姓名郭衡 班级电科1704 学号17419002064 指导教师谭会生 成绩 2020年5 月20 日
快速傅里叶变换FFT 的设计与实现 一、研究项目概述 非周期性连续时间信号x(t)的傅里叶变换可以表示为：= )(?X dt t j e t x ? ∞ ∞ --1 )(?，式中计算出来的是信号x(t)的连续频谱。但是，在实际的控制系统中能够式中计算出来的是信号x(t)的连续频谱。但是，在实际的控制系统中能够算信号x(t)的频谱。 有限长离散信号x(n)，n=0，1，…，N-1的DFT 定义为： ∑-=-=-==1 02,1.....10)()(N n N j N kn N e W N k W n x K X π、、。 可以看出，DFT 需要计算大约N2次乘法和N2次加法。当N 较大时，这个计算量是很大的。利用WN 的对称性和周期性，将N 点DFT 分解为两个N ／2点的DFT ，这样两个N ／2点DFT 总的计算量只是原来的一半，即(N ／2)2+(N ／2)2=N2／2，这样可以继续分解下去，将N ／2再分解为N ／4点DFT 等。对于N=2m 点的DFT 都可以分解为2点的DFT ，这样其计算量可以减少为(N ／2)log2N 次乘法和Nlog2N 次加法。图1为FFT 与DFT-所需运算量与计算点数的关系曲线。由图可以明显看出FFT 算法的优越性。 图1 FFT 与DFT 所需乘法次数比 较
X[1] 将x(n)分解为偶数与奇数的两个序列之和，即x(n)=x1(n)+x2(n)。 x1(n)和x2(n)的长度都是N ／2，x1(n)是偶数序列，x2(n)是奇数序列，则 ∑∑=--=-=+2 )12(120 2)1.....,0()(2)(1)(N n k n N N n km N N k W n x W n x K X 所以)1...,0()(2)(1)(12 22120 -=+=∑∑-=-=N k W n x W W n x K X N n km N k N km N N n 由于km N N j km N j km N W e e W 2/2 /2222===--ππ ，则 )1.....,0)((2)(1)(2)(1)(12 2/120 2/-=+=+=∑∑-=-=N k k X W k X W n x W W n x K X k N N n km N k N N n kn N 其中X1(k)和X2(k)分别为x1(n)和x2(n)的N ／2点DFT 。由于X1(k)和X2(k)均以N ／2为周期，且WNk+N/2=-WNk ，所以X(k)又可表示为： )12/....,1,0)((2)(1)(-=+=N k k X W k X K X k N )12/....,1,0)((2)(1)2/(-=-=+N k k X W k X N K X k N
function x=MyIFFT_FB(y) %MyIFFT_TB:My Inverse Fast Fourier Transform Time Based %按频率抽取基2-傅里叶逆变换算法 %input: % y -- 傅里叶正变换结果,1*N的向量 %output: % x -- 逆变换结果，1*N的向量 %参考文献： % https://www.doczj.com/doc/624147954.html,/view/fea1e985b9d528ea81c779ee.html N=length(y); x=conj(y); %求共轭 x=MyFFT_FB(x);%求FFT x=conj(x);%求共轭 x=x./N;%除以N end %% 内嵌函数====================================================== function y=MyFFT_FB(x,n) %MYFFT_TB:My Fast Fourier Transform Frequency Based %按频率抽取基2-fft算法 %input: % x -- 输入的一维样本 % n -- 变换长度，缺省时n=length(x) 当n小于x数据长度时，x数据被截断到第n个数据% 当n大于时，x数据在尾部补0直到x 含n个数据 %output: % y -- 1*n的向量，快速傅里叶变换结果 %variable define: % N -- 一维数据x的长度 % xtem -- 临时储存x数据用 % m,M -- 对N进行分解N=2^m*M,M为不能被2整除的整数 % two_m -- 2^m % adr -- 变址，1*N的向量 % l -- 当前蝶形运算的级数 % W -- 长为N/2的向量，记录W(0,N),W(1,N),...W(N/2-1,N) % d -- 蝶形运算两点间距离 % t -- 第l级蝶形运算含有的奇偶数组的个数 % mul -- 标量，乘数 % ind1,ind2 -- 标量，下标 % tem -- 标量，用于临时储存 %参考文献： % https://www.doczj.com/doc/624147954.html,/view/fea1e985b9d528ea81c779ee.html %% 输入参数个数检查
时域信号 弧频率表示的 傅里叶变换 注释 1 线性 2 时域平移 3 频域平移, 变换2的频域对应4 如果值较大，则会收缩 到原点附近，而会扩 散并变得扁平. 当 | a | 趋向 无穷时，成为 Delta函数。 5 傅里叶变换的二元性性质。通过 交换时域变量和频域变量 得到. 6 傅里叶变换的微分性质 7 变换6的频域对应 8 表示和的卷积—这
9 矩形脉冲和归一化的sinc 函数 10 变换10的频域对应。矩形函数是理想的低通滤波器，sinc 函数是这类滤波器对反因果冲击的响应。 11 tri 是三角形函数 12 变换12的频域对应 13 高斯函数 exp( ? αt 2) 的傅里叶变换是他本身. 只有当 Re(α) > 0时，这是可积的。 14 15 16 a>0 17 变换本身就是一个公式
18 δ(ω) 代表狄拉克δ函数分布. 这 个变换展示了狄拉克δ函数的重要 性：该函数是常函数的傅立叶变换 19 变换23的频域对应 20 由变换3和24得到. 21 由变换1和25得到，应用了欧拉公 式: cos(at) = (e iat + e?iat) / 2. 22 由变换1和25得到 23 这里, n是一个自然数. δ(n)(ω) 是狄拉克δ函数分布的n阶微分。这 个变换是根据变换7和24得到的。 将此变换与1结合使用，我们可以变 换所有多项式。 24 此处sgn(ω)为符号函数；注意此变 换与变换7和24是一致的. 25 变换29的推广. 26 变换29的频域对应. 27 此处u(t)是单位阶跃函数; 此变换 根据变换1和31得到.
哈尔滨商业大学 DSP课程设计报告 题目快速傅立叶变换（FFT）算法专业电子信息工程 班级 08级02班 姓名学号王玉辉200810930172 李砚秋200810930062 杨兴臻200810930292 尤琳200810930052 指导教师姜海涛 日期2011年12月12日
目录 1.设计目的 .... 错误！未定义书签。 1.1. 设计目的..................................................................... 错误！未定义书签。 1.2. 使用设备..................................................................... 错误！未定义书签。 2.设计任务与要求错误！未定义书签。 3.原理与分析 .. 错误！未定义书签。 4.实验步骤 .... 错误！未定义书签。 5.软件设计 .... 错误！未定义书签。 6.系统仿真及调试错误！未定义书签。 7.完成结果或效果错误！未定义书签。 8.心得体会 .... 错误！未定义书签。 9.参考文献 .... 错误！未定义书签。
1. 设计目的 1.1. 设计目的 1．掌握用窗函数法设计FFT 快速傅里叶的原理和方法； 2．熟悉FFT 快速傅里叶特性； 3．了解各种窗函数对快速傅里叶特性的影响。 1.2. 使用设备 PC 兼容机一台，操作系统为Windows2000(或Windows98，WindowsXP ，以下默认为Windows2000)，安装Code Composer Studio 2.0 软件。 2. 设计任务与要求 按原程序仿真完成后，修改参数，观察波形变化。 3. 原理与分析 1． FFT 的原理和参数生成公式 )()()()()(212 12 22 12 1k X W k X W r x W W r x k x k N rk N N r k N rk N N r +=+=∑∑-=-= 公式（1）FFT 运算公式 FFT 并不是一种新的变换，它是离散傅立叶变换（DFT ）的一种快速算法。由于我们在计算DFT 时一次复数乘法需用四次实数乘法和二次实数加法；一次复数加法则需二次实数加法。每运算一个X （k ）需要4N 次复数乘法及2N+2（N-1）=2（2N-1）次实数加法。所以整个DFT 运算总共需要4N^2 次实数乘法和N*2(2N-1)=2N(2N-1)次实数加法。如此一来，计算时乘法次数和加法次数都是和N^2 成正比的，当N 很大时，运算量是可观的，因而需要改进对DFT 的算法减少运算速度。 根据傅立叶变换的对称性和周期性，我们可以将DFT 运算中有些项合并。 我们先设序列长度为N=2^L ，L 为整数。将N=2^L 的序列x(n)(n=0,1,……，N-1)，按N 的奇偶分成两组，也就是说我们将一个N 点的DFT 分解成两个N/2 点的DFT ，他们又重新组合成一个如下式所表达的N 点DFT ： 一般来说，输入被假定为连续的。当输入为纯粹的实数的时候，我们就可以利用左右对称的特性更好的计算DFT 。 我们称这样的RFFT 优化算法是包装算法：首先2N 点实数的连续输入称为“进包”。其次N 点的FFT 被连续被运行。最后作为结果产生的N 点的合成输出是“打开”成为最
快速傅里叶变换 FFT是离散傅立叶变换的快速算法，可以将一个信号变换到频域。有些信号在时域上是很难看出什么特征的，但是如果变换到频域之后，就很容易看出特征了。这就是很多信号分析采用FFT变换的原因。另外，FFT可以将一个信号的频谱提取出来，这在频谱分析方面也是经常用的。 虽然很多人都知道FFT是什么，可以用来做什么，怎么去做，但是却不知道FFT之后的结果是什意思、如何决定要使用多少点来做FFT。 现在圈圈就根据实际经验来说说FFT结果的具体物理意义。一个模拟信号，经过ADC采样之后，就变成了数字信号。采样定理告诉我们，采样频率要大于信号频率的两倍，这些我就不在此罗嗦了。 采样得到的数字信号，就可以做FFT变换了。N个采样点，经过FFT之后，就可以得到N个点的FFT结果。为了方便进行FFT运算，通常N取2的整数次方。 假设采样频率为Fs，信号频率F，采样点数为N。那么FFT之后结果就是一个为N点的复数。每一个点就对应着一个频率点。这个点的模值，就是该频率值下的幅度特性。具体跟原始信号的幅度有什么关系呢？假设原始信号的峰值为A，那么FFT的结果的每个点（除了第一个点直流分量之外）的模值就是A的N/2倍。而第一个点就是直流分量，它的模值就是直流分量的N倍。而每个点的相位呢，就是在该频率下的信号的相位。第一个表示直流分量（即0Hz），而最后一个点N的再下一个点（实际上这个点是不存在的，这里是假设的第N+1个点，也可以看做是将第一个点分做两半分，另一半移到最后）则表示 采样频率Fs，这中间被N-1个点平均分成N等份，每个点的频率依次增加。例如某点n所表示的频率为：Fn=(n-1)*Fs/N。由上面的公式可以看出，Fn所能分辨到频率为为Fs/N，如果采样频率Fs为1024Hz，采样点数为1024点，则可以分辨到1Hz。1024Hz的采样率采样1024点，刚好是1秒，也就是说，采样1秒时间的信号并做FFT，则结果可以分析到1Hz，如果采样2秒时间的信号并做FFT，则结果可以分析到0.5Hz。如果要提高
第七章 傅里叶变换 1．求下列函数的傅氏变换： (1)1,10, ()1, 01,0,; t f t t --<? 解: (1)[()]()j t F f t f t e dt ω+∞--∞ =? 1 101 10 1 1 22sin cos | 2(1cos ).j t j t j t j t e dt e dt e dt e dt j i tdt t j ωωωωωωω ωω -----=-+=-+=-= =- -????? (2) ()()j t F f t e dt ωω+∞--∞ =? 0(1)(1)0 11|.11t j t j t j t e e dt e dt e j j ωωωωω ---∞ -∞ --∞====--?? 6．求下列函数的傅氏变换 (1) 1,0,sgn 1,0;t t t -? (2) ()sin(5).3f t t π =+ 解: (1)已知 1 [()](),[1]2(),F u t F j πδωπδωω = +=由sgn 2()1t u t =-有 12[sgn ]2( ())2().F t j j πδωπδωωω =+-= (2) 由于 1()sin(5)sin 5cos5,322f t t t t π=+=+ 故 [()][(5)(5)](5)(5)].2j F f t πδωδωδωδω= +--++- 7．已知00()[()()]F ωπδωωδωω=++-为函数()f t 的傅氏变换，求().f t
#include
#include ————第四章———— 快速傅里叶变换FFT 所谓的快速算法，就是根据原始变换定义算法的运算规律及其中的某些算子的特殊性质，找出减少乘法和加法运算次数的有效途径，实现原始变换的各种高效算法。一种好的快速算法可使变换速度提高几个数量级。 由于快速算法很多，而且还在不断研究和发展。较成熟的算法都有现成的程序。所以，通过教材中介绍的四种快速算法，主要学习研究快速算法的基本思想和减少运算量的途径，熟悉各种快速算法的特点、运算效率和适用情况。为今后研究新的快速算法和合理选用快速算法打好基础。 4.1 学 习 要 点 4.1.1 直接计算N 点DFT 的运算量 对于 ()(),1 0∑-==N n kn N W n x k X 1,,1,0-=N k 复数乘法次数： 2 N M c = 复数加法次数： ()1-=N N A c 当1>>N 时，复数乘法和加法次数都接近为2 N 次，随着N 增大非线性增大。 4.1.2 减少运算量的基本途径 DFT 定义式中只有两种运算：()n x 与kn N W 的乘法相加。所以，kn N W 的特性对乘法运算 量必有影响。 （1）根据的对称性、周期性和特殊值减少乘法运算次数。 ①对称性：k N N k N W W -=+ 2 ，()k k N N W 12-=，()k N k N N W W =* - ②周期性：k N lN k N W W =+。 ③kn N W 的特殊值(无关紧要旋转因子)： 1;;124 -===±N N N N N W j W W 。对这些因子不能进行乘法运算。 （2）将较大数N 点DFT 分解为若干个小数点DFT 的组合，减少运算量。这正是FFT 能大量节省运算量的关键。 4.1.3 四种快速算法的基本思想及特点 根据上述减少运算量的途径，巧妙地在时域或频域进行不同的抽取分解与组合，得到不 傅里叶算法的采样电流计算 ******* 广西大学******* 摘要：微机继电保护是用数学运算的方法实现故障的测量、分析和判断的。通过全波傅立叶算法可用于求出各次谐波分量的幅值和相角，并具有一定的滤波作用。本文探讨了傅氏算法在电力系统中的应用。介绍了全波傅立叶算法的基本原理。通过仿真验证了该算法的实用性。 关键词：微机继电保护；电力系统；算法 引言 在微机保护装置中，首先要对反映被保护设备的电气量模拟量进行采集，然后对这些采集的数据进行数字滤波，再对这些经过数字滤波的数字信号进行数学运算、逻辑运算，并进行分析判断，最终输出跳闸命令、信号命令或计算结果，以实现各种继电保护功能。这种对数据进行处理、分析、判断以实现保护功能的方法称为算法。目前广泛采用全波傅氏算法和最小二乘法作为电力系统微机保护提取基波分量的算法。 傅立叶算法可用于求出各谐波分量的幅值和相角，所以它在微机保护中作为计算信号幅值的算法被广泛采用。实际上，傅立叶算法也是一种滤波方法。分析可知，全周傅氏算法可有效滤除恒定直流分量和各正次谐波分量。 傅里叶算法原理 一个周期函数满足狄里赫利条件，就可以将这个周期函数分解为一个级数，最为常用的级数是傅里叶级数，傅氏算法的基本思路来自傅里叶级数， 即一个周期性函数可以分解为直流分量、基波分量及各次谐波的无穷级数，如 ∑∞ =+=011)()]sin()cos([n n n t t nw a t nw b i （1.1） 式中1w 表示基波角频率；n a 和n b 分别是各次谐波的正弦和余弦的幅值， 其中比较特殊的有：0b 表示直流分量，11,b a 表示基波分量正、余弦项的幅 值。根据傅氏级数的原理，可以求出n a 、n b 分别为 ?=T t n dt t nw i T a 0 1)()sin(2 （1.2） ?=T t n dt t nw i T b 0 1)()cos(2 （1.3） 于是n 次谐波电流分量可表示为 )sin()cos()(11t nw a t nw b t i n n n += （1.4） 据此可求出n 次谐波电流分量的有效值和相角为 ???????=+=n n n n n n a b a b a I arctan 222 （1.5） 其中n a 、n b 可用梯形积分法近似求出为快速傅里叶变换FFT.
傅里叶算法的采样电流计算
相关主题
文本预览