电子科技大学通信与信息工程学院标准实验报告
（实验）课程名称DSP设计与实践
电子科技大学教务处制表
电 子 科 技 大 学
实 验 报 告
学生姓名： 学 号 指导教师：实验地点： 实验时间： 一、实验室名称： 科B341 二、实验项目名称：快速傅立叶变换（FFT ）的实现 三、实验学时：4 四、实验原理：
基—2按时间抽取FFT 算法
对于有限长离散数字信号{x[n]}，0 ≤ n ≤ N-1，其离散谱{x[k]}可以由离散付氏变换（DFT ）求得。DFT 的定义为
可以方便的把它改写为如下形式：
不难看出，W N 是周期性的，且周期为N ，即
W N 的周期性是DFT 的关键性质之一。为了强调起见，常用表达式W N 取代W 以便明确其周期是N 。
由DFT 的定义可以看出，在x[n]为复数序列的情况下，完全直接运算N 点DFT 需要（N-1）2次复数乘法和N （N-1）次加法。因此，对于一些相当大的N 值（如1024）来说，直接计算它的DFT 所作的计算量是很大的。FFT 的基本思想在于，将原有的N 点序列序列分成两个较短的序列，这些序列的DFT 可以很简单的组合起来得到原序列的DFT 。例如，若N 为偶数，将原有的N 点序列分成两个（N/2）点序列，那么计算N 点DFT 将只需要约[(N/2)2 ·2]=N 2/2次复数乘法。即比直接计算少作一半乘法。因子（N/2）2
表示直接计算（N/2）点DFT 所需要的乘法次数，而乘数2代表必须完成两个DFT 。上述处理方法可以反复使用，即（N/2）点的DFT 计算也可以化成两个（N/4）点的DFT
()1
,...,1,0][)2(
1
-==--=∑N k e
n x k X nk N
j N n π()1
,...,1,0][1
0-==∑-=N k W n x k X nk
N
N n ...
2,1,0,))((±±==++l m W W nk
N
lN k mN n N
（假定N/2为偶数），从而又少作一半的乘法。这样一级一级的划分下去一直到最后就划分成两点的FFT 运算的情况。比如，一个N = 8点的FFT 运算按照这种方法来计算FFT 可以用下面的流程图来表示：
x(0)
x(1)
x(2)
x(3)x(4)
x(5)
x(6)
x(7)
X(7)
X(6)
X(5)
X(4)X(3)X(2)X(1)X(0)
关于蝶形结运算的具体原理及其推导可以参照讲义，在此就不再赘述。
实数FFT 运算
对于离散傅立叶变换（DFT ）的数字计算，FFT 是一种有效的方法。一般假定输入序列是复数。当实际输入是实数时，利用对称性质可以使计算DFT 非常有效。 一个优化的实数FFT 算法是一个组合以后的算法。原始的2N 个点的实输入序列组合成一个N 点的复序列，之后对复序列进行N 点的FFT 运算，最后再由N 点的复数输出拆散成2N 点的复数序列，这2N 点的复数序列与原始的2N 点的实数输入序列的DFT 输出一致。
使用这种方法，在组合输入和拆散输出的操作中，FFT 运算量减半。这样利用实数FFT 算法来计算实输入序列的DFT 的速度几乎是一般复FFT 算法的两倍。本实验就用这种方法实现了一个256点实数FFT （2N = 256）运算。
⒈ 实数FFT 运算序列的存储分配
如何利用有限的DSP 系统资源，合理的安排好算法使用的存储器是一个比较重要的问题。参见FFT 实验程序的CMD 文件：
MEMORY {
PAGE 0: IPROG: origin = 0x3080, len = 0x1F80 VECT: origin = 0x3000, len = 0x80 EPROG: origin = 0x38000, len = 0x8000
PAGE 1: USERREGS: origin = 0x60, len = 0x1c BIOSREGS: origin = 0x7c, len = 0x4 IDA TA: origin = 0x80, len = 0xB80 EDATA: origin = 0xC00, len = 0x1400
}
SECTIONS
{
.vectors: {} > VECT PAGE 0
.sysregs: {} > BIOSREGS PAGE 1
.trcinit: {} > IPROG PAGE 0
.gblinit: {} > IPROG PAGE 0
.bios: {} > IPROG PAGE 0
frt: {} > IPROG PAGE 0
.text: {} > IPROG PAGE 0
.cinit: {} > IPROG PAGE 0
.pinit: {} > IPROG PAGE 0
.sysinit: {} > IPROG PAGE 0
.data {} > EDA TA PAGE 1
.bss: {} > IDA TA PAGE 1
.far: {} > IDA TA PAGE 1
.const: {} > IDATA PAGE 1
.switch: {} > IDATA PAGE 1
.sysmem: {} > IDATA PAGE 1
.cio: {} > IDA TA PAGE 1
.MEM$obj: {} > IDA TA PAGE 1
.sysheap: {} > IDATA PAGE 1
}
从上面的连接定位CMD文件可以了解到，程序代码安排在0x3000开始的存储器中。其中0x3000-0x3080存放中断向量表。FFT程序使用的正弦表、余弦表数据（.data段）安排在0xc00开始的地方。变量(.bss段定义)存放在0x80开始的地址中。另外，本256点实数FFT程序的输入数据缓冲为0x2300-0x23ff，FFT后功率谱的计算结果存放在0x2200-0x22ff中。
⒉基二实数FFT运算的算法
该算法主要分为四步：
第一步，输入数据的组合和位倒序
把输入序列作位倒序，是为了在整个运算最后的输出中得到的序列是自然顺序。
首先，原始的输入的2N = 256个点的实数序列复制放到标记有“d_input_addr”的相邻单元，当成N = 128点的复数序列d[n]。奇数地址是d[n]的实部，偶数地址是d[n]的虚部。这个过程叫做组合（n是从0到无穷，指示时间的变量，N是常量）。然后，复数序列经过位倒序，存储在数据处理缓冲器中，标记为“fft-data”。
①如图2，输入实数序列为a[n]，n=0，1，2，3，…，255。分离a[n]成两个序
列，如图3所示。原始的输入序列是从地址0x2300到0x23FF，其余的从0x2200到0x22FF的是经过位倒序之后的组合序列：n=0，1，2，3， (127)
②d[n]表示复合FFT的输入，r[n]表示实部，i[n]表示虚部：
d[n]=r[n]+j i[n]
按位倒序的方式存储d[n]到数据处理缓冲中，如图2。
图2
* 编程技巧：在用C54x进行位倒序组合时，使用位倒序寻址方式可以大大提高程序执行的速度和使用存储器的效率。在这种寻址方式中，AR0存放的整数N是FFT点数的一半，一个辅助寄存器指向一数据存放的单元。当使用位倒序寻址把AR0加到辅助寄存器中时，地址以位倒序的方式产生，即进位是从左到右，而不是从右到左。例如，当AR0 = 0x0060，AR2 = 0x0040时，通过指令：
MAR AR2+0B
我们就可以得到AR2位倒序寻址后的值为0x0010。
下面是0x0060（1100000）与0x0040（1000000）以位倒序方式相加的过程：
1 1 0 0 0 0 0
1 0 0 0 0 0 0
+
0 0 1 0 0 0 0
实现256点数据位倒序存储的具体程序段如下：
bit_rev:
STM #d_input_addr,ORIGINAL_INPUT ；在AR3（ORIGINAL_INPUT）中
；放入输入地址
STM #fft_data,DATA_PROC_BUF ；在AR7（DATA_PROC_BUF）中
；放入处理后输出的地址
MVMM DA TA_PROC_BUF,REORDERED_DA TA ;AR2（REORDERED_DA TA）
;中装入第一个位倒序数据指针STM #K_FFT_SIZE-1,BRC
STM #K_FFT_SIZE,AR0 ;AR0的值是输入数据数目的一半=128
RPTB bit_rev_end
MVDD *ORIGINAL_INPUT+,*REORDERED_DATA+ ;将原始输入缓冲中的数据
;放入到位倒序缓冲中去之
;后输入缓冲(AR3)指针加1
;位倒序缓冲(AR2)指针也加
;一
MVDD *ORIGINAL_INPUT-,*REORDERED_DA TA+ ;将原始输入缓冲中的数据
;放入到位倒序缓冲中去之
;后输入缓冲(AR3)指针减一
;位倒序缓冲(AR2)指针加一
;以保证位倒序寻址正确MAR *ORIGINAL_INPUT+0B ；按位倒序寻址方式修改AR3 bit_rev_end:
注意，在上面的程序中。输入缓冲指针AR3（即ORIGINAL_INPUT）在操作时先加一再减一，是因为我们把输入数据相邻的两个字看成一个复数，在用寄存器间接寻址移动了一个复数（两个字的数据）之后，对AR3进行位倒序寻址之前要把AR3的值恢复到这个复数的首字的地址，这样才能保证位倒序寻址的正确。
第二步，N点复数FFT
在数据处理缓冲器里进行N点复数FFT运算。由于在FFT运算中要用到旋转因子W N，它是一个复数。我们把它分为正弦和余弦部分，用Q15格式将它们存储在两个分离的表中。每个表中有128项，对应从0度到180度。因为采用循环寻址来对表寻址，128 = 27 < 28，因此每张表排队的开始地址就必须是8个LSB位为0的地址。参照前面图1 DES 系统的存储区分配，所以我们必须把正弦表第一项“sine_table”放在0x0D00的位置,余弦表第一项
“cos_table ”放在0x0E00的位置。
① 根据公式
利用蝶形结对d[n]进行N=128点复数FFT 运算，其中
所需的正弦值和余弦值分别以Q15的格式存储于内存区以0x0D00开始的正弦表
和以0x0E00开始的余弦表中。我们把128点的复数FFT 分为七级来算，第一级是计算两点的FFT 蝶形结，第二级是计算四点的FFT 蝶形结，然后是八点、十六点、三十二点六十四点、一百二十八点的蝶形结计算。最后所得的结果表示为： D[k] = F{d[n]} = R[k] + j I[k]
其中，R[k]、I[k]分别是D[k]的实部和虚部。
图3
② FFT 完成以后，结果序列D[k]就存储到数据处理缓冲器的上半部分，如图3
()1
,...,1,0][1
0-==∑-=N k W n d k D nk
N
N n
)2sin()2cos(
)2(
nk N
j nk N e
W
nk N
j nk
N
πππ-==-
所示，下半部分任然保留原始的输入序列a[n]，这半部分将在第三步中被改写。
这样原始的a[n]序列的所有DFT的信息都在D[k]中了，第三步中需要做的就是就
是把D[k]变为最终的2N=256点复合序列，A[k]=F{a(n)}。
* 编程技巧：在实际的DSP的编程中为了节约程序运行时间，提高代码的
效率，往往要用到并行程序指令。比如：
ST B,*AR3
||LD *AR3+,B
并行指令的执行效果是，使原本分开要两个指令周期才能执行完的两条指令
在一个指令周期中中就能执行完。上述指令是将B移位（ASM-16）所决定
的位数，存于AR3所指定的存储单元中，同时并行执行，将AR3所指的单
元中的值装入到累加器B的高位中。由于指令的src和dst都是acc B，所以
存入*AR3中的值是这条指令执行以前的值。
这一步中，实现FFT计算的具体程序如下：
fft:
;---------stage1 ：计算FFT的第一步，两点的FFT
.asg AR1,GROUP_COUNTER ；定义FFT计算的组指针
.asg AR2,PX ；定义AR2为指向参加蝶形运算第一个数据的指针
.asg AR3,QX ；定义AR3为指向参加蝶形运算第二个数据的指针
.asg AR4,WR ；定义AR4为指向余弦表的指针
.asg AR5,WI ；定义AR5为指向正弦表的指针
.asg AR6,BUTTERFLY_COUNTER ；定义AR6为指向蝶形结的指针
.asg AR7,DATA_PROC_BUF ；定义在第一步中的数据处理缓冲指针
.asg AR7,STAGE_COUNTER ；定义剩下几步中的数据处理缓冲指针
pshm st0
pshm ar0
pshm bk ；保存环境变量
SSBX SXM ；开启符号扩展模式
STM #K_ZERO_BK,BK ；让BK=0 使*ARn+0% = *ARn+0
LD #-1,ASM ；为避免溢出在每一步输出时右移一位
MVMM DA TA_PROC_BUF,PX ；PX指向参加蝶形结运算的第一个数的实部（PR）
LD *PX,16,A ；AH := PR
STM #fft_data+K_DATA_IDX_1,QX ；QX指向参加蝶形结运算的第二个数的实部（QR）
STM #K_FFT_SIZE/2-1,BRC ；设置块循环计数器
RPTBD stage1end-1 ；语句重复执行的范围到地址stage1end-1处
STM #K_DA TA_IDX_1+1,AR0 ；延迟执行的两字节的指令（该指令不重复执行）
SUB *QX,16,A,B ；BH := PR-QR
ADD *QX,16,A ；AH := PR+QR
STH A,ASM,*PX+ ；PR':= (PR+QR)/2
ST B,*QX+ ；QR':= (PR-QR)/2
||LD *PX,A ；AH := PI
SUB *QX,16,A,B ；BH := PI-QI
ADD *QX,16,A ；AH := PI+QI
STH A,ASM,*PX+0% ；PI':= (PI+QI)/2
ST B,*QX+0% ；QI':= (PI-QI)/2
||LD *PX,A ；AH := next PR
stage1end:
; -------Stage 2 ：计算FFT的第二步，四点的FFT
MVMM DA TA_PROC_BUF,PX ；PX 指向参加蝶形结运算第一个数据的实部（PR）STM #fft_data+K_DATA_IDX_2,QX ；QX 指向参加蝶形结运算第二个数据的实部（QR）STM #K_FFT_SIZE/4-1,BRC ；设置块循环计数器
LD *PX,16,A ；AH := PR
RPTBD stage2end-1 ；语句重复执行的范围到地址stage1end-1处
STM #K_DA TA_IDX_2+1,AR0 ；初始化AR0以被循环寻址
;以下是第二步运算的第一个蝶形结运算过程
SUB *QX,16,A,B ；BH := PR-QR
ADD *QX,16,A ；AH := PR+QR
STH A,ASM,*PX+ ；PR':= (PR+QR)/2
ST B,*QX+ ；QR':= (PR-QR)/2
||LD *PX,A ；AH := PI
SUB *QX,16,A,B ；BH := PI-QI
ADD *QX,16,A ；AH := PI+QI
STH A,ASM,*PX+ ；PI':= (PI+QI)/2
STH B,ASM,*QX+ ；QI':= (PI-QI)/2
; 以下是第二步运算的第二个蝶形结运算过程
MAR *QX+ ；QX中的地址加一
ADD *PX,*QX,A ；AH := PR+QI
SUB *PX,*QX-,B ；BH := PR-QI
STH A,ASM,*PX+ ；PR':= (PR+QI)/2
SUB *PX,*QX,A ；AH := PI-QR
ST B,*QX ；QR':= (PR-QI)/2
||LD *QX+,B ；BH := QR
ST A, *PX ；PI':= (PI-QR)/2
||ADD *PX+0%,A ；AH := PI+QR
ST A,*QX+0% ；QI':= (PI+QR)/2
||LD *PX,A ；AH := PR
stage2end:
; Stage 3 thru Stage logN-1：从第三步到第六步的过程如下
STM #K_TWID_TBL_SIZE,BK ；BK = 旋转因子表格的大小值
ST #K_TWID_IDX_3,d_twid_idx ；初始化旋转表格索引值
STM #K_TWID_IDX_3,AR0 ；AR0 = 旋转表格初始索引值
STM #cos_table,WR ；初始化WR 指针
STM #sine_table,WI ；初始化WI 指针
STM #K_LOGN-2-1,STAGE_COUNTER ；初始化步骤指针
ST #K_FFT_SIZE/8-1,d_grps_cnt ；初始化组指针
STM #K_FLY_COUNT_3-1,BUTTERFLY_COUNTER ；初始化蝶形结指针
ST #K_DA TA_IDX_3,d_data_idx ；初始化输入数据的索引
stage: ；以下是每一步的运算过程
STM #fft_data,PX ；PX 指向参加蝶形结运算第一个数据的实部（PR）LD d_data_idx, A
ADD *(PX),A
STLM A,QX ；QX 指向参加蝶形结运算第二个数据的实部（QR）MVDK d_grps_cnt,GROUP_COUNTER ；AR1 是组个数计数器:
group: ；以下是每一组的运算过程
MVMD BUTTERFLY_COUNTER,BRC ；将每一组中的蝶形结的个数装入BRC
RPTBD butterflyend-1 ；重复执行至butterflyend-1处
LD *WR,T
MPY *QX+,A ；A := QR*WR || QX*QI
MAC *WI+0%,*QX-,A ；A := QR*WR+QI*WI
ADD *PX,16,A,B ；B := (QR*WR+QI*WI)+PR
; || QX指向QR
ST B,*PX ；PR':=((QR*WR+QI*WI)+PR)/2
||SUB *PX+,B ；B := PR-(QR*WR+QI*WI)
ST B,*QX ；QR':= (PR-(QR*WR+QI*WI))/2
||MPY *QX+,A ；A := QR*WI [T=WI]
; || QX指向QI
MAS *QX,*WR+0%,A ；A := QR*WI-QI*WR
ADD *PX,16,A,B ；B := (QR*WI-QI*WR)+PI
ST B,*QX+ ；QI':=((QR*WI-QI*WR)+PI)/2
; || QX指向QR
||SUB *PX,B ；B := PI-(QR*WI-QI*WR)
LD *WR,T ；T := WR
ST B,*PX+ ；PI':= (PI-(QR*WI-QI*WR))/2
; || PX指向PR
||MPY *QX+,A ；A := QR*WR || QX指向QI
butterflyend:
；更新指针以准备下一组蝶形结的运算
PSHM AR0 ；保存AR0
MVDK d_data_idx,AR0 ；AR0中装入在该步运算中每一组所用的蝶形结的数目MAR *PX+0 ；增加PX准备进行下一组的运算
MAR *QX+0 ；增加QX准备进行下一组的运算
BANZD group,*GROUP_COUNTER- ；当组计数器减一后不等于零时，延迟跳转至group处POPM AR0 ；恢复AR0（一字节）
MAR *QX- ；修改QX以适应下一组的运算
;
；更新指针和其他索引数据以变进入下一个步骤的运算
LD d_data_idx,A
SUB #1,A,B ；B = A-1
STLM B,BUTTERFLY_COUNTER ；修改蝶形结个数计数器
STL A,1,d_data_idx ；下一步计算的数据索引翻倍
LD d_grps_cnt,A
STL A,-1,d_grps_cnt ；下一步计算的组数目减少一半
LD d_twid_idx,A
STL A,-1,d_twid_idx ；下一步计算的旋转因子索引减少一半
BANZD stage,*STAGE_COUNTER-
MVDK d_twid_idx,AR0 ；AR0 = 旋转因子索引（两字节）
popm bk
popm ar0
popm st0 ；恢复环境变量
fft_end:
RET
第三步，分离复数FFT的输出为奇部分和偶部分
分离FFT输出为相关的四个序列：RP、RM、IP和IM，即偶实数，奇实数、偶虚数和奇虚数四部分，以便第四步形成最终结果。
①利用信号分析的理论我们把D[k]通过下面的公式分为偶实数RP[k]、奇实数
RM[k]、偶虚数IP[k]和奇虚数IM[k]：
RP[k] = RP[N-k] = 0.5 * (R[k] + R[N-k])
RM[k] = -RM[N-k] = 0.5 * (R[k] - R[N-k])
IP[k] = IP[N-k] = 0.5 * (I[k] + I[N-k])
IM[k] = -IM[N-k] = 0.5 * (I[k] - I[N-k])
RP[0] = R[0]
IP[0] = I[0]
RM[0] = IM[0] = RM[N/2] = IM[N/2] = 0
RP[N/2] = R[N/2]
IP[N/2] = I[N/2]
②下面的图4显示了第三步完成以后存储器中的数据情况，RP[k]和IP[k]存储在
上半部分，RM[k] 和IM[k]存储在下半部分。
图4
这一过程的程序代码如下所示：
unpack
.asg AR2,XP_k
.asg AR3,XP_Nminusk
.asg AR6,XM_k
.asg AR7,XM_Nminusk
STM #fft_data+2,XP_k ; AR2指向R[k] (temp RP[k])
STM #fft_data+2*K_FFT_SIZE-2,XP_Nminusk ; AR3指向R[N-K] (temp RP[N-K])
STM #fft_data+2*K_FFT_SIZE+3,XM_Nminusk ; AR7指向temp RM[N-K]
STM #fft_data+4*K_FFT_SIZE-1,XM_k ; AR6指向temp RM[K]
STM #-2+K_FFT_SIZE/2,BRC ; 设置块循环计数器
RPTBD phase3end-1 ; 从以下指令到phase3end-1处一直
; 重复执行BRC中规定的次数
STM #3,AR0 ; 设置AR0以备下面程序寻址使用
ADD *XP_k,*XP_Nminusk,A ; A := R[k]+R[N-K] = 2*RP[k]
SUB *XP_k,*XP_Nminusk,B ; B := R[k]-R[N-K]= 2*RM[k]
STH A,ASM,*XP_k+ ; 在AR[k]处存储RP[k]
STH A,ASM,*XP_Nminusk+ ; 在AR[N-K]处存储RP[N-K]=RP[k]
STH B,ASM,*XM_k- ; 在AI[2N-K]处存储RM[k]
NEG B ; B := R[N-K]-R[k] =2*RM[N-K]
STH B,ASM,*XM_Nminusk- ; 在AI[N+k]处存储RM[N-K]
ADD *XP_k,*XP_Nminusk,A ; A := I[k]+I[N-K] = 2*IP[k]
SUB *XP_k,*XP_Nminusk,B ; B := I[k]-I[N-K] =2*IM[k]
STH A,ASM,*XP_k+ ; 在AI[k]处存储IP[k]
STH A,ASM,*XP_Nminusk-0 ; 在AI[N-K]处存储IP[N-K]=IP[k]
STH B,ASM,*XM_k- ; 在AR[2N-K]处存储IM[k]
NEG B ; B := I[N-K]-I[k] =2*IM[N-K]
STH B,ASM,*XM_Nminusk+0 ; 在AR[N+k]处存储IM[N-K]
phase3end:
第四步，产生最后的N = 256点的复数FFT结果
产生2N = 256个点的复数输出，它与原始的256个点的实输入序列的DFT一致。输出驻留在数据缓冲器中。
①通过下面的公式由RP[k]、RM[n]、IP[n]和IM[n]四个序列可以计算出a[n]的
DFT：
AR[k] = AR[2N-k] = RP[k] + cos(kπ/N) * IP[k] - sin(kπ/N) * RM[k]
AI[k] = -AI[2N-k] = IM[k] - cos(kπ/N) * RM[k] - sin(kπ/N) * IP[k]
AR[0] = RP[0] + IP[0]
AI[0] = IM[0] - RM[0]
AR[N] = R[0] - I[0]
AI[N] = 0
其中：
A[k] = A[2N-k] = AR[k] + j AI[k] = F{a(n)}
②实数FFT输出按照实数/虚数的自然顺序填满整个4N=512个字节的数据处理
缓冲器。如图5所示。
图5
这一段程序可以和上面的步骤3的程序合成一个子程序unpack，这一步的程序代码如下：
ST #0,*XM_k- ; RM[N/2]=0
ST #0,*XM_k ; IM[N/2]=0
; 计算AR[0],AI[0], AR[N], AI[N]
.asg AR2,AX_k
.asg AR4,IP_0
.asg AR5,AX_N
STM #fft_data,AX_k ; AR2指向AR[0] (temp RP[0])
STM #fft_data+1,IP_0 ; AR4指向AI[0] (temp IP[0])
STM #fft_data+2*K_FFT_SIZE+1,AX_N ; AR5指向AI[N]
ADD *AX_k,*IP_0,A ; A := RP[0]+IP[0]
SUB *AX_k,*IP_0,B ; B := RP[0]-IP[0]
STH A,ASM,*AX_k+ ; AR[0] = (RP[0]+IP[0])/2
ST #0,*AX_k ; AI[0] = 0
MVDD *AX_k+,*AX_N- ; AI[N] = 0
STH B,ASM,*AX_N ; AR[N] = (RP[0]-IP[0])/2
; 计算最后的输出值AR[k], AI[k]
.asg AR3,AX_2Nminusk
.asg AR4,COS
.asg AR5,SIN
STM #fft_data+4*K_FFT_SIZE-1,AX_2Nminusk ; AR3指向AI[2N-1](temp RM[1])
STM #cos_table+K_TWID_TBL_SIZE/K_FFT_SIZE,COS ; AR4指向cos(k*π/N)
STM #sine_table+K_TWID_TBL_SIZE/K_FFT_SIZE,SIN ; AR5指向sin(k*π/N)
STM #K_FFT_SIZE-2,BRC
RPTBD phase4end-4
STM #K_TWID_TBL_SIZE/K_FFT_SIZE,AR0 ; AR0中存入旋转因子表的大
;小以备循环寻址时使用
LD *AX_k+,16,A ; A := RP[k] ||修改AR2指向IP[k]
MACR *COS,*AX_k,A ; A :=A+cos(k*π/N)*IP[k]
MASR *SIN,*AX_2Nminusk-,A ; A := A-sin(k*π/N)*RM[k] ||
;修改AR3指向IM[k]
LD *AX_2Nminusk+,16,B ; B := IM[k] ||修改AR3指向RM[k]
MASR *SIN+0%,*AX_k-,B ; B := B-sin(k*π/N)*IP[k] ||
;修改AR2指向RP[k]
MASR *COS+0%,*AX_2Nminusk,B ; B := B-cos(k*π/N)*RM[k]
STH A,ASM,*AX_k+ ; AR[k] = A/2
STH B,ASM,*AX_k+ ; AI[k] = B/2
NEG B ; B := -B
STH B,ASM,*AX_2Nminusk - ; AI[2N-K] = -AI[k]= B/2
STH A,ASM,*AX_2Nminusk- ; AR[2N-K] = AR
; [k] = A/2
phase4end:
⒊计算所求信号的功率
由于最后所得的FFT数据是一个复数，为了能够方便的在虚拟频谱仪上观察该信号的特征，我们通常对所得的FFT数据进行处理——取其实部和虚部的平方和，即求得该信号的功率。
power:
.asg AR2,AX
.asg AR3,OUTPUT_BUF
pshm st0 ；保存寄存器的值
pshm ar0
pshm bk
STM #d_output_addr,OUTPUT_BUF ；AR3指向输出缓冲地址
STM #K_FFT_SIZE*2-1,BRC ；块循环计数器设置为255
RPTBD power_end-4 ；带延迟方式的重复执行指令
STM #fft_data,AX ；AR2指向AR[0]
SQUR *AX+,A ；A := AR2
SQURA *AX+,A ；A := AR2 + AI2
STH A,7,*OUTPUT_BUF ；将A中的数据存入输出缓冲中，
ANDM #7FFFH,*OUTPUT_BUF+ ；避免输出数据过大在虚拟示波器
；中显示错误
popm bk ；保存各个寄存器值
popm ar0
popm st0
power_end:
RET
注意，在上面的程序中将数据放回输出缓冲准备输出时使用了指令：
STH A,7,*OUTPUT_BUF
对acc A左移7位是为了让显示的数据值在一个合适的范围内有利于观察显示的图形，由于所有的数据都左移了7位，所以从总体上看整个波形的性质还是一样的。同时，由于有的数据太大，为了避免显示数据的溢出而导致在虚拟示波器中观察到的波形错误所以我们使用了指令：
ANDM #7FFFH,*OUTPUT_BUF+
来取出有效的数据位数。
五、实验目的：
在数字信号处理系统中，FFT作为一个非常重要的工具经常使用，甚至成为DSP运算能力的一个考核因素。FFT是一种高效实现离散付氏变换的算法。离散付氏变换的目的是把信号由时域变换到频域，从而可以在频域分析处理信息，得到的结果再由付氏逆变换到时域。本实验的目的在于学习FFT算法，及其在TMS320C54X上的实现，并通过编程掌握C54X 的存储器管理、辅助寄存器的使用、位倒序寻址方式等技巧，同时练习使用CCS的探针和图形工具。另外在BIOS子目录下是一个使用DSP/BIOS工具实现FFT的程序。通过该程序，你可以使用DSP/BIOS提供的分析工具评估FFT代码执行情况。
六、实验内容：
本实验在CCS下完成256点的实数FFT，并通过CCS的图形显示工具观察结果。其主程序为初始化，并通过探针工具读入256点方波数据（在文件fft.dat中，该数据文件可以通过程序fft_data.c修改，但注意数据的绝对值不要超过0x23ff）。FFT的实现由四个子程序代码bit_rev、fft、unpack和power代码完成。实验可以分为以下几步：
（1）启动CCS，在Project菜单相项中打开FFT目录下的fft.mak文件。
（2）用鼠标展开左面项目栏，打开fft.asm源程序。
（3）使用Build选项完成编译、连接，然后使用File菜单中的Load Program将OUT 文件装入
（4）在代码中”wait_input”后的”nop”处和”process”中”b _c_int00”处添加断点
（5）选择View -> Graph -> Time/Frequency菜单打开一个图形显示窗口，两处断点和显示窗口的配置如图所示：
其中输入数据为读入的“FFT.dat”文件，起始位置为0x2300，读入数据总长为256，显示窗口长度也为256，显示窗口的“Start Address”从0x2200开始(6）再打开一个图形显示窗口，显示输入信号。点击Animate按钮，运行程序，观察输出图像。
(7) 右键输入信号的图形显示窗口，设置其显示类型为其fft变换，比较两个窗口图形。
七、实验器材（设备、元器件）：
电脑一台，C54x软件
八、实验步骤：
1．启动CCS，在Project菜单相项中打开FFT目录下的fft.mak文件。用鼠标展开左面项目栏，打开fft.asm源程序。
2．使用Build选项完成编译连接，然后使用File菜单中的Load Program将OUT文件装入3．按照实验内容中所述设置断点，配置参数；并运行程序，观察结果。
4．将输入信号图形显示为其FFT变换，观察结果。
九、实验数据及结果分析：
256点方波输入下的输入信号与其FFT变换的图形如下：
调整输入信号显示窗口，使其进行系统自带的FFT变换之后的图形显示如下：
可以看出，两种FFT变换基本符合。
十、实验结论：
十一、总结及心得体会：
十二、对本实验过程及方法、手段的改进建议：
无
报告评分：指导教师签字：
快速傅里叶变换 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。如果要提高
快速傅里叶变换(FFT)的DSP 实现 (马灿明 计算机学院 计算机应用技术 2110605410) 摘要:本文对快速傅里叶变换(FFT)原理进行简单介绍后，然后介绍FFT 在TMS320C55xx 定 点DSP 上的实现，FFT 算法采用C 语言和汇编混合编程来实现，算法程序利用了CCS 对其结果进行了仿真。 关键字：FFT ，DSP ，比特反转 1.引言 傅里叶变换是将信号从时域变换到频域的一种变换形式，是信号处理领域中一种重要的分析工具。离散傅里叶变换（DFT ）是连续傅里叶变换在离散系统中的表现形式。由于DFT 的计算量很大，因此在很长一段时间内使其应用受到很大的限制。 20世纪60年代由Cooley 和Tukey 提出了快速傅里叶变换（FFT ）算法，它是快速计算DFT 的一种高效方法，可以明显地降低运算量，大大地提高DFT 的运算速度，从而使DFT 在实际中得到了广泛的应用，已成为数字信号处理最为重要的工具之一。 DSP 芯片的出现使FFT 的实现变得更加方便。由于多数的DSP 芯片都能在单指令周期内完成乘法—累加运算，而且还提供了专门的FFT 指令（如实现FFT 算法所必需的比特反转等），使得FFT 算法在DSP 芯片上实现的速度更快。本节首先简要介绍FFT 算法的基本原理，然后介绍FFT 算法的DSP 实现。 2.FFT 算法的简介 快速傅里叶变换（FFT ）是一种高效实现离散傅里叶变换（DFT ）的快速算法，是数字信号处理中最为重要的工具之一，它在声学，语音，电信和信号处理等领域有着广泛的应用。 2.1离散傅里叶变换DFT 对于长度为N 的有限长序列x(n)，它的离散傅里叶变换（DFT ）为 1,1,0, )()(1 0-==∑-=N k W n x k X n n nk N （1） 式中， N j N e W /2π-= ，称为旋转因子或蝶形因子。 从DFT 的定义可以看出，在x(n)为复数序列的情况下，对某个k 值，直接按（1） 式计算X(k) 只需要N 次复数乘法和（N-1）次复数加法。因此，对所有N 个k 值，共需要N 2 次复数乘法和N(N-1)次复数加法。对于一些相当大有N 值（如1024点）来说，直接计算它的DFT 所需要的计算量是很大的，因此DFT 运算的应用受到了很大的限制。 2.2快速傅里叶变换FFT 旋转因子W N 有如下的特性。 。对称性： 2/N k N k N W W +-= 。周期性： N k N k N W W += 利用这些特性，既可以使DFT 中有些项合并，减少了乘法积项，又可以将长序列的DFT
kn N W N N 第四章 快速傅里叶变换 有限长序列可以通过离散傅里叶变换(DFT)将其频域也离散化成有限长 序列.但其计算量太大,很难实时地处理问题,因此引出了快速傅里叶变换 (FFT). 1965 年，Cooley 和 Tukey 提出了计算离散傅里叶变换（DFT ）的快 速算法，将 DFT 的运算量减少了几个数量级。从此，对快速傅里叶变换（FFT ） 算法的研究便不断深入，数字信号处理这门新兴学科也随 FFT 的出现和发 展而迅速发展。根据对序列分解与选取方法的不同而产生了 FFT 的多种算 法，基本算法是基２DIT 和基２DIF 。FFT 在离散傅里叶反变换、线性卷积 和线性相关等方面也有重要应用。 快速傅里叶变换（FFT ）是计算离散傅里叶变换（DFT ）的快速算法。 DFT 的定义式为 N ?1 X (k ) = ∑ x (n )W N R N (k ) n =0 在所有复指数值 W kn 的值全部已算好的情况下，要计算一个 X (k ) 需要 N 次复数乘法和 N －1 次复数加法。算出全部 N 点 X (k ) 共需 N 2 次复数乘法 和 N ( N ? 1) 次复数加法。即计算量是与 N 2 成正比的。 FFT 的基本思想：将大点数的 DFT 分解为若干个小点数 DFT 的组合， 从而减少运算量。 W N 因子具有以下两个特性，可使 DFT 运算量尽量分解为小点数的 DFT 运算： （1） 周期性： ( k + N ) n N = W kn = W ( n + N ) k （2） 对称性：W ( k + N / 2 ) = ?W k N N 利用这两个性质，可以使 DFT 运算中有些项合并，以减少乘法次数。例子： 求当 N ＝4 时，X(2)的值
#include 实验四 快速傅里叶变换（FFT ） 4.1实验目的 1）加深对快速傅里叶变换（FFT ）基本理论的理解； 2）了解使用快速傅里叶变换（FFT ）计算有限长序列和无限长序列信号频谱的方法； 3）掌握用MATLAB 语言进行快速傅里叶变换时常用的子函数。 4.2实验原理 1）用MATLAB 提供的子函数进行快速傅里叶变换 从理论学习可知，DFT 是唯一在时域和频域均为离散序列的变换方法，它适用于有限长序列。尽管这种变换方法是可以用于数值计算的，但如果只是简单的按照定义进行数据处理，当序列长度很大时，则将占用很大的内存空间，运算时间将很长。 快速傅里叶变换是用于DFT 运算的高效运算方法的统称，FFT 只是其中的一种。FFT 主要有时域抽取算法和频域抽取算法，基本思想是将一个长度为N 的序列分解成多个短序列，如基2算法、基4算法等，大大缩短了运算的时间。 MATLAB 中提供了进行快速傅里叶变换（FFT ）的子函数，用fft 计算DFT ，用ifft 计算IDFT 。 2）用FFT 计算有限长序列的频谱 基本概念： 一个序号从1n 到2n 的时域有限长序列()x n ，它的频谱()j X e ω定义为它的离散时间傅里叶变换，且在奈奎斯特（Nyquist ）频率范围内有界并连续。序列的长度为N ，则211N n n =?+。计算()x n 的离散傅里叶变换（DFT ）得到的是()j X e ω的N 个样本点()k j X e ω。其中数字频率为 k 2πω()d ωk k N == 式中：d ω为数字频率的分辨率；k 取对应－(N －1)/2到(N －1)/2区间的整数。 在实际使用中，往往要求计算出信号以模拟频率为横坐标的频谱，此时对应的模拟频率为 s s 2π2π?ω/T ()()T k k k k kD N L ==== 式中：D 为模拟频率的分辨率或频率间隔；T s 为采样信号的周期，Ts ＝1/Fs ；定义信号时域长度L ＝N T s 。 第一章 快速傅里叶变换（FFT ） 4.1 填空题 （1）如果序列)(n x 是一长度为64点的有限长序列)630(≤≤n ，序列)(n h 是一长度为128点 的有限长序列)1270 (≤≤n ，记)()()(n h n x n y *=（线性卷积），则)(n y 为 点的序列，如果 采用基FFT 2算法以快速卷积的方式实现线性卷积，则FFT 的点数至少为 点。 解：64+128-1＝191点; 256 (2)如果一台通用机算计的速度为：平均每次复乘需100s μ，每次复加需20s μ，今用来计算N=1024点的DFT )]([n x 。问直接运算需（ ）时间，用FFT 运算需要（ ）时间。 解：①直接运算：需复数乘法2 N 次，复数加法） （1-N N 次。 直接运算所用计算时间1T 为 s s N N N T 80864.12512580864020110021==?-+?=μ）（ ② 基2FFT 运算：需复数乘法 N N 2log 2 次，复数加法N N 2log 次。 用FFT 计算1024点DTF 所需计算时间2T 为 s s N N N N T 7168.071680020log 100log 2 222==?+?=μ。 (3)快速傅里叶变换是基于对离散傅里叶变换 和利用旋转因子k N j e π2-的 来减少计算量，其特点是 _______、_________和__________。 解：长度逐次变短；周期性；蝶形计算、原位计算、码位倒置 （4）N 点的FFT 的运算量为复乘 、复加 。 解：N N L N mF 2log 2 2== ；N N NL aF 2log == 4.2 选择题 1．在基2DIT —FFT 运算中通过不断地将长序列的DFT 分解成短序列的DFT ，最后达到2点DFT 来降低运算量。若有一个64点的序列进行基2DIT —FFT 运算，需要分解 次，方能完成运算。 A.32 B.6 C.16 D. 8 解：B 2．在基2 DIT —FFT 运算时，需要对输入序列进行倒序，若进行计算的序列点数N=16，倒序前信号点序号为8，则倒序后该信号点的序号为 。 A. 8 B. 16 C. 1 D. 4 解：C 3．在时域抽取FFT 运算中，要对输入信号x(n)的排列顺序进行“扰乱”。在16点FFT 中，原来x(9) 快速傅里叶变换(FFT)的原理及公式 原理及公式 非周期性连续时间信号x(t)的傅里叶变换可以表示为 式中计算出来的是信号x(t)的连续频谱。但是，在实际的控制系统中能够得到的是连续信号x(t)的离散采样值x(nT)。因此需要利用离散信号x(nT)来计算信号x(t)的频谱。 有限长离散信号x(n)，n=0，1，…，N-1的DFT定义为： 可以看出，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算法的优越性。 将x(n)分解为偶数与奇数的两个序列之和，即 x1(n)和x2(n)的长度都是N／2，x1(n)是偶数序列，x2(n)是奇数序列，则 其中X1(k)和X2(k)分别为x1(n)和x2(n)的N／2点DFT。由于X1(k)和X2(k)均以N／2为周期，且WN k+N/2=-WN k，所以X(k)又可表示为： 上式的运算可以用图2表示，根据其形状称之为蝶形运算。依此类推，经过m-1次分解，最后将N点DFT分解为N／2个两点DFT。图3为8点FFT的分解流程。 FFT算法的原理是通过许多小的更加容易进行的变换去实现大规模的变换，降低了运算要求，提高了与运算速度。FFT不是DFT的近似运算，它们完全是等效的。 关于FFT精度的说明： 因为这个变换采用了浮点运算，因此需要足够的精度，以使在出现舍入误差时，结果中的每个组成部分的准确整数值仍是可辨认的。为了FFT的舍入误差，应该允许增加几倍log2(log2N)位的二进制。以256为基数、长度为N字节的数 FFT原理与实现 2010-10-07 21:10:09| 分类：数字信号处理 | 标签：fft dft |举报|字号订阅 在数字信号处理中常常需要用到离散傅立叶变换(DFT)，以获取信号的频域特征。尽管传统的DFT算法能够获取信号频域特征，但是算法计算量大，耗时长，不利于计算机实时对信号进行处理。因此至DFT被发现以来，在很长的一段时间内都不能被应用到实际的工程项目中，直到一种快速的离散傅立叶计算方法——FFT，被发现，离散傅立叶变换才在实际的工程中得到广泛应用。需要强调的是，FFT并不是一种新的频域特征获取方式，而是DFT的一种快速实现算法。本文就FFT的原理以及具体实现过程进行详尽讲解。 DFT计算公式 本文不加推导地直接给出DFT的计算公式： 其中x(n)表示输入的离散数字信号序列，WN为旋转因子，X(k)为输入序列x(n)对应的N个离散频率点的相对幅度。一般情况下，假设x(n)来自于低通采样，采样频率为fs，那么X(k)表示了从-fs/2率开始，频率间隔为fs/N,到fs/2-fs/N截至的N个频率点的相对幅度。因为DFT计算得到的一组离散频率幅度值实际上是在频率轴上从成周期变化的，即X(k+N)=X(k)。因此任意取连续的N个点均可以表示DFT的计算效果，负频率成分比较抽象，难于理解，根据X(k)的周期特性，于是我们又可以认为X(k)表示了从零频率开始，频率间隔为fs/N,到fs-fs/N 截至的N个频率点的相对幅度。 N点DFT的计算量 根据(1)式给出的DFT计算公式，我们可以知道每计算一个频率点X(k)均需要进行N次复数乘法和N-1次复数加法，计算N各点的X(k)共需要N^2次复数乘法和N*(N-1)次复数加法。当x(n)为实数的情况下，计算N点的DFT需要2*N^2次实数乘法，2*N*(N-1)次实数加法。 旋转因子WN的特性 1.WN的对称性 2.WN的周期性 3.WN的可约性 根据以上这些性质，我们可以得到式(5)的一系列有用结果 基-2 FFT算法推导 假设采样序列点数为N=2^L，L为整数，如果不满足这个条件可以人为地添加若干个0以使采样序列点数满足这一要求。首先我们将序列x(n)按照奇偶分为两组如下： 于是根据DFT计算公式(1)有： kn N W N N 第四章 快速傅里叶变换 有限长序列可以通过离散傅里叶变换(DFT)将其频域也离散化成有限长 序列.但其计算量太大,很难实时地处理问题,因此引出了快速傅里叶变换 (FFT). 1965 年，Cooley 和 Tukey 提出了计算离散傅里叶变换（DFT ）的快 速算法，将 DFT 的运算量减少了几个数量级。从此，对快速傅里叶变换（FFT ） 算法的研究便不断深入，数字信号处理这门新兴学科也随 FFT 的出现和发 展而迅速发展。根据对序列分解与选取方法的不同而产生了 FFT 的多种算 法，基本算法是基２DIT 和基２DIF 。FFT 在离散傅里叶反变换、线性卷积 和线性相关等方面也有重要应用。 快速傅里叶变换（FFT ）是计算离散傅里叶变换（DFT ）的快速算法。 DFT 的定义式为 N ?1 X (k ) = ∑ x (n )W N R N (k ) n =0 在所有复指数值 W kn 的值全部已算好的情况下，要计算一个 X (k ) 需要 N 次复数乘法和 N －1 次复数加法。算出全部 N 点 X (k ) 共需 N 2 次复数乘法 和 N ( N ? 1) 次复数加法。即计算量是与 N 2 成正比的。 FFT 的基本思想：将大点数的 DFT 分解为若干个小点数 DFT 的组合， 从而减少运算量。 W N 因子具有以下两个特性，可使 DFT 运算量尽量分解为小点数的 DFT 运算： （1） 周期性： ( k + N ) n N = W kn = W ( n + N ) k （2） 对称性：W ( k + N / 2 ) = ?W k N N 利用这两个性质，可以使 DFT 运算中有些项合并，以减少乘法次数。例子： 求当 N ＝4 时，X(2)的值 快速傅里叶变换的原理及其应用 摘要: 快速傅氏变换（FFT），是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。它对傅氏变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。傅里叶变换的理论与方法在“数理方程”、“线性系统分析”、“信号处理、仿真”等很多学科领域都有着广泛应用,由于计算机只能处理有限长度的离散的序列,所以真正在计算机上运算的是一种离散傅里叶变换. 虽然傅里叶运算在各方面计算中有着 重要的作用，但是它的计算过于复杂，大量的计算对于系统的运算负担过于庞大，使得一些对于耗电量少，运算速度慢的系统对其敬而远之，然而，快速傅里叶变换的产生，使得傅里叶变换大为简化，在不牺牲耗电量的条件下提高了系统的运算速度，增强了系统的综合能力，提高了运算速度，因此快速傅里叶变换在生产和生活中都有着非常重要的作用，对于学习掌握都有着非常大的意义。 关键词：快速傅氏变换；图像处理；matlab 前言： 傅里叶变换在信号处理中具有十分重要的作用，但是基于离散时间的傅里叶变换具有很大的时间复杂度，根据傅里叶变换理论，对一个有限长度且 长度为N的离散信号，做傅里叶变换的时间复杂度为) O，当N很大时τ， (2 N 其实现的时间是相当惊人的（比如当N为4 10 10时，其完成时间为τ8 (τ为计算机的时钟周期）），故其实现难度是相当大的，同时也严重制约了DFT在信号分析中的应用，故需要提出一种快速的且有效的算法来实 现。正是鉴于DFT极其复杂的时间复杂度，1965年..JWCooley 和..JWTukey巧妙地利用 NW因子的周期性和对称性，提出了一个DFT的快速算法，即快速傅里叶变换（FFT），从而使得DFT在信号处理中才得到真正的广泛应用。 傅立叶变化的原理； （1）原理 《测试信号分析及处理》课程作业 快速傅里叶变换 一、程序设计思路 快速傅里叶变换的目的是减少运算量，其用到的方法是分级进行运算。全部计算分解为M 级，其中N M 2log =；在输入序列()i x 中是按码位倒序排列的，输出序列()k X 是按顺序排列；每级包含2N 个蝶形单元，第i 级有i N 2 个群，每个群有12-i 个蝶形单元； 每个蝶形单元都包含乘r N W 和r N W -系数的运算，每个蝶形 单元数据的间隔为12-i ，i 为第i 级； 同一级中各个群的系数W 分布规律完全相同。 将输入序列()i x 按码位倒序排列时，用到的是倒序算法——雷德算法。 自然序排列的二进制数，其下面一个数总比上面的数大1，而倒序二进制数的下面一个数是上面一个数在最高位加1并由高位向低位仅为而得到的。 若已知某数的倒序数是J ，求下一个倒序数，应先判断J 的最高位是否为0，与2 N k =进行比较即可得到结果。如果J k >，说明最高位为0，应把其变成1，即2 N J +，这样就得到倒序数了。如果J k ≤，即J 的最高位为1，将最高位化为0，即2N J -，再判断次高位；与4N k =进行比较，若为0，将其变位1，即4 N J +，即得到倒序数，如果次高位为1，将其化为0，再判断下一位……即从高位到低位依次判断其是否为1，为1将其变位0，若这一位为0，将其变位1，即可得到倒序数。若倒序数小于顺序数，进行换位，否则不变，防治重复交换，变回原数。 注：因为0的倒序数为0，所以可从1开始进行求解。 二、程序设计框图 （1）倒序算法——雷德算法流程图 (2)FFT算法流程 快速傅里叶变换FFT的C语言算法彻底研究LED音乐频谱显示的核心算法就是快速傅里叶变换，FFT的理解和编程还是比较难的，特地撰写此文分享一下研究成果。 一、彻底理解傅里叶变换 快速傅里叶变换（Fast Fourier Transform)是离散傅里叶变换的一种快速算法，简称FFT，通过FFT可以将一个信号从时域变换到频域。模拟信号经过A/D转换变为数字信号的过程称为采样。为保证采样后信号的频谱形状不失真，采样频率必须大于信号中最高频率成分的2倍，这称之为采样定理。假设采样频率为fs，采样点数为N，那么FFT结果就是一个N点的复数，每一个点就对应着一个频率点，某一点n(n 从1开始)表示的频率为：fn=(n-1)*fs/N。举例说明：用1kHz的采样频率采样128点，则FFT结果的128个数据即对应的频率点分别是0，1k/128，2k/128，3k/128，…，127k/128 Hz。这个频率点的幅值为：该点复数的模值除以N/2（n=1时是直流分量，其幅值是该点的模值除以N）。 二、傅里叶变换的C语言编程 1、对于快速傅里叶变换FFT，第一个要解决的问题就是码位倒序。假设一个N 点的输入序列，那么它的序号二进制数位数就是t=log2N.码位倒序要解决两个问题：①将t位二进制数倒序；②将倒序后的两个存储单元进行交换。如果输入序列的自然顺序号i用二进制数表示，例如若最大序号为15，即用4位就可表示n3n2n1n0，则其倒序后j对应的二进制数就是n0n1n2n3，那么怎样才能实现倒序呢？利用C语言的移位功能！ 程序如下，我不多说，看不懂者智商一定在180以下！复数类型定义及其运算#define N 64 //64点 #define log2N 6 //log2N=6 /*复数类型*/ typedef struct { float real; float img; }complex; complex xdata x[N]; //输入序列 /*复数加法*/ complex add(complex a,complex b) { complex c; c.real=a.real+b.real; c.img=a.img+b.img; return c; } /*复数减法*/ complex sub(complex a,complex b) { complex c; c.real=a.real-b.real;实验四 快速傅里叶变换(FFT)
快速傅里叶变换(FFT)试题
快速傅里叶变换(FFT)的原理及公式
快速傅里叶变换FFT原理与实现
详解FFT(快速傅里叶变换FFT
快速傅里叶变换的原理及其应用
快速傅里叶变换(FFT)原理及源程序
快速傅里叶变换fft变换
相关主题
文本预览