DS原理基础及其设计 1

DSP原理与应用PrincipleandApplicationofDSP长沙师范学院信息科学与工程学院08第八章TMS320C6748在数字信号处理中的应用学容DSP应用特点与FFT算法FIR与IIR滤波器设计123习内卡尔曼滤波与工程实现引言——DSP在数字信号处理中的应用特点数字信号处理的核心任务：把输入的采样信号序列通过一定的运算，变换为具有一定频率特性的输出序列，如滤波、变换、增强、压缩等。DSP处理器处理数字信号的突出优势：改进的哈佛结构——程序与数据总线分开，可同时取指与取数；硬件乘法器与深度流水线——单周期完成一次乘累加（MAC）；浮点运算单元——动态范围大、精度高，算法溢出风险小；专用指令与高速缓存——满足实时信号处理的时限要求。引言——DSP在数字信号处理中的应用特点

8.1快速傅里叶变换（FFT）概述快速傅里叶变换（FFT，FastFourierTransform）是计算离散傅里叶变换（DFT，DiscreteFourierTransform）的一种高效算法，广泛应用于频谱分析与滤波等领域。FFT利用DFT中的对称性、周期性和可约性等性质，把运算复杂度从O(N²)降至O(N·log₂N)，从而大幅提高计算效率。按抽取方式分类：按时间抽取法（DIT-FFT，Decimation-In-Time）；按频率抽取法（DIF-FFT，Decimation-In-Frequency）。两者的分解思路互为对偶，运算量相同；本讲以基2按时间抽取（DIT）FFT为主线，讲清分解思想、蝶形运算及其在C6748上的程序实现。8.1快速傅里叶变换（FFT）

8.1.1按时间抽取的基2FFT算法

8.1.1按时间抽取的基2FFT算法基2时间抽选法蝶形运算单元：上支路X₁(k)直通输出；下支路X₂(k)先乘旋转因子W_N^(k)，再与上支路进行“加”“减”运算；一次蝶形运算的运算量：1次复数乘法、2次复数加法。蝶形运算信号流图（1）8点DFT分解为两个4点DFT以长度N=8的序列x(n)为例，按奇偶分组后：•偶序列x(0)、x(2)、x(4)、x(6)作4点DFT，得到X₁(0)~X₁(3)；•奇序列x(1)、x(3)、x(5)、x(7)作4点DFT，得到X₂(0)~X₂(3)。两个4点DFT的输出经4个蝶形（旋转因子W₈⁰、W₈¹、W₈²、W₈³）组合，得到8点DFT的全部输出。8.1.1按时间抽取的基2FFT算法（2）4点DFT进一步分解为2点DFT每个4点DFT仍可继续分解为两个2点DFT：•偶数元素x(0)、x(4)构成一组2点DFT，x(2)、x(6)构成另一组，两组输出经蝶形组合得到X₁(0)~X₁(3)；•奇数元素x(1)、x(5)与x(3)、x(7)同理，经蝶形组合得到X₂(0)~X₂(3)。2点DFT不必再分——它本身就是一次蝶形运算。分解到底，N=8的完整FFT只需3级（log₂8=3）蝶形运算。8.1.1按时间抽取的基2FFT算法8.1.2快速傅里叶逆变换（IFFT）

8.1.3按时间抽取的基2的FFT算法实例【例8-1】任务：实现基2时间抽取快速傅里叶变换（FFT）和快速傅里叶逆变换（IFFT）处理功能。原始信号：包含150Hz和350Hz两种频率分量，信号描述为s(t)=5·sin(2π·150t)+15·sin(2π·350t)采样：输入序列x(n)是对原始信号以Fs=1kHz采样得到的采样序列，抽样点数Tn=1024（2的整数幂，满足基2FFT要求）。处理流程：①产生测试信号并装载为复数序列（虚部置0）；

②生成N/2个旋转因子并按比特位反转排列；

③调用DSPF_sp_cfftr2_dit()完成FFT；④求模得到幅度谱，在150Hz、350Hz处观察谱峰；

