现代数字信号处理5齐齐哈尔大学PPT课件

1、现代数字信号处置齐齐哈尔大学通讯学院何鹏第五章第五章 自顺应滤波器自顺应滤波器60年代以后才出现，开展很快。所谓自顺应DF：利用前一时辰已获得的滤波器参数等结果，自动地调理现时辰的滤波器参数，以顺应信号与噪声未知的或随时间变化的统计特性，从而实现最优滤波。这个概念是从仿生学中引伸出来的，生物能以各种有效的方式顺应生存环境，生命力极强。二最小均方误差LMS自顺应DF：以均方误差最小为准那么，能自动调理单位脉冲呼应h(n),以到达最优滤波的时变最正确DF也即：参数会变，随着外界参数变化自动调理，使滤波器效果最正确。19571966年美国通用公司运用于天线，为了抑制旁瓣而提出。奠定自顺应滤波器的人是

2、：美国B.Windrow 及Hoff：提出自顺应DF算法，主要用于随机信号处置。三、目的设计自顺应DF，可以不用预先知道信号与噪声的自相关函数。在滤波过程中，即使信号与噪声的自相关函数随时间缓慢变化，DF也能自动顺应，自动调理到满足均方误差最小的要求。四、自顺应开展前景11、广泛用于系统模型识别如系统建模：其中自顺应滤波器作为估计未知系统特性的模型。2、通讯信道的自顺应平衡如：高速modem采用信道平衡器：用它补偿信道失真，modem必需经过具有不同频响特性而产生不同失真的信道有效地传送数据，那么要求信号平衡器具有可调系数，据信道特性对这些系数进展优化，以使信道失真的某些量度最小化。又如：数字

3、通讯接纳机：其中自顺应滤波器用于信道识别并提供码间串扰的平衡器。四、自顺应开展前景23、雷达与声纳的波束构成如自顺应天线系统，其中自顺应滤波器用于波束方向控制，并可在波束方向图中提供一个零点以便消除不希望的干扰。4、消除心电图中的电源干扰如：自顺应回波相消器，自顺应噪声对消器：其中自顺应滤波器用于估计并对消预期信号中的噪声分量。5、噪声中信号的滤波、跟踪、谱线加强以及线性预测等。五、目前常见的自顺应滤波器由于设计简单、性能最正确，自顺应DF是目前数字滤波器领域是活泼的分支，也是数字滤波器研讨的热点。主要自顺应滤波器：递推最小RLS滤波器，最小均方LMS滤波器，格型滤波器、无限冲激呼应IIR滤波

