数字信号处理实验基于matlab

数字信号处理实验基于matlab内容摘要：

般周期信号 将 DT周期信号先展开为 DTFS，再进行逐项变换；      020  mcTXex mdNmdFTnTjmNm md     截断的影响 对无限长时间 DT 信号进行计算时，必须进行截断； 时域截断模型：以窗口函数乘以时间函数    nPnx m2 频域对应卷积：         TT dddd dXWTWX //2  21   TTNW d  i n i n 离散取样函数 N=2M+1 截断效果 图 使单频率展宽，出现主瓣（高 L=2a、宽 4π/L）和旁瓣（高 、宽 2π/L）； 对于有限带宽信号，截断导致带外泄露（能量）和纹波现象； L 越小，上述效应越显著。 矩形窗截断必然导致 Gibbs 现象 图 Nyquist 采样定理 采样：利用冲激串相乘使连续时间信号离散化            mdFTns TmXTXnTttxtx  21 恢复采样信号的条件： 3 带限信号   WforX   0 存在最高频率 2/max Wf  4 采样频率 m a x2/1 fTfs  （ Nyquist rate） 满足上述条件时，可通过频域的低通滤波（截止频率为 2/sf ）分离出原始频谱，恢复信号； 该操作等效于时域的理想内插恢复（在每个采样点插入连续取样函数）。 频率混叠问题 不满足条件 2 时，可能产生部分混叠（高频），低频信号可恢复 ； 不满足条件 1 时，总是产生混叠； 当 CT 信号能量有限时，可以增大 T 使得混叠影响足够小； 对信号进行预先低通滤波是消除混叠的有效方式。 时限带限理论 任何非零连续信号都不可能即为时限又为带限； 利用数字技术处理连续信号时必然需要截断，必然产生误差； 在一定误差范围内，有限能量的连续信号可以近视看作为时限带限信号，并利用数字技术处理。 音频信号的实际处理过程 模拟信号 采样及零阶保持 ADC数字编码信号 数字编码信号 DAC零阶保持信号 低通滤波 模拟信号 第四章 DFT和 FFT频谱的计算 绝大多数信号不能采用解析方式进行频谱分析，只能采用数值计算方法计算频谱； DFT 离散付氏变换 只对有限长度的序列定义，由 DTFT 频谱在一个周期内进行 N 点采样得出 定义式 22     1,....1,01 10/2  NnemXNnx NmNmnjd      1,....1,010/2   NmenxmX NnNmnjd  特点：将有限信号的无限宽频谱压缩到一个周期内（有限宽度） 将离散信号的连续频谱采样为有限序列 设 NjeW /2 则     1,....1,01 10  NnWmXNnx Nmmnd     1,....1,010 NmWnxmX Nnmnd 重要性质 时间序列和频谱序列均以 N为周期；可扩展为周期信号；可以从任何一个周期中计算； 例 DFT和 DTFT的对比计算 DFT（离散序列）是 DTFT（连续频谱）的采样 — 140 例 DFT与 DTFS 的关系： 将有限长时间序列看作无限长时间序列的一个周期，就可以进行 DTFS 变换   dmd NcmX  只相差常数 N 反变换的时域混叠问题 当时域序列为有限长 m 时，作 N 点的 DFT N≥m 不会产生混叠 DFT 和 DTFT 反变换结果相同； Nm 会产生混叠 DFT 和 DTFT 反变换结果不同； 当时域序列为有限长 m 时，作 N 点的 DFT 当时域序列为无限长时，作 N 点的 DFT 的反变换一定会产生混叠；但若时域序列为绝对可和，则当 N 足够大时，混叠可以忽略。 DFT的性质 （与 DTFS 和 DTFT类似） 对于时域序列，将有限长序列扩展为无限长周期序列进行讨论。 奇偶性 周期时移 —圆周移动 DFT解决频率计算问题的方式： 对有限时间序列的连续频谱，只计算其中离散的采样点； DFT存在的问题： 计算量问题： 对于 N 点 DFT， 需进行 N2次复数乘法， N（ N1）次复数加法； N的增加导致运算量庞大， N 的减少导致误差增大； FFT 快速付氏变换 利用 DFT 的奇偶性质及周期性质，对计算公式进行分解，降低计算量。 23     1,....1,010 NmWnxmX nmNnd 要点：  nmNjnm eW /2 以 N 为周期，均匀分布于单位圆上 nmNnm WW  nmNnm WW  2/ 例： N=8 8 点序列的 DFT 分析           mmmd WxWxWxxmX 32 3210          mmmm WxWxWxWx 7654 7654          mmm WxWxWxx 642 6420          mmmm WxWxWxWx 753 7531          mmm WxWxWxx 642 6420  1X         mmmm WxWxWxxW 642 7531  2X     7,....1,021  mmXWmX m 利用 nmW 的性质可以得到    mXmX 11 4     mXmX 22 4  所以可以得出       3,2,1,021  mmXWmXmX md       3,2,1,04 21  mmXWmXmX md 注意：  mX1 是以 0， 2， 4， 6 点构成序列 nx1 的 DFT  mX2 是以 1， 3， 5， 7 点构成序列 nx2 的 DFT 将 N 点序列拆分为两 个 N/2 点序列进行计算： 24 运算结构：蝶形运算（ butterfly equation） 左边 m 只取 4 个值，右边 m 取 8 个值 继续对  mX1 和  mX2 进行分析：             mmm WWxxWxxmX 2441 6240  11X 12X     3,2,1,012211  mmXWmX m 利用 nmW 的性质可以得到    mXmX 1111 2     mXm 1212 2  所以可以得出第 2 层蝶形运算关系：       1,0122111  mmXWmXmX m       1,02 122111  mmXWmXmX m  mX11 是以 0， 4 点构成序列 nx11 的 DFT  mX12 是以 2， 6 点构成序列 nx12 的 DFT 将 N/2 点序列拆分为两个 N/4 点序列进行计算 进一步可以得到第 3 层蝶形运算关系：      40011 xxX       40111 xxX  上述运算过程可以由蝶形运算图表示为： 25 可以看到，在上图中采用了蝶形运算作为运算基本单元，一个蝶形运算的单元如下图所示： 输入 3 个复数 A ， B， C 输出 2 个复数 X1， X2 输入 /输出关系为： X1=A+B*C X2=AB*C 1 个蝶形运算涉及 1 次复数乘法（ 4 次实数乘， 2次实数加）， 2 次复数加（ 4次实数加）。 FFT 的蝶形运算方式可以极大地减少 DFT 的运算工作量；在 的表 N 的增加， DFT 与 FFT 运算量的比较； 在 FFT 的蝶形运算图中可以看出，各层蝶形运算可以构成流水线处理形式，一次蝶形运算完成，数据提供给下一层后，就可以对新的数据进行处理；因此， FFT 对数据的平均处理时间可以压缩到 1 次蝶形运算的时间。 蝶形运算在软件实现时可以将其设置为函数，在运算中调用；在硬件实现时可将其设计为基本功能单元（目前已有专用的 FFT 集成器件）； MATLAB 的 FFT 函数定义 FFT 函数 fft（ x， N） ：计算 N 点时间序列 x的 DFT     1,....1,010 NmWnxmX nmNn 由时间序列 nx 计算频谱序列  mX ：都为 N 点序列，下标排列都为 0N1。 对称区间排布函数 shift（ fft（ x， N））：参见 可以使  NTm /2  的取值位于 ]/,/( TT  的对称区间内； 用 FFT 进行频谱计算 26 有限时间序列的频谱计算 从 0 开始的序列： 直接对该序列进行计算 ( 例 ) 不是从 0 开始的序列： 将该序列进行周期性扩展后，选取从 0 开始的周期进行计算 % program N=4。 T=。 D=2*pi/(N*T)。 %设置序列点数 N，时域采样周期 T，频域采样周期 D x=[2 1 1 1]。 % 给出时间信号 序列 X=fftshift(fft(x,N))。 % 选取对称区间进行 FFT，得出离散频率序列 m=floor((N1)/2):floor((N1)/2)。 % 设置离散频率坐标向量 w=2*pi::2*pi。 % 设置准连续频率坐标范围及分辨率 X1=2exp(j**w)+exp(j*w)+exp(j**w)。 % 根据定义写出连续频率函数 subplot(1,2,1),plot(w,abs(X1),m*D,abs(X),39。 o:39。 ),title(39。 (a)39。 )。 %幅频特性图 subplot(1,2,2),plot(w,angle(X1),m*D,angle(X),39。 o:39。 ),title(39。 (b)39。 )。 %相频特性图 当序列很短时，频域采样点太少，频率分辨率低； 提高频率分辨率的手段： 为时域序列补 0，增加采样点数 N 从 0 开始的序列： 在该序列之后补 0 可以通过直接改写程序中的 N 值实现，当 N 大于 x序列长度时， fft（ x， N）自动为 x补 0； 直接对该序列进行计算 ( 例 ) 不是从 0开始的序列： 将该序列进行尾端补 0后 ，再进行周期性扩展，然后选取从 0开始的周期进行计算；（ 例 ） 若只是计算幅频特性，也可以将序列起点直接移到 0点进行计算（时移不影响幅频特性）； 无限时间序列的频谱计算 无限序列通常不存在分辨率问题，但必须进行截断才能计算； 截断必然导致误差，计算时需要考虑计算精度问题； 方案： 寻找最小的 a，使得采用 N=2a点序列和 N/2 点序列计算之差在指定误差范围（幅频特性峰值的百分比）内； 要点： 只对幅频特性进行比较； 若时间序列为实序列，则幅频特性为偶函数，只对正区间比较； 两序列频谱采样 点密度不同，只能在公共的采样点上进行比较； 例 要点：两序列比较时，长序列密度比短序列大，因此长序列隔位与短序列逐位比较； % program T=。 a=1。 b=100。 beta=1。 while bbeta N1=2^a。 n1=0:N11。 x1=.^n1。 X1=fft(x1)。 N2=2*N1。 n2=0:N21。 x2=.^n2。 X2=fft(x2)。 m1p=0:N1/2。 d=max(abs(abs(X1(m1p+1))abs(X2(2*m1p+1))))。 mm=max(abs(X1(m1p+1)))。 b=d/mm*100。 a=a+1。 end N2,b 连续时间信号的频谱计算 27 当连续时间信号不能表达为闭合形式时，只能采用数值计算； 连续信号的计算必须先在长度为 L 的区间内经时间采样（周期 T）成为 N 点序列，才能进行计算； 本节只考虑绝对可积的连续时间信。