《DSP技术及应用》课件第6章

第6章汇编程序设计6.1算法设计方法6.2

FIR滤波器的设计6.3

IIR滤波器设计6.4快速傅立叶变换的DSP实现

6.1算法设计方法

一个复杂的算法，如果不考虑实时实现，可以考虑调用现有的函数实现。如求y=sin(x)，可以直接调用C语言的sin函数。该函数直接采用浮点运算，其内部用泰勒级数公式实现计算，输入和输出参数也都为浮点数。调用C语言的sin函数不仅简单、快捷，而且精度高，但运算量大，要求较多的运行时间。如果时间不够，可考虑用以下方法改进。6.1.1查表法

如果输入变量与输出变量可以用简单的映射关系来描述，就可以用查表法实现相关函数。如sin函数(这里泛指任意一个非线性函数)，若在它一个周期内先均匀采样N个点，制作成一个表，然后将输入变量x变换成表的偏移地址，即可实现，如图6.1所示。图6.1

sin函数一个周期的采样值(N=16)表的偏移地址的计算方法是，将x=［0~2π)弧度变换到0~N－1之间的整数index。

index=floor(x′)

(6.2)则y=sin(x)近似等价为

y=sin_table(index)(6.3)查表法速度快，但表中两点之间的值需用靠近它们的值代替，故精度不高。要提高精度，有多种方法。如果存储单元数富有余地，可采用加大表的长度的方法提高精度;如果存储单元数不足，可采用线性插值的方法提高精度。设采用线性插值计算x′和index，如图6.2所示。图6.2线性插值示意图则插值点x对应的y值近似等价于

y≈y′=y(index)+［y(index+1)-y(index)］(x′-index)(6.4)式中，y(index+1)-y(index)可单独制作一张表，而不必每次都查表计算。采用更高阶插值的方法可以进一步提高数据精度，但运算量也随之增加。6.1.2数学变换法

C54x的直接运算只有加、减、乘、移位和逻辑运算等，要实现其他运算，需采用数学变换将其他运算变换成

这类运算才可以实现，最常见的数学变换为级数展开法。如正弦用级数展开，取前5项为(6.5)式(6.5)具有迭代性。若系数用Q15格式表示，从最里面的括号开始计算，均为乘加(减)运算，可用C54x的汇编专用指令快速实现。余弦等其他函数也可用类似的方法处理。

正弦函数也可用迭代公式(6.6)计算，但不推荐使用，因为它容易产生累加误差。

sin(nθ)=2cos(θ)sin［(n－1)θ］－sin［(n－2)θ］(6.6)例6.1用查表法产生线性调频信号。

线性调频信号(LinearFrequencyModulation，LFM)的理论公式为(6.7)其中式中,K为调频斜率，f1为起始频率，f2为结束频率，T为持续时间。若用数字系统实现调频信号，需对信号进行离散化处理。设采样频率为fs=1/Ts，Ts为采样周期，则t=nTs，n=0,1,2,…,N－1，N=T/Ts为一个LSF信号周期的采样总点数。另设一个正弦周期的离散相位数为M，离散化式(6.7)，得(6.9)式中，u(n)应是u(nTs)，因Ts为常数，书写习惯上一般省略。公式虽然简单，但运算量却不小。要在一个采样周期内算出一个采样值u(n)，需要足够快的CPU。一般将运算公式变换成运算量小的加减运算，尽量避免运算量大的乘除运算。

令其中(6.10)(6.11)(6.12)当n>0时，f(n)可写成递推公式的形式

f(n)=f(n－1)+2Kcn－Kc+Kp

(6.13)

用两个累加器实现：(6.14)其中φ(0)=0，φ(0)=0。从式(6.7)可以看出，u(n)仍为正弦信号，只要按式(6.13)求出n时刻的f(n)，查正弦表即可求出u(n)(直接计算正弦值的运算量很大，一般不采用)。

如f1=450Hz，f2=150Hz，T=5s，设采样频率fs=8kHz，

一个正弦周期的离散相位数M=1024,可算出Kc=－4.8×10－4,Kp=57.6。

6.2

FIR滤波器的设计

数字滤波器可分为IIR滤波器和FIR滤波器。从设计理

论来讲，设计数字滤波器是比较复杂的过程，但借助于

MATLAB等现代设计软件，可使问题大大简化。如果利用这些工具设计出来的结果不能满足实践要求，则需要根据数字滤波器的设计理论进行计算，再转化为DSP来实现。