4、器。第二节最小均方误差LMS)自顺应DF的根本原理一、均方误差 用统计方法，大量数求平均，提出均方误差最小准那么，即输出信号与进展信号之间误差最小。其定义为：22)( )()(nsnsEnE丈量数据越多，那么越准确。h(n)x(n)=s(n)+w(n)( )(nsny其中s(n)信号可以是随机信号或规那么信号。10)()()()()(Nmmnxnhnxnhny输出：二、自顺应DF根本原理1.自顺应DF的原理框图自顺应数字滤波器参考输入-+d(j)(j)原始输入x(j)y(j)x(j)表示j时辰的参考输入，y(j)表示j时辰的输出呼应;d(j)表示j时辰的原始输入信号，即所期望的输出呼应；(j)

5、为误差信号=d(j)-y(j);2、自顺应DF的原理(1)自顺应DF的h(n)单位脉冲呼应受(j)误差信号控制。(2)根据(j)的值而自动调理，使之适宜下一刻(j+1)的输入x(j+1)，以使输出y(j+1)更接近于所期望的呼应d(j+1),直至均方误差E2 (j)到达最小值.(3)y(j)最正确地逼近d(j)，系统完全顺应了所参与的两个外来信号，即外界环境。留意： x(j)和d(j)两个输入信号可以是确定的，也可以是随机的，可以是平稳的随机过程，也可以是非平稳的随机过程。从图中可见：自顺应DF是由普通DF+相关抵消回路构成。3、ADF实现可以由FIR DF或IIR DF实现。但由于收敛性及稳

6、定性，目前用得多为FIR DF 实现。FIR滤波器构造有：横向型构造(直接型)Transveral Structure)对称横向型构造(Symmetric Transveral Structure)格形构造(Lattice Structure)4、FIR ADF实现 假设FIR DF的单位脉冲呼应长度为，那么其输出为 可见:(1)是个如今或过去输入值的加权和.(2)加权系数就是h(m)。(3)在自顺应DF中，这个加权系数常用符号wj表示,时间用j表示.(4)那么输出可表示为： (5)自顺应DF可变成自顺应线性组合器。10)()()(NnmnxmhnyNiiijxwjy1)()(5、FIR AD

7、F的框图也即自顺应线性组合器自顺应算法.x1jx2jxNj+-d(j)(j)y(j)w1w2wN 假设设x1j, x2j , x3j xNj ,为同一信号的不同延时组成的延时线抽头方式，即所谓横向FIR构造。它是最常见的一种自顺应DF构造方式。 普通来讲x1j, x2j , x3j xNj , 可以是恣意一组输入信号，并不一定要求当时x1j = xj, x2j= x(j-1),x3j= x(j-2) ，xNj= x(j-N+1) ，即并不要求各xi(j)是由同一信号的不同延时组成.6、横向FIR ADF的构造自顺应算法.x(j)x(j-1) x(j-N+1)+-d(j)(j)y(j)w1w2w

8、N 假设设x(j), x(j-1) , x(j-2) x(j-N+1)j ,为同一信号的不同延时组成的延时线抽头方式，即为横向FIR构造。它是最常见的一种自顺应DF构造方式。 AFy(j)(j)x(j)简化符号为7、由横向FIR AF组成的自顺应系统.x1(j)x2(j)xN(j)+-d(j)(j)y(j)w1w2wNAFAFAF 当所处置的输入信号x1(j), x2(j) , x3(j) xN(j)来自不同的信号源时，它实践上就等于自顺应线性组合器。三、寻觅E2(j)=min时的各wi值自顺应DF的关键在于按照(j)和各xi(j)的值，经过某种算法寻觅出E2(j)=min时的各wi值，从而可

9、自动地调理各wi值。1.写出均方误差的式子首先我们推导出自顺应线性组合器均方误差E2(j)与加权系数wi的关系式。 )()()()()()()()()(,321)()()()(2110jXWjdjyjdjjxjxjxjXwwwWjXWWjXjxWjyTNTNnTii求均方误差：此处大写代表矩阵式中：写成矩阵形式：2.x(j)信号与d(j)信号的自相关函数的自相关函数是的自相关矩阵称为输入信号）（令)()0()()()0()0(0)0()0()0()0()0()()(221222111211221222111221jdjdEjxRxxxxxxxxxxxxxxEjXjXERddxNxNxNxxNx

10、xxxxxNxxxxxjNjNjjNjjjjjjNjjjT3.x(j)信号与d(j)信号的相互关函数0)()0(),0(),0()()()()()()()()(12121mmmjxjdjxjdjxjdEjXjdEPdxTxNddxdxN为时间差，同一时间相关函数，与期待输出的信号的互为输入信号在随机过程中令：4.求出E2(j)与wi的关系 单个值单个值单个值行方阵列列行单个值均方误差为：代入式中求得，将)0( 2)()()0()()()()()()()()()(2)()()()()(2)()()()(2222222ddTTddTTTTTTTWRWWPjdEjEjdEjXjXERjXjdEPWj

11、XjXWEWjXjdEjdEWjXjXWjXWjdjdEjXWjdEjE5.求出自顺应滤波器的E2(j)与wi的关系11111122211122)(2)()0()()0(2)0()0()(1)0(2)0()0( 2)()(NiNmNidixxmiddxdxxddNiNmNidxixxmiddTTiWmiWWjEFIRWWjENWWWWRWWPjdEjEimi数字滤波器有横向不难证明，对于自适应时，只有一个信号当由于均方误差为：第二节性能函数E2(j)及其梯度一、研讨E 2(j)与W的关系 11111122)(2)()0(2)()(NiNmNidixxmiddTTiWmiWWWRWWPjdEjE

12、由于均方误差为：看出：均方误差E 2(j)是加权系数W的二次函数，它是一个中间上凹的超抛物形曲面，是具有独一最小值的函数。二、E 2(j)与W的关系曲线)(2jE)(22jE)(12jEAB调理加权系数W使均方误差最小，相当于沿超抛物形曲面下降到最小值。w)(2jW)(1jW三、梯度法在数学上，可用梯度法沿着该曲面调理权矢量 的各元素 得到这个均方误差E 2(j)的最小值。1.均方误差梯度将对上式 均方误差对权矢量的各wi进展求导，得到均方误差梯度：NdwjdEdwjdEj)()()(212)(2jE2.求最正确权矢量用w*表示 1对均方误差梯度求导 00 , 0 , 12012)(00 ,

13、0 , 100 , 0 , 1010100 , 0 , 10120)(0)(1212WRPdwjdERRRWRRWRWRWWRPdwjdEjTTTTTTTTT为对称方阵，即：求最佳权矢量，则令2求出均方误差梯度表示式 2222)(22)(22)(22)(00 , 0 , 12012,222211211PWRWRPjWRxddwjdEWRxddwjdEWRxddwjdEWRWRxdPxdxdPNjNjNjjjjjjTjNjjjT均方误差梯度为：以此类推：个元素第第二元素第一元素第一元素3维纳-霍夫方程霍夫方程这是著名的维纳或达最佳值。为最小，时，当均方误差梯度为：)(0)( 2 2 2 2)(1

14、*2PRWPWRWjEjPWRWRPj(4)最小均方误差算法,)()(*221*WPRPWjdEjEPRW即可以求出已知最小误差为：最小权矢量为： 实践上，设计自顺应DF无需知道R和P。自顺应DF与维纳平稳随机过程DF比较，其差别在于添加了一个识别控制环节，将输出y(j)与所期望的呼应d(j)比较，看能否一样，假设有误差(j)，用(j)去控制w,使w为E2(j)=min时的W*. 因此，关键：找到LMS算法，寻觅一个W的递推式，由W=W0，起始值开场，沿着趋于W*的正确方向逐渐递推，直至W=W*，E2(j)=min为止。这就是最小均方误差算法，简称LMS算法。第三节LMS递推算法寻觅一个W的递

15、推式，由W=W0，起始值开场，沿着趋于W*的正确方向逐渐递推，直至W=W*，E2(j)=min为止一、LMS算法递推式LMS递推算法是Windrow与Hoff两个提出的。设w(j)是j时辰的权矢量，w(j+1)是j+1时辰的权矢量；那么LMS算法的递推公式为：式中0, 是一个控制稳定性与收敛速度的参数。由于E2(j)是权矢量W的二次方程，即E2(j)与W的关系在几何上是一个“碗形的多维曲面。时刻的均方误差梯度。是jj :)()()()1(jjWjW指引正确的递推方向。导向参数是,:)(jj二、自顺应过程的物理意义)(2jE)(22jE)(12jEAB)(2jW)(1jW点。即碗底：。，去寻找“

16、碗”的底点连续地调节，据递推式自适应的物理意义：根*2, 0)()()() 1(wwdwjdEWjjWjWW 为了简单，设W是一维的，那么E2(j)与W的关系成为一个抛物线。三、自顺应递推算法的递推过程1、步骤1右边。必在则时，当左边。必在则时，设当WjWdwjdEjWWWjWdwjdEjWWjwwjww)(, 0)()()(, 0)()()(22)(21212.步骤2.) 1() 1)(, 0)()()() 1(, 0) 122*2*2)(2211*12WjWjWWjWWWjWdwjdEjWWWjWjWwWjWWjww（应为，更接近于（值为了使下一个右边。必在即时，如果设，更接近（值为了使下

17、一个3.步骤3-合并TNjWWdWjdEdWjdEdWjdEjjWjjWjWjWjWjWWWjWdWjdEjWjW)()()()()()()()1()(.)()1()(,0)()()1(,22212*)(2，可用矩阵表示：是多维的情况时，梯度当表示：用梯度更接近于比值都能使下一个左边或右边，在不论原来式中得：将上两式合并4.步骤4-结论*2*)()1(0)()()1()()()()()1(WjWjWjjjWWWLMSjEjWWjjWjW，当确定。靠扰的步距由的方向向碗底靠扰。将沿着时，根据上式，当。算法也称为最陡下降法因此下降最快的方向。代表因为这一点。可以找到由：四、LMS自顺应滤波器递推公

18、式1LMS算法如何实时处置及实现的无偏估计是，得出的数学期望代入上式，得：得又的估值，有：作为均方误差的梯度（用单样本怎样处理呢？参数，不能实时计算，是个集合平均的，可知：由：)()()()()()()(2)()()(,)(,)(),()()()(,)(,)()(2)(,)(,)()()()()()()() 1(21)(21)(222122jjjjEjjXjjjXdWjddWjddWjdjXWjdjdWjddWjddWjdjdWjddWjddWjdjjjjjjjWjWTNTTjWWNTjWWN2 LMS自顺应滤波器递推公式的结构图。也可由此设计硬件实现型，软件实现的编程数学模它们就是自适应）可

19、按上式递推计算，（，）给定输入（可任意设置。）（）起始条件（式中：自适应滤波器递推公式得出一组代入方程中将上式：由误差方程：DFjxjxjxWWWLMSjXjjWjWjXjjWjWjXjjjjWjWNNNNN3)(),()(2)0(,),0(,01.)()(2)() 1()()(2)() 1()()(2)()()() 1(2121111五、自顺应滤波器的主要结论1低。误差最大，估计精度最时不相关时，互相关系数与当)()(, 00)()()()(2min21*jdEjEPRWjXjdEPjdjX 2)()(22WRWWPjdEjETT因为均方误差为：五、自顺应滤波器的主要结论2高。误差减少，估计

20、精度提则相关时，与当)()()()()(2*221*jdEPWjdEjEPRWjdjXT 2)()(22WRWWPjdEjETT因为均方误差为：五、自顺应滤波器的主要结论3高。，最准确，估计精度最均方误差为则时，即阶数完全相关时，与当00)(0)()()()()()(1)()(1)()(21*jEjyjdjjdjxjyPRWRPjdjXNjdjX 2)()(22WRWWPjdEjETT因为均方误差为：五、自顺应滤波器的主要结论4自顺应数字滤波器是个线性系统，时变，服从叠加原理。第四节 自顺应数字滤波器的运用 自顺应滤波器最重要特性：能有效地在未知环境中跟踪时变的输入信号，使输出信号到达最优。因

21、此在电信，雷达，声纳，实时控制以及图象处置等领域都有胜利的运用。一、运用引见自顺应数字滤波器的运用非常广泛，这里引见四种。 一、自顺应噪声抵消器； 二、自顺应陷波滤波器； 三、自顺应预测系统。二自顺应噪声抵消器1、自顺应噪声抵消器引入 固定参数的数字滤波器利用本身的传输特性来抑制信号中的干扰成分，消除干扰的效果遭到很大的限制。假设知道干扰信号的来源，就可利用干扰源的输出去抵消信号中的混杂的干扰。但直接利用干扰源的输出去抵消干扰的做法是危险的，由于由于延迟的影响，不仅不能减小信号中的干扰，反而有能够使干扰加强。在自顺应噪声抵消器中，是利用干扰源的输出，经过一个数字滤波器，最正确地估计出干扰值，从

22、而从混有干扰的输入中减去干扰估值，实现了干扰与信号相当完善的分别。2、自顺应噪声抵消器的原理框图信号源噪声源自顺应滤波+-)(ne)( ns)(0nV)(1nV)()(0nVnS原始输入互不相关。、与且程，为零均值的平稳随机过、设定误差最小。最佳估值，即两者均方为相关的噪声参考输入为与图中：原始输入为：)()()()()()(:)()()()()()(101000100nVnVnSnVnVnSnVnVnVnVnVnS3、自顺应噪声抵消器的输出1噪声抵消。是相关，这样才能进行、要求。时，。误差最小求自适应滤波是使均方无关，、与程为零均值的平稳随机过、出自适应噪声抵消器的输)()(min)(0)(

23、)(min)(0)()()()()()(,)()()()()()(2)()()()()()()(2)()()()()()()()()(00220020010100020022002002200nVnVneEnVnVEneEnVnVnsEnVnVnSnVnVnSnVnVnsEnVnVEnSEneEnVnVnSnVnVnSnenVnVnSnene3、自顺应噪声抵消器的输出2相等则最佳。（一般为噪声）相关，必须与被抵消信号参考输入信号务条件：完成自适应噪声抵消任看出：自适应滤波器要。这时自适应滤波器关闭号，什么信号，出来什么信也表明：原始信号进来那么：不相关，则、若)()(000)()()()()(

24、2)()()()(0)()()()(0120202002020220000nVnVnVEnVEnSEnVnVEnVEnVEnSEneEnVnVEnVnV4、自顺应噪声抵消器的运用 1胎心心电图中消去母体干扰1由于胎儿心电图的研讨，处理妇产科难胎位、单胎位、双胎位、分娩期间心率能否正常及优生学方面孕妇怀孕的中，后期预测胎儿在的生理情况。胎儿心电只需母体的1/10.胎儿心电图在胎儿腹壁丈量。称腹壁胎儿心电图，简称心电图。1胎心心电图中消去母体干扰2输出腹部电极原始输入胸参部输电考极入自顺应胎儿心电图测试仪其中原始输入a(t)=f(t)+m(t)+n(t)f(t):胎儿心脏产生信号m(t):母亲心脏

25、产生信号n(t):噪声干扰信号主要由肌肉起的，有时称“肌肉噪声。 采用自顺应噪声抵消器消除胎儿心电图中母体心脏信号干扰。普通采用：四个普通胸导每路信号一样记录母亲心跳，作为参考输入信号。经过自顺应噪声抵消器处置后，母亲心脏干扰信号被显著消弱，胎儿心声可辨。1胎心心电图中消去母体干扰3预滤波带宽335Hz.胎儿脉冲谐波很多抽样频率：为256Hz.2噪声抵消器其他运用语音信号的镇噪、飞机、汽车，船舱内大量噪声的抑制。天线旁瓣干扰的消除，以及消除50Hz纹波等。三、自顺应陷波器假设信号中的干扰是单频的正弦波，设频率为w0那么消除这种干扰的正确方法是运用陷波器。1、陷波器理想频率特性)(jHw02、自

26、顺应陷波器的优点与普通陷波器比较，有两 大优点：1可以自顺应地准确跟踪干扰频率。2容易控制带宽，且值越高，带宽越窄。为带宽。，其中单频干扰频率为QffffQ00,3、单频干扰陷波器框图90LMS算法+-参考输入原始输入)( jjy)cos(0twC)cos()(0tAts)( jd)(1jx)(2jx频率相同相位与干扰频率不同但幅度参考输入陷波要去掉信号其中：原始输入：)cos(:)cos()(00tctAts4.多麦克风降噪参考麦克风输入主麦克风输入噪声源MATLAB 处置同窗们好模型S 信号n0 噪声,它是n1经过滤波以后得到的噪声,即n0 =filter(w, n1)n1 参考噪声y n

27、1 经过滤波以后得到的信号(滤波器w的系数未知)假设: 它们都是零均值的, S与n0及n1不相关, 但是n0 和n1相关Singal + Noise0FilterNoise1s+n0n1yRecovered Signalz+-主麦克风输入参考麦克风输入模型目的: 对n1进展自顺应滤波,滤波器为w, 从主麦克风得到的信号中减去y, 得到 “干净的信号z,从而到达降噪 的目的。Singal + Noise0FilterNoise1s+n0n1yRecovered Signalz+-主麦克风输入参考麦克风输入算法)(20nyE 20min22minnyEsEzE代价函数战略战略: 调整滤波器系数调整

28、滤波器系数, 使使 最小最小, 这也意味着这也意味着 最小最小. 这个意义下这个意义下, 可以看成是误差可以看成是误差 信号信号z2zE 2202022nyEsEnysEsEzE滤波器系数的调整滤波器的阶数是worder, 滤波器的输出其中调整滤波器系数的LMS算法Tinworderini)(,),1()(111n)()()(1iiiyTnw)()() 1()(1iziiinww滤波器系数的调整调整滤波器系数的RLS算法)()() 1()(iziiikww)() 1()(1)() 1()(11111iiiiiiTnPnnPk) 1()()() 1()(111iiiiiTPnkPP语音00.51

29、1.522.533.544.55x 104-0.6-0.4-0.6signal噪声环境中的语音00.511.522.533.544.55x 104-10-50510primary microphone signalLMS滤波效果00.511.522.533.544.55x 104-4-2024primary microphone signal00.511.522.533.544.55x 104-101filtered output0123456x 1040102030Mean Squared Error主麦克风信号滤波输出信号均方误差RLS滤波效果00.511.522.53

30、3.544.55x 104-505primary microphone signal00.511.522.533.544.55x 104-0.500.5filtered output0123456x 104024Mean Squared Error主麦克风信号滤波输出信号均方误差线性调频信号0200400600800100012001400160018002000-1-0.8-0.6-0.4-0.60.81signalLMS滤波效果主麦克风信号滤波输出信号均方误差0200400600800100012001400160018002000-505020040060080010

31、0012001400160018002000-1010500100015002000250001020RLS滤波效果主麦克风信号滤波输出信号均方误差0200400600800100012001400160018002000-5050200400600800100012001400160018002000-101050010001500200025000510第六章第六章 现代谱估计现代谱估计6.1 离散随机过程与非参数化谱估计离散随机过程与非参数化谱估计6.1.1 离散随机过程离散随机过程x(n)是延续随机过程是延续随机过程x(t)的均匀采样的均匀采样.121211221212(,) ()()

32、()()()()(,)()()xxxxxxxxxxRnnEx nx nCEx nnx nnRnnnn 先采样再计算还是先计算再采样？ 离散随机过程x(n)是广义平稳的，假设其均值为常数，自相关函数只取决于时间 kn1-n2 ，即2( ) ( ) ()( ) |xxxxxRkE x n x nkCk 离散随机过程x(n)和y(n)是广义结合平稳的，假设二者均为广义平稳，且相互关函数只取决于时间 kn1-n2 ，即2( ) ( ) ()( ) |xyxyxRkE x n y nkCk 平稳离散随机过程x(n)的功率谱为自相关函数的Fourier变换，即2( )( ),jkTxxxxkPRk eT周

33、期为 其反变换 Fourier 级数为1( )( ),/2jkTxxxxRkPedT ( )( )jkTxyxykPRk e 互功率谱可定义为6.1.2 非参数化功率谱估计直接法先计算Fourier变换和间接法先计算自相关函数。1直接法12201011( )|( )|( )|( )( )( )NjnTxNnNjnTNnPXx n eNNXx n ex n（是否为整个序列的幅值谱密度?）2间接法1*01( )()( ),0,1,.,( )( )NxxnMjkTxxxkMRkx nk x n kMMNNPRk e ，11201220( )1( )|( ) ( )|11| ( )|( )|2NjnT

34、xnNnx nPx n c n eNWWc nCdNN 由于将视为周期函数(幅值谱离散，功率谱连续)，所以称为周期图法，得到的功率谱估计为有偏估计。 为了减小偏差，可以采用窗函数对周期图进行平滑。 第一种窗函数直接加给样本数据，修正后的周期图为 另一种( )( ) ( )MjkTxxBTkMPRk w k e窗函数是加给样本自相关函数(Blackman-Tukey法)，功率谱为 6.2 平稳ARMA过程假设离散随机过程xt服从线性差分方程11( )( )()()qpijijx ne nbx nia x njZ变换后可得A(z)X(z)=B(z) E(z).1111( )1.( )1.ppqqA

35、 za za zB zbzb z 记为ARMA(p,q).定义 6.2.1(因果稳定性) 一个由A(z)X(z)=B(z)E(z)定义的ARMA过程是因果稳定的，或称x(n)是e(n)的因果稳定函数，假设存在一常数序列满足以下条件：00|,( )()iiiihx nhe ni且定理 6.2.1 令x(n)是一个A(z)和B(z)无公共零点的ARMA(p,q)过程，那么x(n)是因果稳定的，当且仅当A(z)的零点都在单位圆内.定义 6.2.2(可逆性) 一个由A(z)X(z)=B(z)E(z)定义的ARMA过程是可逆的，假设存在一常数序列i满足以下条件：00|,( )()iiiie nx ni且

36、定理 6.2.2 令x(n)是一个A(z)和B(z)无公共零点的ARMA(p,q)过程，那么x(n)是可逆的，当且仅当B(z)的零点都在单位圆内.定理 6.2.3 假设对一切|z|1有A(z) 0，那么ARMA过程A(z)X(z)=B(z)E(z)具有独一的平稳解1( )( )(),( )iiiiiB zx nhe ni hhzA z由确定定理 6.2.4 任何一个具有有限方差的ARMA或MA过程都可以表示为独一的、阶数能够为无穷大的AR过程；同样，任何一个ARMA或AR过程也可以表示为一个阶数能够为无穷大的MA过程。6.3 平稳ARMA过程的功率谱密度6.3.1 ARMA过程的功率谱密度定理

37、 6.3.1 令y(n)是一具有零均值的离散平稳随机过程，其功率谱密度为Py ()。假设x(n)是由2( )()|, ( )( ) |()|( ),()()|jiiiiijxyjjjiiz eix nh y nihhx nPH ePH eeH eh z 其中 是绝对可求和的，即则也是一个零均值的离散平稳过程，且其功率谱密度由下式给出：式中，是一个的多项式：定理 6.3.2 令x(n)是一个满足差分方程11( )( )()()qpijijx ne nbe nia x nj的平稳ARMA(p,q)过程，e为均值为零方差为2的白噪声，那么x(n)的功率谱密度为2211()( )()( )1(), (

38、 )1()j Txj TpqjijiB ePA eA za x nj B zbx ni ( )( )( )( )( )WoldPPPPP分解：一平稳过程的功率谱密度式中为有理式功率谱密度，为线谱。6.3.2 功率谱等价 功率谱等价不同ARMA模型得到的信号能够具有一样的功率谱。 设ARMA模型为：11211(1) ( )(1) ( ), ( )(0,)| 1,1()| 1,()| 1,1()| 1,()pqijijiiiizx nze n e nNprqsirririssiq 且 个极点中有 个位于单位圆内，其余位于单位圆外；个零点中有 个位于单位圆内，其余的在单位圆外。因果部分反因果部分最小

39、相位部分最大相位部分 ARMA过程x(n)的功率谱密度为：2222121_1111_1111_*| (1) |( )()( )| (1) |ARMA( )( )( ) ( )( )(1)(1)( )(1)(1)1/,iqjiixpjkkpriiiirprkkkksieB zPA zeA z x nB z e nA zzzB zzz 将所有在单位圆外的零、极点全部反演到单位圆内，得到一新的过程式中其中，_*1,;1/,1, .kkirpksq222221122111112 ( )|(1)|(1)|( )( )( )|(1)|(1)|(1)| |(1)| |1| |1|(1( )xxqsjjiii

40、i sprjjkkkk rjjjjjkkkkkkARMAx neeB zPA zeeeeeeeP 新的过程的功率谱密度为21 21 211121 21 2111)|( )|(1)|qqqjiiiii si sxpppjkkkkk rk rePe 结论：假设系统是非因果的或者是非最小相位的，利用功率谱密度，只能辨识出|H(ej)|，而不能辨识出H(ej).可利用互功率谱密度或高阶矩统计量辨识此类系统。6.4 ARMA谱估计问题:利用N个知的观测数据x(0),x(1),x(N-1)估计出ARMA过程x(n)的功率谱密度。直接运用式(3.3.6)估计时，需求辨识出整个ARMA模型及鼓励噪声的方差。M

41、A参数的估计需求解非线性方程。6.4.1 ARMA功率谱估计的两种线性方法1.Cadzow谱估计 问题：知条件是什么？112x110( )()( )()( )( )()( )()( )( )piiiB z B zN zN zPzA z A zA zA zN zpN zn z其中是一个 阶多项式，11100( )()( )()( )( )( )()0.5( ),0( )()( ),kkkxkkkxxN zN zzzA zA zzC k zk zk zC k kkkC kx是的多项式，是 的多项式.利用功率谱密度是自协方差函数的傅立叶变换，可得P其它00000(3.4.2)(3.4.5)( )(

42、)( )(),0,1,. ()piikipiiiipiiipiin zN zk zA za za zakikpk将与相比较，可令两边同乘以，可得n卷积关系优点：防止了MA阶数和参数确实定以及鼓励白噪声方差的估计。2. Kaveh谱估计1-2x11-21-( ) ()( )( )( ) ()( ) ()( ) ()qkkkqlxlqkkkqc zB z B zP zC l zA z A zA z A zB z B zc z 由第二个等式可得，1-( ) ()( )kkqklkxkqlccc zA z A zC l z 左端为实数，因此有由的第三个等式，可得问题：知条件是什么？00-x2i 1(

43、-),0,1,.,( )1jwppkijxijqkkkqpiiz ecaa C k ij kqKavehARMAc zPaz谱估计为6.4.2 修正Yule-Walker方程 在Cadzow谱估计和Kaveh谱估计中，都需求知AR阶数和参数。00000 ( )( )( ) ()( ) ( ) ()( ) ()( ) ()( ) ( )ixikikARMAx nx nh i e niRE x n x nEh i e nih k e nkh i h k 因果过程具有唯一的平稳解其相关函数() ()E e ni e nk2200( )() ()0,( )( ) ()ARMA()xipinie nki

44、E e ni e nkRh i h iah nib由于是白噪声，故有，其它由过程的定义，可得200020()( )()( )ppiiikik lka R lih kah klih k b 由可得1( , )0,(3.4.18)( )()0,(3.4.19)Yule-Walker(3.2.1)AR(3.4.19)( )(ipxixixixARMA p qMAbiqlqR la R lilqpR la R l 对于一个过程，其参数，因此当时，式恒等于零，有以上方程称为修正方程。由式可以给出更简洁、更直观的理解。对于一个( )过程,可简化为1)0,0(3.4.20)Yule-Walkerpiil 以

45、上方程称为方程。定理6.4.1 (AR参数的可辨识性) 假设ARMA(p,q)模型的多项式A(z)和B(z)无对消因子,且ap0，那么该模型的参数a， ap可由以下p个修正Yule-Walker方程确定 1()( ),1,.p(3.4.21)(1)(2) .( )(2)(3) .(1)(3.4.22 )( )(1)(1).pixxixxxxxxxxxaR l iR l lqqR qpR qpR qRR qpR qpR qaR qR qR qp 证明:令1,.,(3.4.22 )(1),.,()(3.4.22 )Yule-Walker-pxxaaabrR qR qpcRar则修正方程可以写成(3

46、.4.23)(1)(2).( )(2)(3).(1)( )(1).(1)(1)(2).()xxxxxxxxxxxxR qpR qpR qR qpR qpR qRR qR qR qpR qR qR qp 令112-21121-112,1,.,1()( ),1,.,1(2)(ppppppppipixxixxaaaaaaRppuipu RliRllq qqpRqpRqR rrrrrrrr111由 修 正 Yule-Walker方 程 可 得进 一 步 可 得的 秩 由 其 前行 决 定 。假 定 前行 线 性 相 关 ， 则 存 在 不 全 为 零 的， 满 足令3).(1)( )(1).(1)(1

47、)(2).()xxxxxxxpRqRqRqRqpRqRqRqp2-11111211111 1,1,.,1()( ),1,.,1(1)( )(1)(2)(1)(2)ipixxipxpxpxppxxxRv ipvR liR l lqqpuR qvuR qv uR qvuR qpu R qpvu R qpv 假设 不满秩，则存在不全为零的，满足2 11 112111-11(3)( )()(1)(2)(1)1()(1)(1),1,.,1()( ),1,2,.(3.xpxxxxpxxxpxipixxiu R qpvu R qR qpv R qpv R qpvR qpR qpu R qpuR qu ipu

48、 R liR l lqq上式右端最后一行等于前行之和，因此有重复上述过程可得存在不全为零的满足与4.19)0,paR比较，可得与题设矛盾，故矩阵 可逆。 由定理6.4.1可知，当ARMA(p,q)过程x(n)的AR阶数p和自相关函数Rx()知时，只需求解p个Yule-Walker方程，便可辨识出AR参数。 然而在实践运用中， AR阶数p也是未知的，因此还需先确定p. 3.4.1(1)( ).(1)(2)(1).(2)(3.4.25 )()(1)().,().xexexeeexexexeexexexeeeeeeeeR qR qR qpRR qR qR qpaR qMR qMR qMpMp pp

49、qqqpqprank Rp 命题 令若且 则6.4.3确定AR阶数的奇特值分解方法 11223.4.3 A,(3.4.36)Hhhm nm mUn nVAA UVm n定理 是一个复数矩阵，则存在一个酉矩阵 和酉矩阵使得 可以分解为是一个对角阵，其主对角线元素是非负实数，并按下列顺序排列(3.4.37)min( , ).hm n11111/2( )222211_111),1(3.4.41).12)/, 1kkkFhhFFkkkkAAkkhAAAFrobeniusTkTkAkh 两种确定矩阵 有效秩的方法：定义（ ）为矩阵 的范数定义略小于 的数 作为阈值，（ ）的最小整数 定义为矩阵 的有效秩