⑤调用DSPF_sp_icfftr2_dif()完成IFFT并除以Tn，与原始输入逐点比对（误差小于1e-6判定成功）。例8-1代码分析——TIDSPLIB函数接口（1）cossp()/sinsp()：TI浮点运算数学函数，以弧度表示角度，返回值范围[−1.0,+1.0]。（2）DSPF_sp_cfftr2_dit(float*x,float*w,unsignedshortn)：计算输入数组的时域抽取（DIT）基2FFT。•x——指向复数输入数据的指针，包含N个复数，以连续的实数、虚数对形式排列；•w——指向N/2个已按比特位反转排列的复数旋转因子的指针；•n——FFT的复样本长度，必须是2的幂，且满足n≥32。（3）DSPF_sp_icfftr2_dif(float*x,float*w,unsignedshortn)：计算长度为n（2的幂）的复数序列的频域抽取（DIF）基2IFFT。•x——长度为n的输入/输出序列，实部虚部交错存储；输入为比特位反转排列，输出为正常顺序；•w——指向N/2个比特位反转排列的复数旋转因子；n——复样本长度，为2的幂。注：函数的源码见TI官方参考文档《TMS320C67xDSPLibraryProgrammer'sReferenceGuide》。从频谱分析到滤波器设计FFT解决的是“信号里含有哪些频率成分”的问题——频域分析。新问题：如果只想要其中的某些频率成分、把不需要的成分滤掉，该怎么做？——引出数字滤波器。数字滤波器：把输入的采样序列通过一定的运算变换为具有一定频率特性的输出序列，基本运算单元仍是加法器、延时单元和常数乘法器。•8.2有限冲激响应滤波器（FIR）——无反馈、全局稳定、易实现线性相位；•8.3无限冲激响应滤波器（IIR）——含反馈、阶数低、实时性好；•8.4卡尔曼滤波器（拓展）——从“滤频”走向“状态最优估计”。8.2有限冲激响应滤波器（FIR）

8.2有限冲激响应滤波器（FIR）FIR滤波器的横截型结构如图所示，它由三部分组成：①输入序列x(n)的延时线路，包含N−1个延时单元z⁻¹；②N个乘法器——把x(n)及各延时单元的输出与滤波器系数h(0)~h(N−1)分别相乘；③一个加法器——把N个乘法器的输出累加，得到当前输出y(n)。工程意义：该结构的一次输出就是N次“乘法—累加”（MAC）运算，恰好对应DSP硬件乘法器的强项；配合循环寻址与深度流水线，可在C6748上用极少的时钟周期完成。8.2.1FIR滤波器设计——窗函数法基本思想

窗函数法是一种简单高效的FIR滤波器设计方法，其基本思想是：在时域用一个有限长的窗口函数去截断理想滤波器的无限长脉冲响应。设计三步：①选取一种合适的理想频率选择滤波器——其单位冲激响应h_d(n)为非因果、无限长序列；②用一个有限长的窗口函数w(n)截取理想的冲激响应：h(n)=h_d(n)·w(n)；③得到线性相位、因果的有限长序列h(n)，其频率响应即所要设计的FIR滤波器频率响应。由卷积定理可知：加窗截取后的FIR滤波器频率响应，是理想滤波器频率响应与窗函数频率响应的卷积。与理想特性相比，FIR滤波器频率响应出现失真，主要表现在两个方面：一是过渡带加宽；二是过渡带两侧肩峰附近出现起伏振荡（吉布斯现象）。8.2.1常用的五种窗函数（1）矩形窗（RectangularWindow）：w(n)=1，0≤n≤N−1；（2）汉明窗（HammingWindow）：w(n)=0.54−0.46·cos(2πn/(N−1))；（3）汉宁窗（HannWindow）：w(n)=0.5·[1−cos(2πn/(N−1))]；（4）布莱克曼窗（BlackmanWindow）：w(n)=0.42−0.5·cos(2πn/(N−1))+0.08·cos(4πn/(N−1))；（5）凯泽窗（Kaiser）：由零阶贝塞尔函数I₀(·)构造，含可调参数β，。说明：以上窗函数均关于中心点对称——这是加窗后FIR滤波器能够获得严格线性相位（群延时为(N−1)/2个采样周期）的根本原因，也是FIR相较IIR的一个突出优势。8.2.1五种窗函数的时域序列（N=31）五张图分别给出N=31时矩形窗、汉明窗、汉宁窗、布莱克曼窗和凯泽窗的时域序列图像。•五种窗的序列长度均为31点（有限长），均以中点（第16点）为对称中心；•矩形窗在窗内恒为1，两端“硬截断”；其余四种窗两端平滑衰减到0，故称“平滑窗”；•窗函数越平滑（旁瓣越小），对应的时域波形两端越发平缓——这正是频域旁瓣被压低在时域上的体现。矩形窗汉宁窗汉明窗布莱克曼窗凯泽窗8.2.2例8-2——低通FIR滤波器设计任务【例8-2】任务：设计一个低通滤波器。设计指标：•采样频率Fs=1000Hz；•通带截止频率为150Hz；•阻带截止频率为350Hz；•选用汉宁窗（Hann）；•滤波器阶数为N=31。设计步骤：