FIR滤波器的系统函数为(6.15)对z作反变换，可得其时域输入输出关系为(6.16)这是一个典型的乘加运算，可采用直接型结构实现，如图6.3所示。图6.3

FIR滤波器的直接型实现

FIR的实现结构还有级联型、频率抽样型、快速卷积结构和线性相位FIR滤波器(也称系数对称FIR滤波器)结构等多种

形式。如果系数具有对称性，可采用线性相位FIR滤波器结构，以节省运算时间。在DSPLIB库函数中，一般采用级联型结构，因为级联型结构更容易实现代码的重复使用。

要实现FIR滤波器，需合理使用存储器，以获得较高的运算效率。6.2.1直接型FIR滤波器实现

直接型FIR滤波器延迟单元(z－1)的实现不同，对存储器的要求也不同，主要有线性缓冲区法和循环缓冲区法，用来实现数据延迟和更新。

1.线性缓冲区法

如果在数据存储器定义的缓冲区是双访问随机存储器(DRAM)，在一个周期内既能读一个操作数又能写一个操作数，则可用线性缓冲区法实现。

在程序存储器中定义滤波系数h向量，首地址设为pmad。在数据存储器定义一个缓冲区用于存放x序列值，最新的数据存放在缓冲区顶部，最“老”的数据放在底部，将辅助寄存器指针ARx初始化并指向底部，如图6.4所示。图6.4线性缓冲区法的数据存储设源累加器src初始化为0，pmad指针每计算一次就自动加1，而ARx指针初始化为缓冲区底部，每计算一次就自动减1，且数据缓冲区的数据每计算一次就自动移到下一个地址(高地址)的单元上，其运算如下。第1次运算：src=src+h(N-1)×x［n-(N-1)］，x［n-(N-1)］移出到缓冲区外；

第2次运算：src=src+h(N-2)×x［n-(N-2)］，x［n-(N-2)］移到原x［n-(N-1)］单元；

…

第N次运算：src=src+h(0)×x(n)，x(n)移到原x(n－1)单元；接下来输入x新的采样数据x(n+1)，然后重复上述计算。一次滤波运算需N次乘加运算。C54x有专门的硬件结构和MACD指令来实现算法。

MACD指令格式和执行内容如下：

MACDSmem，pmad，src

；src=src+Smem*pmad,(Smem+1)=Smem,T=Smem其中，pmad绝对寻址在重复执行时会自动加1并指向下一个单元。例6.2

FIR线性缓冲区法滤波子函数示例。

_fir:

PSHMST1;压栈保护ST1

SSBXSXM

SSBXFRCT

SSBXOVM;有符号小数运算,Q15状态

;设当前采样值x(n)已载入到A中，更新x(n)

STM#data_buffer,AR2;data_buffer为数据缓冲区x的首地址

STLA,*AR2;更新x(n)

STM#(data_buffer+N-1),AR2;从高地址开始计算

RPTZ A,#N-1

MACD*AR2-,COFF,A;COFF为系数向量h的首地址SFTA A,-16;仅输出高16位，函数返回值必须存放在A的低16位

RET

2.循环缓冲区法

如果在数据存储器定义的缓冲区是单访问随机存储器(SARAM)，在一个周期内就只能读一个操作数或写一个操作数，那么线性缓冲区法就不能实现，可用循环缓冲区法实现。如图6.5(a)、(b)所示，在数据存储器开辟2个存储块存放滤波系数向量h和输入向量x，设h首地址(低地址，以下同)为h，用ARy指针指示系数表位置，x首地址为x，用ARx指示缓

冲区位置。在滤波运算过程中,x缓冲区的数据位置不变，仅用最新的数据更新最老的数据，如图6.5(c)所示，用x(n+1)更新x［n－(N－1)］，降低了对存储器的要求。图6.5循环缓冲区示意图运算过程如下所述。

初始化：设置循环缓冲大小寄存器BK=R(缓冲区大小)，AR0=1(循环寻址地址增量)，ARy用系数向量的首地址(h)初始化，ARx用输入数据的首地址(x)初始化，且h和x的地址必须

满足循环寻址要求，即它们低N位地址为0，且2N>R，运算之前源累加器初始化为0。第1次运算：src=src+h(0)×x(n)，运算后指针未超边界，ARy→h+1,ARx→x+1；

第2次运算：src=src+h(1)×x(n-1)，运算后指针未超边界，ARy→h+2,ARx→x+2；

…