50、。定义1_1(3.4.42)kkTTkA定义接近于零的正数 作为阈值，的最大正数 定义为矩阵 的有效秩。6.5 ARMA模型辨识 在一些运用中，不仅希望得到 AR阶数参数，而且还希望获得MA阶数和参数。6.5.1 MA阶数确定1( )()0(3.5.2).pxixiMAR laR l il 阶数是使成立的最大整数上述方法虽然简单，但是在数据比较短的情况下，数值稳定性不好。6.5.2 MA参数估计0222020 111202(3.4.10)()()(3.5.8)211.qqqqqbbcbbb bcbbcqq由可得上述非线性方程组共有个未知数，但是只有个方程，可以假定 00101,0,1,(3.5

51、.9)(1) 1,(3.5.10 ),qkii kkiTqTqfbbckqqbbbafffbf定义拟合误差函数及维向量0000101(3.5.10 )(1)(1)(3.5.11)qqqqqbqqfffbbbfffbbbTfFb维矩阵=01011010(1)( )( )1( )(0)(0000Newton-Raphson3.5.2 Newton-Raphson3.4.12MA,0,1, ,qqqqqiiiikibbbbbbbbbbbbckqbcbFbbFf=利用方法，可得算法（算法）(1)利用式()计算谱系数并令初始值及0)0,1,iq；( )( )(1)(2)3.5.9,0,1,3.5.12(3)3.5.13M A(4)M A1,ikiifkqiiFb由 式 ()计 算 拟 合 误 差 函 数并 用 式 ()计 算；利 用 式 ()更 新参 数 估 计 向 量；检 验参 数 估 计 向 量 是 否 收 敛 ， 若 收 敛 则 停 止 ，否 则转 (1).实验二实验二 ：( )0.2 (1)0.23 (2)0.06 (3)( )0.7 (1)0.12 (2)x nx nx nx ne ne ne n的ARMA过程。 根本要求：将系统视为灰箱可以利用方程的构造信息，p=3, q=2，根据系统的输入输出关系确定当e(n)为N(0,1)