①由理想低通滤波器求非因果、无限长的h_d(n)；②用汉宁窗截取得到有限长、因果的h(n)（共31点）；③检查幅度谱的过渡带宽度与阻带衰减是否满足指标。例8-2设计结果(a)理想低通滤波器的离散冲激序列（31个参数）——可见其为非因果、关于中心对称的sinc形序列；(b)用汉宁窗截取后N=31的FIR滤波器冲激序列——两端被平滑压至零；(c)理想低通滤波器的幅度谱——理想矩形特性，过渡带为零；(d)截取后FIR滤波器的幅度谱——过渡带被展宽，通带与阻带出现起伏振荡，与理论分析一致。(a)(b)(c)(d)8.3无限冲激响应滤波器（IIR）

8.3.1IIR滤波器的结构——直接Ⅰ型结构：直接Ⅰ型结构：•第一部分由输入序列及其延时部分构成，M个延时单元的输出经加权（系数b_i）后相加；•第二部分由输出序列及其延时部分构成，N个延时单元的输出经加权（系数−a_j）后相加；•最后将两个和式相加，即得到输出y(n)。

特点：结构直观、与差分方程一一对应，便于理解；但延时单元（存储单元）使用最多，硬件开销最大。8.3.1IIR滤波器的结构——直接Ⅱ型

8.3.1IIR滤波器的结构——级联型

8.3.1IIR滤波器的结构——并联型

8.3.2IIR滤波器的设计（一）：冲激响应不变法

基本思想：使数字滤波器的单位冲激响应序列h(n)逼近模拟滤波器的单位冲激响应h_a(t)——对h_a(t)进行等间隔抽样，使h(n)=h_a(nT)，其中T是采样周期。设计流程（设模拟滤波器系统函数H_a(s)的极点均为单阶，且极点个数大于零点个数）：①用部分分式法将H_a(s)展开；②对式展开式作拉普拉斯逆变换，得到模拟滤波器的冲激响应h_a(t)；③对h_a(t)进行等间隔抽样，得到h(n)；④对h(n)作Z变换，即得数字滤波器的系统函数H(z)。

优点：①能够很好地模仿模拟滤波器的时域特性，使数字滤波器在时域上保持与模拟滤波器相近的性能；②数字频率ω与模拟频率Ω之间呈现线性关系，保证了频率域上性能的一致性。缺点：会产生频率混叠（模拟滤波器频响的周期延拓相互重叠），因此仅适用于低通/带通滤波器。8.3.2IIR滤波器的设计

（二）：双线性变换法