第N次运算：src=src+h(N-1)×x［n-(N-1)］，运算后指针超边界，ARy→h,ARx→x。

接下来输入x新的采样数据x(n+1)更新x［n-(N-1)］。由于数据指针ARx→x(即当前x(n)所在单元)，因此指针需后退一步才可以更新x［n-(N-1)］，然后重复上述计算。

一次滤波运算需N次乘加运算。C54x有专门的硬件结构和指令来实现MAC算法。MAC指令格式和执行内容如下：MAC *AR2+0%,*AR3+0%,A例6.3

FIR循环缓冲区法滤波程序示例。

.globalstart,_FIR

.mmregs

N.SET 6;

.bss

h,N,1,8；定义一个8字的边界,8为大于6的最小整数次幂

.bssx,N,1,8

.bssy,1

.text

start:

STM#1,AR0

STM#h,AR2

STM#x,AR3

STM#y,AR4

STM#N,BK

CALL_FIR;

B

start;

_FIR:

PSHMST1

SSBXSXM

RSBXFRCT

SSBXOVM

;当前采样值在AL中

STL

A,*AR3;更新采样值

RPTZA,N－1

MAC*AR2+0%,*AR3+0%,A;双间接寻址只能使用AR2~AR5

SFTA

A,－16;若为函数输出，只能放在AL

STL

A,*AR4

MAR*AR3-0%;指针调整

POPM

ST1

RET

.end6.2.2系数对称FIR滤波器实现

线性相位FIR滤波器的系数是对称的，即h(n)=±h(N－1－n)，n=0,1,…,N－1，利用对称性可以进一步提高运算速度。下面以N=8的偶对称为例来说明这类FIR滤波器在C54x

上的实现。展开其输出方程，可得

y(n)=h(0)×［x(n)+x(n－7)］+h(1)×［x(n－1)+x(n－6)］

+h(2)×［x(n－2)+x(n－5)］+h(3)×［x(n－3)+x(n－4)］

C54x有一条专用指令(FIRS指令)来完成系数对称FIR滤波器计算。格式为

FIRSXmem,Ymem,pmad;B+A(32~16)×pmad→B,((Xmem)+(Ymem))<<16→A

;如果重复计算，pmad会自动加1，指令变成单字指令

FIRS可以高效地完成系数对称的FIR滤波运算，但由于指令流水线是先读后写，执行B+A(32~16)×pmad→B时对累加器A为读，其中A(32~16)为原来A中的值，而执行((Xmem)+(Ymem))<<16→A时对累加器A为写，故有一个时间错位，本次计算的A不能在本次B+A(32~16)×pmad→B中使用。如图6.6所示，在程序存储器安排系数表，首地址为pmad。在数据存储器开辟2个缓冲区，大小均为N/2=4个单元，首地址分别定义为x_new和x_old，均需满足循环寻址要求。这里N/2=4，存储器边界应取8,数据的安排如图6.6所示。图6.6

FIRS数据存储示意图运算过程如下所述。

初始化：AR2用x_new初始化,AR3用x_old初始化,AR0=1(地址增量);BK=N/2(循环缓冲区的大小)；

第1次：先计算x(n)+x(n-7)并存储到AH，AR2和AR3循环寻址，地址加1;

第2次：执行FIRS指令，B=B+［x(n)+x(n－7)］×h(0)，x(n－1)+x(n－6)→AH；AR2和AR3同上；第3~5次：与第2次类似，最后把B的结果输出;

第6次：调整AR2和AR3指针，将new缓冲区最老的数

据复制给old缓冲区最老的数据单元,本次将x(n－3)复制到

x(n－7)；

第7次：将最新的数据(本次为x(n+1))载入new缓冲区最老的数据单元(本次为x(n－3))并调整指针，为下次FIR滤波作准备，如图6.7所示。图6.7更新后的FIRS数据存储示意图例6.4系数对称的FIR滤波程序示例。

.globalstart,_FIR,_input

.mmregs

.data

COFF:

.word

1,2,3,4

N.SET4;

.bssx_new,8,1,8;边界必须比N大

.bssx_old,8,1,8

.bss_input,1;数据输入

.bssy,1

.text

start:

STM#1,AR0

STM#x_new,AR2

STM#x_old,AR3

STM#y,AR4

STM#N,BK

_FIRS:

PSHMST1

SSBXSXM

RSBXFRCT

SSBXOVM

;当前采样值在变量_input中

LD*(_input),A

STLA,*AR2;更新采样值x(n)