基本思路：通过非线性频率映射，把模拟滤波器的频率响应映射到数字滤波器，同时保持稳定性、彻底避免混叠。两步映射：①频率压缩：把整个s平面的jΩ轴压缩到s₁平面的横带−π/T~π/T之内；②标准变换：再通过标准变换关系把这横带映射到z域的单位圆之内。由此建立起s平面与z平面一一对应的单值关系，彻底消除混叠。③变换常数c的选择（式8-39~8-43）：为使模拟滤波器与数字滤波器在不同频率点处有对应关系，引入待定系数c；若要求低频处有较确切的对应关系（Ω≈Ω₁），通常取c=2/T。8.3.3低通IIR滤波器设计实例【例8-3】任务：从模拟滤波器设计IIR数字滤波器。设计指标：模拟低通滤波器逼近类型：巴特沃斯（Butterworth）型低通逼近滤波器；采样频率为1000Hz；通带截止频率为150Hz；信号频率100Hz，噪声频率450Hz；滤波器阶数N=5。设计流程：确定指标→频率预畸变求模拟域指标→设计巴特沃斯型H_a(s)→双线性变换得H(z)→提取分子NUM与分母DEN系数→在C代码中固化为数组。8.4卡尔曼滤波器（KF）

卡尔曼滤波（KalmanFilter，KF）是一种递归最优估计算法：采用信号与噪声的状态空间模型，利用前一时刻的估计值和当前时刻的观测值，来更新对状态变量的估计，求出当前时刻的估计值。1960年由美国科学家鲁道夫·卡尔曼（RudolfE.Kalman）提出，具有计算高效、估计性能优、递归实时性强等优点，被广泛应用于航空航天、导航和控制、信号处理等诸多领域。典型应用：雷达目标跟踪——预测和估计目标的位置、速度等运动状态；GPS/INS组合导航——利用卡尔曼滤波估计的系统误差状态去修正系统，能够提高组合导航系统位置和姿态的估计精度；信号处理——信道估计、信号去噪，从受噪声干扰的接收信号中恢复出原始信号。与FIR/IIR的区别：FIR/IIR属于“频率选择”型滤波器，按频率范围取舍信号；卡尔曼滤波属于“状态最优估计”型滤波器，按统计模型在噪声中恢复有用信息。8.4.1卡尔曼滤波原理（一）：系统模型与状态预测

8.4.1卡尔曼滤波原理（二）：状态更新

8.4.2卡尔曼滤波算法步骤(1)初始化：给定初始状态估计X_(0|0)和初始均方误差P_(0|0)；(2)计算k时刻的状态向量估计值；(3)计算一步预测均方差矩阵；(4)计算滤波增益矩阵；(5)由最新的量测向量，更新状态估计；(6)更新估计均方误差矩阵；(7)在下一个k+1时刻，重复步骤(2)~(6)。8.4.3一维卡尔曼滤波器【例8-4】任务：实现采样电压信号的卡尔曼滤波。关键点：一维条件下的公式退化（Φ=1，H_k=1，无控制输入项）：•状态一步预测：V_pred=V；•均方差预测：P_pred=P+Q；•滤波增益：K=P_pred/(P_pred+R)；•状态更新：V=V_pred+K·(input−V_pred)；•均方差更新：P=(1−K)·P_pred。/*一维卡尔曼滤波器结构体*/typedefstruct{floatV;/*状态估计，电压值*/floatP;/*预测均方差*/floatQ;/*过程噪声协方差*/floatR;/*测量噪声协方差*/floatK;/*卡尔曼滤波增益*/}Kalman1D;/*卡尔曼滤波器初始化*/voidKalman_Init(Kalman1D*kf_1d,floatinit_V,floatinit_P,floatinit_Q,floatinit_R){kf_1d->V=init_V;kf_1d->P=init_P;kf_1d->Q=init_Q;kf_1d->R=init_R;kf_1d->K=0.0f;}floatKalman_Filter(Kalman1D*kf_1d,floatinput){/*状态一步预测*/floatV_pred=kf_1d->V;/*状态预测*/floatP_pred=kf_1d->P+kf_1d->Q;/*均方差预测*//*状态更新*/