ADD*AR2+0%,*AR3+0%,A;AH=x(n)+x(n-7)

RPTZB,#N-1

FIRS*AR2+0%,*AR3+0%,COFF

STHB,*AR4

MAR*+AR2(-2)%;指针调整

MAR*+AR3(-1)%

MVDD*AR2,*AR3+0%

;用new最新的数据更新old最老的数据

POPMST1

BFIRS;

.end6.2.3

FIR滤波系数的MATLAB辅助设计

利用MATLAB函数或图形界面fdatool可以方便地计算FIR滤波系数。

(1)函数fir1：用窗函数设计法计算FIR滤波系数，调用格式有以下三种：

B=fir1(N,wc)

B=fir1(N,wc,′highorstop′)

B=fir1(N,wc,win)

B=fir1(N,wc)和B=fir1(N,wc,win)表示设计一个具有线性相位的N阶FIR低通或带通滤波器，截止频率wc为对奈奎斯特频率进行归一化的值。如果是带通，wc为低端和高端截止频率组成的向量。窗函数用变量win指明，默认为汉明窗(hamming)。

如果是高通或带阻，用B=firl(N,wc.′highorstop′)格式,输入参数加字符串‘high’或‘stop’。

(2)函数fir2：用频率抽样法计算FIR滤波系数，调用格式有以下两种：

B=fir2(N,F,A)

B=fir2(N,F,A,win)

用频率抽样法设计一个具有线性相位的N阶FIR滤波器，其频率响应由向量F和A指定，默认为汉明窗(hamming)。如果改用其他窗函数，可用变量win指明。

(3)界面fdatool：利用MATLAB的滤波器设计与分析图形界面(FilterDesignandAnalysisTool)可以快速地设计各类滤波器。使用方法是在MATLAB命令窗输入fdatool命令即可打开设计界面，设置相关选项和输入相关参数，即可得到滤波器参数。6.2.4

DSPLIB库中的FIR滤波函数

DSPLIB库有几个现成的FIR滤波函数,可以直接调用。(1)fir：直接型实FIR滤波，调用格式为

fir(DATA*x,DATA*h,DATA*r,DATA*dbuffer,ushortnh,ushortnx)

算法为变量说明：

*x为输入数据实向量x的指针，x为DATA数据(Q15格式数表示为整型数据)，nx为向量大小，x的存储需满足循环寻址要求；

*h为滤波系数h的指针，h=［b0,b1,b2,…］为DATA数据，nh为向量大小，h的存储需满足循环寻址要求；

*r为滤波输出r的指针，r为DATA数据；

*dbuffer为滤波计算缓冲区，dbuffer的存储需满足循环寻址要求，第一次滤波时需初始化为0。

(2)firdec：直接型带输入抽取运算的实FIR滤波，调用格式为

firdec(DATA*x,DATA*h,DATA*r,DATA**dbuffer,ushortnh,ushortnx,ushortD)

算法为

D为抽取因子，其他变量说明均同fir函数。

(3)firinterp：直接型带输入插值运算的实FIR滤波，调用格式为

firinterp(DATA*x,DATA*h,DATA*r,DATA**dbuffer,ushortnh,ushortnx,ushortI)

算法为

I为抽取因子，其他变量说明均同fir函数。

(4)firs：系数对称的实FIR滤波，调用格式为

firs(DATA*x,DATA*r,DATA**dbuffer,ushortnh2,ushortnx)

算法为

(5)firs2：通用系数对称的实FIR滤波，调用格式为

firs2(DATA*x,DATA*h,DATA*r,DATA**dbuffer,ushortnh2,ushortnx)

算法为

同firs函数相比，系数向量不必存储在程序存储器中，适合多滤波器程序，但执行效率较低。6.2.5调用DSPLIB库文件的方法

设54xdsp.lib库文件在C:＼ti＼c5400＼dsplib下，其相关头文件在C:＼ti＼c5400＼dsplib＼include中,则CCS软件中BuildOptions(构建选项)的设置方法如下。

(1)在Compiler选项卡的Preprocessor中的IncludeSearchPath(-i)中输入头文件地址,如图6.8所示。图6.8设置包含文件搜索途径

(2)在Compiler选项卡的Basic中的LibarySearchPath(-i)中输入库文件地址，IncludeLibraries(-l)中输入工程所包含的库文件名,如图6.9所示。

如果所建工程正在库文件所在的文件夹，上述设置可以简化。

在Preprocessor中的IncludeSearchPath(-i)中输入引号中的符号“..＼..＼”即可，同样，在Basic中的LibarySearchPath(-i)中也是输入“..＼..＼”。图6.9设置库文件搜索途径及包含的库文件名6.3

IIR滤波器设计

6.3.1

IIR滤波器传递函数及实现结构

一个IIR数字滤波器的系统函数为(6.17)式中,一般有M≤N，称为N阶IIR系统。令z=ejω，可得其频域特性为据此可求出其幅度特性和相位特性。(6.1８)式(6.18)的IIR滤波器时域输入输出关系为(6.19)

IR滤波器有多种实现结构，图6.10(a)所示即为一个二阶IIR滤波器直接I型结构，直观明了，但需要4个存储单元存储x和y过去的值。若需要减少存储单元，可使用直接II型结构，如图6.10(b)所示。对于一个二阶IIR滤波器来说，只需2个存储单元。图6.10二阶IIR滤波器的直接I型和II型结构若要避免直接型结构的缺点，可采用级联型或并联型结构,如图6.11所示。它将系统函数分解为若干个二阶系统函数(称为一节)并以级联或并联的形式实现，这里把一阶系统函数看做是二阶的特例。由于各节结构相同，容易实现，同时又相互独立，因此系数调整对其他各阶没有影响。图6.11

IIR滤波器的结构6.3.2

IIR滤波系数的MATLAB辅助设计

例6.5设计一个数字切贝雪夫带通滤波器，给定指标为(1)带内波动δ1≤2dB，当200Hz≤f≤400Hz；

(2)带外衰减δ2≥20dB，当f≤100Hz或f≥600Hz；

(3)抽样频率fs=2kHz。解：(1)求传递函数。

在MATLAB输入：

fs=8e3;

%fs=8kHz

wp=500/(fs/2);

%Hz

ws=1000/(fs/2);%Hz

rp=3;

rs=50;

［n,wn］=buttord(wp,ws,rp,rs);

［b,a］=butter(n,wn);

freqz(b,a,512);画出频率特性运行MATLAB后，可以求出所设计的系统函数为(6.20)其中其频率特性如图6.12所示，其中横坐标为归一化频率，乘以奈奎斯特频率(fs/2)即为实际频率，可以验证其是否符合设计要求。图6.12例6.5的IIR滤波器频率特性

(2)系数归一化。要解决DSP的实现问题最简单的方法

就是直接采用浮点运算，但其运算量较大。C54x中，采

用Q15数据格式作乘加运算，可大大提高运算速度。为此，可考虑将输入x(n)、输出y(n)及系数向量b和a都统一为Q15

格式。因此，对式(6.20)的传递函数作适当变换，分子分母都除以略大于分子分母绝对值最大的数,用MATLAB命令处理如下：

MAX=max(abs(［b,a］))*1.01；

b2=b/MAX;

a2=a/MAX;

处理后的系数向量b2和a2的每个系数绝对值都小于1，可用Q15表示。当输入x(n)和系数向量[WTBZ]b2和a2[WTBZ]都用Q15格式数据表示后，滤波后的输出y(n)是否能保证用Q15格式数据表示，还需要研究，可用MATLAB进行仿真验证。有关MATLAB代码如下：

fc=400;fs=fc*8

n=0:fs*8;

x=sin(2*pi*fc/fs*n)/2;

y2=filter(b2,a2,x);用调整后的系数向量b2和a2计算输出。如果y(n)的每个绝对值都小于或等于1，则没有溢出;如果有溢出，考虑到IIR系统是线性系统，可将输入x(n)都统一缩小一定的幅度，以保证输出不溢出，最简单的方法就是将x(n)的值缩小一半，x(n)算术右移一位。

(3)转化为DSP数据格式。由于CCS开发环境的汇编器

不接受小数数据，需将Q15格式数据转化为整型数据。用MATLAB生成一个系数表，复制到一个空的.inc文件中。

B=round(b2*2∧15);A=round(a2*2∧15);

(4)必要时进行滤波器模型变换。使用级联型IIR结构，可以方便、独立地调整传递函数的零极点，且所用的存储单元最少。在MATLAB中，将传递函数转化为二阶分式级联模型的命令为

［sos,g］=tf2sos(num,den);式中，num、den分别为传递函数分子、分母多项式的系数，sos为二阶分式级联模型的系数，如下所示：对式(6.20)进行函数传递，可求出二阶分式级联模型的系数为同上,对sos进行归一化，即可得到DSP编程所需的系数。6.3.3

DSPLIB库的IIR滤波器程序

根据算法，可以自编汇编程序实现IIR滤波，也可以调用DSPLIB库的IIR滤波函数来实现。IIR滤波函数主要有以下几个函数。

(1)iircas4：二阶级联的IIR滤波函数，每节滤波器采用直接II型，4个系数。

shortiircas4(DATA*x,DATA*h,DATA*r,DATA**dbuffer,

ushortnbiq,ushortnx)

(2)iircas5：二阶级联的IIR滤波函数，每节滤波器采用直接II型，5个系数。

shortiircas5(DATA*x,DATA*h,DATA*r,DATA**dbuffer,ushortnbiq,ushortnx)

(3)iircas51：二阶级联的IIR滤波函数，每节滤波器采用直接I型，5个系数。

shortiircas51(DATA*x,DATA*h,DATA*r,DATA**dbuffer,ushortnbiq,ushortnx)

(4)iir32：二阶级联的IIR滤波函数，32位数据，每节滤波器采用直接II型，5个系数。

iir32(DATA*x,LDATA*h,DATA*r,LDATA**dbuffer,ushortnbiq,ushortnx)

6.4快速傅立叶变换的DSP实现

6.4.1基2复数FFT算法的原理

有限长序列的频谱是周期性的连续函数，在数字频率的一个周期ω=0~2π内，频谱有无数个值。用DFT(离散傅立叶变换)可计算在一个数字频率周期内的N个抽样值，其计算公式为(6.21)式中,数字抽样频点为

FFT是离散傅立叶变换(DFT)的一种快速算法。由于我们在计算DFT时,一次复数乘法需用四次实数乘法和两次实数加法,一次复数加法则需两次实数加法。这样，每计算一个X(k)需要4N次复数乘法和2N+2(N－1)=2(2N－1)次实数加法,因此整个DFT运算总共需要4N2次实数乘法和N×2(2N－1)=4N2－2N次实数加法。设复数序列x(n),n=0,1,…,N－1的长度为N=2L，L为整数，若长度N不是2L形式，则用数字0补足长度。将N=2L的复数序列x(n)按N的奇偶性分成x1(m)和x2(m)两组m=0,1,…,

N/2－1，分别作DFT变换，得X1(k)和X2(k)，也就是说将一个N点的DFT分解成两个N/2点的DFT，然后按式(6.22)重新组合成一个N点的DFT。(6.22)(6.2３)同直接计算DFT的运算量相比，将一个N点的DFT分解成两个N/2点DFT，其运算量大约减少一半。如果继续这个过程，运算量将以指数级下降。一个N点的DFT运算最终分解成N/2个2点的DFT运算，整个运算只需求（N／2）log2N次复数乘法运算，运算量比原来的N2次的复数乘法次数大大降低了(特别是N较大时)，这种算法称为基2复数FFT算法。基2复数FFT算法将一个N点DFT运算每一级都分解成N/2个2点DFT运算，共有log2N级，但每级的加权系数WkN的个数不同。在log2N级(最后一级)，共有N/2个系数WkN,0≤k≤N/

2－1，每往前一级，系数个数减少一半，但系数的指数为前一级系数指数的2倍。第一级时，系数只有1个，为W0N=1，这时复数乘法就可以省去。6.4.2实序列FFT算法(RFFT)原理

当输入为纯实数的时候，可以利用左右对称的特性，更好地计算DFT，称这样的RFFT优化算法为包装算法：首先将2N点实序列的连续输入包装成复数，称为“进包”，然后运行N点复数FFT，最后将N点复数FFT输出并合成2N点实序列的FFT输出，称为“开包装”。

2N点实序列FFT算法如下：(6.24)把x(n)分成偶数点和奇数点，可得到(6.25)

n为偶数时，令n=2r;n为奇数时，令n=2r+1，则有(6.26)将X(k)分解为实部和虚部，即

X(k)=XR(k)+jXI(k)

(6.27)式中(6.28)(6.29)由FFT的性质可知道，实序列的DFT是共轭对称的，即

X(k)=X*(2N－k)

(6.30)

或XR(k)=XR(2N－k),XI(k)=－XI(2N－k)。

下面通过N点复数FFT运算得到上面的X(k)。令N点复序列为y(r)，通过如下方法构造

y(r)=x(2r)+jx(2r+1)，r=0,1,2,…,N－1(6.31)则y(r)的DFTY(k)为(6.32)将Y(k)分解为实部和虚部，得