现代信号处理论文_AR模型的功率谱估计BURG算法的分析与仿真

AR模型谱估计方法研究及其应用摘要数字信号处理(DSP)重要的应用领域之一，是建立在周期信号和随机信号基础上的功率谱估计。

在实际应用中往往不能获得具体信号的表达式,需要根据有限的数据样本来获得较好的谱估计效果，因而谱估计被广泛的应用于各种信号处理中[1]。

本论文研究了功率谱估计的几种常用的方法，包括经典谱估计和现代谱估计的各种方法，且对每种方法的估计质量做了数学推导，并给出仿真程序及仿真图。

经典法主要包括周期图法、自相关法，但这两种方法都存在缺陷，即认为观测数据之外的数据都为零，所以对经典法中的周期图法进行了加窗、平均等修正，因此提出了周期图法的改进方法；现代谱估计的方法分类比较多，AR模型法,MA 模型法和ARMA模型法是现代功率谱估计中最主要的参数模型，本论文着重讨论了AR模型参数法[2]。

同时论文将通过对经典谱估计和现代谱估计的实现方法及仿真图的比较，得出经典功率谱估计方法的方差性较差，分辨率较低，而现代谱估计的目标正是在于努力改善谱估计的分辨率，因此能得到较好的谱估计效果，为此应用更为广泛[3]。

关键字：数字信号处理，功率谱估计，周期图法，自相关法，AR模型法ABSTRACTDigital signal processing (DSP) important application of one of the field. Actually, we can’t get the expression of a specific signal, so we need to estimate the power spectral of a signal according to some sample data sequences.So spectrum estimation which is widely used in various signal processing.In this thesis, some common methods of Power Spectral Estimation, such as classical spectral estimation and modern spectral estimation, are studied. The quality of each estimation method is derived, simulation program and simulation figure is given. Classical methods of Power Spectral Estimation mainly include the Periodogram and the BT method. But both of them have a common drawback: the data sequences, beyond the area of the observed sequences, are all presumed to zero. So the Windows and the average method are introduced to improve the quality of the Periodogram. Therefore the improvement of The Periodogram estimation method is proposed. The classification of modern spectral estimation methods are more , AR，MA, and ARMA is the most important parameters of modern spectral estimation. This thesis will focus on discussion of AR model parameters method. At the same time , It can be seen from the comparison and realization of classical spectral estimation and modern spectral estimation, classical power spectrum estimation variance is poor, low resolution .The goal of modern spectral estimation is working to improve the resolution of spectral estimation, better results of the estimation of the power spectrum can be obtained, so it is applied more widely.Keywords: digital signal processing, Power Spectrum Estimation, The Periodogram, the BT methods，AR model目录摘要 (I)1 绪论 (1)1.1功率谱简介 (1)1.2经典谱估计 (2)1.3现代谱估计 (3)1.4功率谱估计应用及用途 (4)2 谱估计简介 (5)2.1随机信号简介 (5)2.2平稳随机信号 (7)2.3估计质量的评价标准 (10)3 现代谱估计 (12)3.1平稳随机信号的参数模型 (12)3.2 AR模型的正则方程与参数计算 (13)3.3 MA模型谱估计 (16)3.4 ARMA模型谱估计 (17)3.5 AR模型功率谱估计实验 (19)3.5.1、实验内容 (19)3.5.2、实验分析 (20)3.5.3、实验结果及分析 (20)3.5.4、实验思考 (24)3.5.5、实验源代码 (25)3.6 AR模型的应用 (30)3.7 小结 (36)4论文总结 (36)参考文献 (38)1 绪论1.1功率谱简介1.功率谱估计技术渊源流长，在过去的几十年获得了飞速的发展。

AR 模型的功率谱估计BURG 算法的分析与仿真钱平（信号与信息处理 S101904010）一．引言现代谱估计法主要以随机过程的参数模型为基础，也可以称其为参数模型方法或简称模型方法。

现代谱估计技术的研究和应用主要起始于20世纪60年代，在分辨率的可靠性和滤波性能方面有较大进步。

目前，现代谱估计研究侧重于一维谱分析，其他如多维谱估计、多通道谱估计、高阶谱估计等的研究正在兴起，特别是双谱和三谱估计的研究受到重视，人们希望这些新方法能在提取信息、估计相位和描述非线性等方面获得更多的应用。

现代谱估计从方法上大致可分为参数模型谱估计和非参数模型谱估计两种。

基于参数建摸的功率谱估计是现代功率谱估计的重要内容，其目的就是为了改善功率谱估计的频率分辨率，它主要包括AR 模型、MA 模型、ARMA 模型，其中基于AR 模型的功率谱估计是现代功率谱估计中最常用的一种方法，这是因为AR 模型参数的精确估计可以通过解一组线性方程求得，而对于MA 和ARMA 模型功率谱估计来说，其参数的精确估计需要解一组高阶的非线性方程。

在利用AR 模型进行功率谱估计时，必须计算出AR 模型的参数和激励白噪声序列的方差。

这些参数的提取算法主要包括自相关法、Burg 算法、协方差法、 改进的协方差法，以及最大似然估计法。

本章主要针对采用AR 模型的两种方法：Levinson-Durbin 递推算法、Burg 递推算法。

实际中，数字信号的功率谱只能用所得的有限次记录的有限长数据来予以估计，这就产生了功率谱估计这一研究领域。

功率谱的估计大致可分为经典功率谱估计和现代功率谱估计，针对经典谱估计的分辨率低和方差性能不好等问题提出了现代谱估计，AR 模型谱估计就是现代谱估计常用的方法之一。

信号的频谱分析是研究信号特性的重要手段之一，通常是求其功率谱来进行频谱分析。

功率谱反映了随机信号各频率成份功率能量的分布情况，可以揭示信号中隐含的周期性及靠得很近的谱峰等有用信息，在许多领域都发挥了重要作用。

实验报告2——AR过程的线性建模与功率谱估计一、实验目的1.理解AR过程的产生机理，复习实验1估计自相关序列的方法。

2.利用估计出的自相关序列来求解信号的功率谱，即用周期图法来估计功率谱。

3.分别采用自相关法（Yule-Walker法），协方差法，Burg法，修正协方差法来估计功率谱，并与周期图法进行比较，分析性能孰优孰劣。

4.学习matlab在数字信号处理中的应用。

二、实验过程和分析AR过程的线性建模与功率谱估计。

考虑AR过程：=-+-+-+-+()(1)(1)(2)(2)(3)(3)(4)(4)(0)()x n a x n a x n a x n a x n b v n v n是单位方差白噪声。

()(a) 取b(0)=1, a(1)=2.7607, a(2)=-3.8106, a(3)=2.6535, a(4)=-0.9238，产生x(n)的N=256个样点。

用randn(1,N)产生单位方差高斯白噪声v(n)，用v(n)激励滤波器产生AR(4)过程，即用x=filter(b,a,v)产生x(n)，b是滤波器分子系数，这里为b(0)=1，a是滤波器分母系数，a=[1 -2.7607 3.8106 -2.6535 0.9238]。

(b) 计算其自相关序列的估计ˆ()r k，并与真实的自相关序列值相比较。

x结论：真实自相关序列与估计出的序列对比如上图所示，两者在100点之前的形状相似，幅度有一定差异，而且估计出的自相关序列有较大的波动，这是因为估计的点数较少，使得估计精度不够，另外，估计自相关序列的下标越大，用来估计的点数就越少，因而后面的估计值是很不精确的。

因为假设n>= 500 处的数据为0，所有估计自相关不精确，误差较大。

(c) 将ˆ()x rk 的DTFT 作为x (n )的功率谱估计，即： 1211ˆˆ()()|()|N j jk j xx k N P e rk e X e Nωωω--=-+==∑。

毕业设计（论文）题目：基于BURG算法的谱估计研究及其MATLAB实现目录一、毕业设计（论文）开题报告二、毕业设计（论文）外文资料翻译及原文三、学生“毕业论文（论文）计划、进度、检查及落实表”四、实习鉴定表XX大学XX学院毕业设计（论文）开题报告题目：基于BURG算法的谱估计研究及其MATLAB实现机电系电子信息工程专业学号：学生：指导教师：（职称：讲师）（职称：）XXXX年XX月X日英文原文Parametric spectral estimation on a single FPGAABSTRACTParametric, model based, spectral estimation techniques can offer increased frequency resolution over conventional short-term fast Fourier transform methods, overcoming limitations caused by the windowing of sampled, time domain, input data. However, parametric techniques are significantly more computationally demanding than the Fourier based methods and require a wider range of arithmetic functionality; for example, operations such as division and square-root are often necessary. These arithmetic processes exhibit communication bottleneck and their hardware implementation can be inefficient when used in conjunction with multipliers. A programmable, bit-serial, multiplier/divider, which overcomes the bottleneck problems by using a data interleaving scheme, is introduced in this paper. This interleaved processor is used to show how the parametric Modified Covariance spectral estimator can be efficiently routed on a field programmable gate array for real-time applications.1. INTRODUCTIONDue to its ease of hardware and software implementation the shortterm fast Fourier transform(STFFT)is widely used for spectral estimation and is known as the conventional method. However, the technique has drawbacks in terms of spectral resolution and accuracy caused by the finite length of the input data sequence used. Windowing of input data causes spectral broadening and Gibb’s phenomenon of spectral leakage can mask the weaker frequency components of the true power spectral density(PSD)[1]. These unwanted effects can be reduced by using longer data sequence lengths, so that the transformed signal becomes a better representation of the infinite data sequence, but in real life this usually is not feasible as the characteristics of the input data may change with time. Over short periods of time the data signals can often be assumed to exhibit wide-sense stationarity, where the signal characteristics are assumed approximately constant but the spectral resolution is therefore limited. In attempts to improve the PSD estimation, windowing functions, Bartlett or Hanning for example, can be used to reduce side-lobe levels but these lower spectral resolution by broadening the main lobe of the PSD[2].Model based, parametric spectral estimation techniques can alternatively be used, where the unrealistic assumption that data is zero outside the window of interest is dropped[1]. Either knowledge of the underlying process or reasonableassumptions about the nature of the unobserved data are used to improve frequencyresolution over the conventional approaches. The computational burden of suchprocessors is however much higher than the STFFT and arithmetic functions such asdivision and square-root often become necessary. In the division and square-rootnon-restoring algorithms there is an inherent dependency that the result bits mustbe computed in a most significant bit(MSB)first manner, with the computation of abit directly dependent upon the result of the previous one[3]. This interdependencymakes it difficult to efficiently realize such arithmetic functions in hardware,and implementations are usually much slower than other basic functions such asmultiplication, addition and subtraction. Communication bottlenecks can thereforeeasily occur in systolic arrays where different types of processors are interconnected.The difficulties with hardware implementation of parametric spectral estimatorshave led to a preference of software implementation on homogeneous DSP networks[4].However, high levels of processing capacity have not been fully reflected in systemthroughput since the increased communication incurred as a result of parallelismis constrained by communication bus performance. This restricts the range ofproblems that can be computed in realtime and the software approach may sometimesbe inadequate for real-time spectral estimation.In this paper, hardware implementation of a parametric spectral estimator isaddressed. A bit-serial processor capable of division and inner product stepcomputation is developed by combining separate processors for these functions. Thedesign uses a high level of pipelining so that division can be computed at a highrate and multiplication is performed on a MSB first data stream, eliminating thebottleneck problem. The high level of pipelining allows many independentcomputations to be performed simultaneously or interleaved. The use of theinterleaving scheme is demonstrated by implementing the design of a ModifiedCovariance type of parametric spectral estimator, to produce a field programmablegate array(FPGA)based system for the spectral analysis of Doppler signals fromultrasonic blood flow detectors.2. MODIFIED COVARIANCE SPECTRAL ESTIMATIONThe model order p=4 Modified Covariance(MC)spectral estimator, proven to beoptimally cost efficient for the blood flow application where mean velocity and flowdisturbance are of interest[5], involves solving the following linear system ofcovariance matrix equations:(1) where each element j i C,is obtained from:∑∑-=--=+•++-•--=110,)()1(21N p n p N n ji j n i n j n i n N C(2) for a window of length N data samples. The k aˆ filter parameter estimates are obtained by solution of the linear system(1), using the Cholesky, forward elimination and back substitution algorithms. The signal white noise varianceestimate, 2ˆσis calculated as: k p k k c a c ,010,02ˆˆ•+=∑=σ (3)and the power spectral density(PSD), )(ˆnMC f P , is obtained from: 2122221ˆ)(ˆ)(ˆ∑=-=+==p k f j k k n n MC ne z z af A f P πσσ(4)Hence, the MC spectral estimator may be partitioned onto four different programming modules:•CMR-calculation of the elements of the covariance matrix and right-hand side vector, 5N multiply accumulates taking into account matrix symmetry.•Cholesky-solution of the linear system of equations, 6 divisions and 10 inner step products for non-square-root Cholesky, 4 divisions and 12 inner step products for solving triangular systems.•WNV-calculation of the white noise variance, 4 multiply accumulates.•PSD-computation of the power spectral density, 4N inner step products for a zero padded DFT, N multiplications to find absolute value of DFT and N/2 divisions for the PSD.The number of samples, over the fixed time duration window of 10ms, is required to be either 64, 128, 256 or 512 depending on Doppler signal conditions. Implementation of the algorithm in Matlab software proved to be in excess of a factor of 310times too slow for real-time operation and that a performance of up to 13.5 MFLOPS/s is required[4]. Execution times of MC algorithm implementation using various topologies of Texas Instruments TMS320C40 DSPs with T8 transputers as routers have also fallen short of the real-time requirements, where processing time is over 150ms too long in the worst case[4][6]. Use of a single DSP in a PC hosted system has been shown sufficient for the smaller N but the specification of N=512 could not be achieved[4], thus prompting consideration of the hardware approach.3. BIT-SERIAL INTERLEAVED PROCESSORStudy of word-parallel systolic implementations of the MC method has shown the method to provide more than adequate throughput for the specified real-time blood flow application but the cost of such a system is very high in terms of arithmetic units, communication burden and control complexity[7]. For example, a systolic array processor for non-square root Cholesky decomposition[8] requires 13 processing elements(PEs), each PE having 2 to 6 ports of either m(single precision)or 2m(double precision)lines, and control is necessary to reverse data streams before back substitution. An alternative way to approach the hardware design involves consideration of bit-serial processing techniques.The nature of multiplication algorithms normally involve the computation of least significant bits(LSBs)first and bit-serial multipliers reflect this in their output ordering. Conversely, division algorithms such as non-restoring are MSB first in nature[9]. Computation of each quotient bit can be performed from m controlled add subtract(CAS)operations, the decision on whether to add or subtract being taken given the result of the previous bit computed(except on the first operation where the signs of the input operands are used to decide). Allowing carries to ripple through therefore leads to a propagation delay greater than m CAS cells. In a bit-serial multiplier, the delay between successive bits being output is likely to be around a single full adder(FA)delay, leading to a maximum clock frequency approximately m times higher and a communication bottleneck with the divider. The clockrate of the divider can be increased to a similar maximum rate as the multiplier by pipelining the carries in each individual CAS stage. However, this means that each output bit is then available only once in every m clock cycles. There is also the problem that data streams must be reversed between multipliers and dividers. One possibility is to use registers and extra control logic to reorder the bit stream from the divider but the operation time is still limited.The efficiency of the divider with the pipelined carry can be greatly improved by using the redundant slots between the output of successive bits to perform other separate divisions. The bit-serial/word-parallel divider shown in[3]allows m+1 individual divisions to be performed simultaneously or interleaved. This decreases the mean division operation time to achieve similar performance to a bit-serial multiplier but there is still the problem of data stream matching when interfacing such devices. One way to tackle this problem is to redesign the multiplier so that it works on a MSB first data stream, rather than storing and reordering the divider outputs which increases latency and control requirement[10]. MSB first multiplication, first demonstrated by McCanny et al.[11], shows it is possible to perform multiplication on positive numbers by summing partial products(PPs)inreverse order to the norm. This also requires inclusion of an MSB first addition unit to ensure that output carries from the PPs are added into the final product. Larsson-Edefors and Marnane[12], extend the concept of MSB first multiplication to the two’s complement number system and show bit-serial architectures for this application. In order to match the divider bit-streams exactly to the multiplier bit-streams it is then just a matter of inserting extra delays along the FA sum pipeline so that the addition of PPs from a number of different multiplications can be performed simultaneously as shown by Bellis et al[13].Study of the bit-serial interleaved divider and multiplier reveals that both architectures show a large degree of similarity. Both work in load/operational phases; the loading networks for the divisor and multiplier both consist of m+1 delay feedback SISO registers and the FA sum/carry pipelines are alike. Both designs also require MSB first, half adder(HA)cell, addition stages; the divider requires m PEs, for 1’s complement error correction which occurs for negati ve dividends, and the multiplier requires m-1 HA PEs to add the output carries from the PPs. Therefore, it is possible to combine the two designs to make a programmable bitserial device which allows m+1 computations to be simultaneously interleaved, as shown in figure 1.The processor has two mode selection inputs DIVi and SUBi, which control four modes of operation i i Y Z X /0±= or i i i Y X Z Z •±=0 where i Z and 0Z are both double precision. Ldi is the load/operational mode select signal for the storage of i Y and i Z over the first m(m+1) clock cycles. Ldi switches into operational mode over the next m(m+1) clock cycles where the remaining data is input and the bulk of the computation is performed in the FA array. All control signals are fully pipelined similarly to the data, allowing the shortest possible block pipeline period of 2m(m+1) clock cycles and continuous input/output of data(i.e. while one block set of m+1 computations are being output, the next block set may be loaded in). The pipeline also allows independent functionality between each of the separate interleaves and on the same interleave a division may immediately follow an inner step product computation and vice-versa.4. INTERLEAVED PROCESSOR BASED MODIFIED COVARIANCE SYSTEMCost-bene fit analysis on systolic array implementation of the CMR and Cholesky sections of the MC spectral estimator shows that a 12 bit fixed point word-length is suf ficient for these computations[7]. Using the bit-serial processor with a 12 bit word-length results in the capacity for interleaving 13 computations.On interleaves 0 to 4 the CMR multiplications are performed over N consecutive block sets, such that the products i n n x x +• are produced on interleave)40(≤≤i i andblockset )10(-≤≤N n n . A bit-serial systolic array provides the correct input data sequencing from consecutive Doppler signal samples and a separate MSB first double precision accumulator, whose architecture is similar to that of HA section in figure1, computes the covariance matrix elements, which are then stored in RAM. The system for computing the CMR calculation is shown in figure 2.The entire Cholesky, forward elimination, back substitution and WNV computations are performed on interleave 5 on the system shown in figure 3. Here division and inner product step computation are necessary. Once the covariance matrix elements are stored in the dual port RAM after block set N the Cholesky decomposition can commence on interleave 5 while in parallel the CMR computation on the next set of data can be processed on interleaves 0 to 4. A ROM block controls the addressing of the dual port RAM for retrieval of stored data to go onto the processor inputs and storage of the processor results. To achieve good dynamic resolution for the low word-length used, a systolic array scaling module is included between the RAM and the processor, whose scaling factors are also produced by the ROM controller along with the mode control. Overall timing in the system is controlled by three counters, qi(range 0 to 12),qb(range 0 to 23) and qw(range 0 to N)corresponding to the interleaves, bit-position and input word.A zero padded point DFT is computed on interleaves 6, 7, 8 and 9. This is basically amatrix vector multiplication and is computed by using the processor in inner product step mode. The system for this section consists of a ROM to provide storage of the twiddle factor matrix n W , another ROM to control the addressing of the twiddle factors for a particular qw and 4 registers which continuously recirculate the filterparameter results(n aˆ)from the Cholesky decomposition stage. On interleave 6 thereal and imaginary parts of the first set of products 1ˆa W i N • are alternately formed.Using a single flip-flop delay the results of these computations are then fed backinto the i Z input of the interleaved processor to be added to the products 2ˆaW i N • and the DFT is built up in this way. The dynamicrange of the PSD computation is quite high compared to the rest of the system, therefore, at this stage a floating point representation of the DFT results is taken using a systolic based conversion circuit. PIPO registers are used to store the 6 bit exponents of the real and imaginary parts of the DFT, whose squares are computed on interleave 10. On interleave 11 the absolute value of the DFT is computed. The maximum of each pair of real and imaginary results from interleave 10 is fed to the i Z input while the other value is piped into the i Y to be appropriately scaled by the difference in the two squared exponents appearing on the i X input. The PSD is then computed on interleave 12, involving N/2 divisions of the WNV formed on interleave 5 with the absolute values from interleave 11. The exponents of the PSD are then easily derived from the exponents of the DFT results.5. CONCLUSIONThis paper has proposed a bit-serial interleaved processor which can be programmed for use in division or inner product step computations. The interleaving idea was introduced in order to perform bit-serial division at the same high clock rate as multiplication without resorting to carry look-ahead schemes to remove the communication bottleneck. The result is a high throughput processor which is cost ef ficient in terms of VLSI implementation, since communication between PEs in the linear array is localised and control is very simple. An application in parametric spectral estimation, namely implementation of theModi fied Covariance spectral estimator, which makes full use of the interleaving scheme, was described. This system has been programmed and simulated using VHDL. Synthesis was targeted to exploit the resources of a Xilinx XC4036EX-2 FPGA. This type of FPGA has dual port RAM capability, where a 16x1 bit dual port RAM can be implemented in a single con figurable logic block (CLB). A dual port RAM cell is an area ef ficient method to implement a 13 or 14 bit SISO register, as used in the interleaving process. Such registers would otherwise have to be implemented using the pairs of flip flops in each CLB, i.e.7 CLBs. CLBs can also be con figured as ROM blocks which are useful for generating the address signals in the Cholesky and PSD modules, and for storage of the DFT twiddle factors. The processor design exhibits mostly localised communication to make use of the fast routing resources between nearest neighbours in the FPGA’s CLB matrix and enable high speed operation. Timing analysis of the FPGA layout shows that the maximum processor clock frequency of 35MHz allows real-time spectral estimation to be performed for the speci fied constraints. There-programmable aspect of the FPGA is also useful; rather than designing control logic to switch between the different values of N, which uses resources and is likely to slow clock speed, a different bit-stream can be downloaded for each N. This idea can also be extended for changing to higher model order estimations where otherwise it would be difficult to paramete rise p in such a system.6. REFERENCES[1] S. M. Kay, Modern Spectral Estimation-Theory & Application. Prentice Hall, 1988.[2] M. Kassam, K. W. Johnston, and R. S. C. Cobbold, “Quantitative estimation of spectral broadening for the diagnosis of carotid arterial disease: Method and in vitro results, ”Ultrasound in Medicine and Biology, vol.11, pp. 425–433, 1985.[3] W .P. Marnane, S. J. Bellis, and P. Larsson-Edefors, “Bitserial interleaved high-speed division, ”Electronics Letters, vol. 33, pp. 1124–1125, June 1997.[4] M. M. Madeira, S. J. Bellis, M. G. Ruano, and W. P. Marnane, “Configurable processing for real-time spectral estimation, ”in Preprints of AARTC98,pp. 209–214, 1998.[5] M. G. Ruano and P. J. Fish, “Cost/benefit criterion for selection of pulsed Doppler ultrasound spectral mean frequency and bandwidth estimation, ”IEE E Transactions on Biomedical Engineering, vol. 40, no. 12, pp. 1338–1341, 1993.[6] M. G. Ruano, D. F. G. Nocetti, P. J. Fish, and P. J. Fleming, “Alternative parallel implementationsof an AR-modified covariance spectral estimator for diagnostic ultrasound blood flow studies, ” Parallel Computing,vol. 19, no. 4, pp. 463476, 1993.[7] S. J. Bellis, P. J. Fish, and W. P. Marnane, “Optimal systolic arrays for real-time implementation of the Modified Covariance spectral estimator, ”Parallel Algorithms and Applications, vol .11, no. 1-2, pp. 71–96, 1997.[8] S. J. Bellis, W. P. Marnane, and P. J. Fish, “Alternative systolic array for non-square-root Cholesky decomposition, ”IEE Proceedings: Computers and Digital Techniques, vol. 144, pp. 57–64, Mar. 1997.[9] K. Hwang, Computer Arithmetic: Principles Architecture and Design. John Wiley & Sons, 1979.中文译文单一FPGA上的参数谱估计摘要参数化谱估计模型技术可以提供针对传统的短期快速傅里叶变换提高频率分辨率的方法，克服了窗口采样，时域，输入数据造成的限制。

基于Burg 算法的AR 模型功率谱估计简介摘要：在对随机信号的分析中，功率谱估计是一类重要的参数研究，功率谱估计的方法分为经典谱法和参数模型方法。

参数模型方法是利用型号的先验知识，确定信号的模型，然后估计出模型的参数，以实现对信号的功率谱估计。

根据wold 定理,AR 模型是比较常用的模型，根据Burg 算法等多种方法可以确定其参数。

关键词：功率谱估计；AR 模型；Burg 算法随机信号的功率谱反映它的频率成分以及各成分的相对强弱, 能从频域上揭示信号的节律, 是随机信号的重要特征。

因此, 用数字信号处理手段来估计随机信号的功率谱也是统计信号处理的基本手段之一。

在信号处理的许多应用中, 常常需要进行谱估计的测量。

例如, 在雷达系统中, 为了得到目标速度的信息需要进行谱测量; 在声纳系统中, 为了寻找水面舰艇或潜艇也要对混有噪声的信号进行分析。

总之, 在许多应用领域中, 例如, 雷达、声纳、通讯声学、语言等领域, 都需要对信号的基本参数进行分析和估计, 以得到有用的信息, 其中, 谱分析就是一类最重要的参数研究。

1 功率谱估计简介一个宽平稳随机过程的功率谱是其自相关序列的傅里叶变换，因此功率谱估计就等效于自相关估计。

对于自相关各态遍历的过程，应有：)()()(121lim *k r n x k n x N N x N N n =⎭⎬⎫⎩⎨⎧++∞→∑-= 如果所有的)(n x 都是已知的，理论上功率谱估计就很简单了，只需要对其自相关序列取傅里叶变换就可以了。

但是，这种方法有两个个很大的问题：一是不是所有的信号都是平稳信号，而且有用的数据量可能只有很少的一部分；二是数据中通常都会有噪声或群其它干扰信号。

因此，谱估计就是用有限个含有噪声的观测值来估计)(jw x e P 。

谱估计的方法一般分为两类。

第一类称为经典方法或参数方法，它首先由给定的数据估计自相关序列)(k r x ，然后对估计出的)(ˆk rx 进行傅里叶变换获得功率谱估计。

已知AR(4)过程：x(n)=(n-1)(n-2)+(n-3)(n-4)+w(n)试用Burg算法估计模型参数、阶数和功率谱。

实验仿真50次，给出估计结果的方差和偏差（用估计结果均值与真值的差来估计偏差）。

用Burg算法对该AR过程估计,50次实验得平均阶数p=，取整数即为P=5，阶数方差var(p)=，阶数偏差为。

AR参数（平均）: [ ]AR参数方差:,AR参数偏差：Burg最大熵谱估计：matlab代码：u=1;while u<=50 %50次仿真实验n=10000;|w=randn(1,n);%生成零均值，方差为的白噪声w=w-mean(w);w=w/std(w);w=sqrt*w;x(1)=1+w(1);x(2)=2+w(2);x(3)=3+w(3);x(4)=4+w(4);%由ar模型生成含噪信号；for i=5:nx(i)=*x(i-1)*x(i-2)+*x(i-3)*x(i-4)+w(i);endP(1+1)=1/n*sum(x.^2);%计算预测误差功率的初始值-g(1,:)=x;f(1,:)=g(1,:);m=1;k=0;l=0;a(1,1)=1;P(1)=1;while abs(P(m+1)-P(m))/P(m)>=for j=m+1:nk=k+f(m,j)*g(m,j-1);l=l+f(m,j)^2+g(m,j-1)^2;endK(m)=-k/(1/2*l);%求反射系数a(1,1)=1;、if m>=2for j=1:m-1%计算前向预测滤波器系数a(m,j)=a(m-1,j)+K(m)*a(m-1,m-j);a(m,m)=K(m);endendP(m+1+1)=(1+abs(K(m))^2)*P(m+1);%计算预测误差功率for j=1:n%计算滤波器输出f(m+1,j)=f(m,j)+K(m)*g(m,n-1);g(m+1,j)=K(m)*f(m,j)+g(m,n-1);end（m=m+1;endM(u)=m-1;%记录每次实验的阶数u=u+1;endfor i=1:50ar=arburg(x,M(i));%估计AR系数ar_coef(i,1:length(ar))=ar;[Px,fx]=pburg(x,M(i));%估计功率谱Pxx(i,1:length(Px))=10*log10(Px');…fxx(i,1:length(fx))=fx';end%%求阶数和AR系数的方差和偏差var_p=var(M);for i=1:max(M)var_ar(i)=var(ar_coef(:,i+1));endbias_p=mean(M)-4;ar_r(max(M))=0;ar_r(1)=;ar_r(2)=;ar_r(3)=;ar_r(4)=;for i=1:max(M)`mean_ar(i)=mean(ar_coef(:,i+1));bias_ar(i)=mean(ar_coef(:,i+1))-ar_r(i);end%%50次实验功率谱平均for i=1:length(Pxx(1,:))mean_Pxx(i)=mean(Pxx(:,i));end%%作图plot(fxx(1,:),mean_Pxx);%作功率谱图title('burg最大熵估计功率谱'); xlabel('频率');ylabel('功率谱密度');》figure(2);plot(var_ar);%AR系数估计方差title('AR系数估计方差');xlabel('AR(n)');ylabel('方差');figure(3);plot(bias_ar);%AR系数估计偏差title('AR系数估计偏差');xlabel('AR(n)');ylabel('偏差');disp('阶数方差：');disp(var_p);disp('阶数偏差：');disp(bias_p);disp('AR系数:');disp(mean_ar);。

基于AR模型的功率谱估计邢务强;钮金鑫【摘要】A self-correlation method and Burg method based on AR model are usually used in the modern PSD for better simulating and computing the power spectrum. Through simulation, PSD curve is acquired by three methods. The results show that Burg method has good performance in the spectral resolution. An appropriate sampling point should be selected to determine the model order. At the same time, the principle of order selection is given.%为了更好地模拟计算功率谱,采用了现代功率谱估计中常用的基于AR模型的自相关法和Burg法,通过仿真给出了估计出的功率谱曲线,结果表明,Burg算法在谱分辨率方面有着良好的性能.选择合适的采样点数来确定模型阶数,给出阶数选择的原则.【期刊名称】《现代电子技术》【年(卷),期】2011(034)007【总页数】3页(P49-51)【关键词】AR模型;自相关法;Burg算法;阶数选择【作者】邢务强;钮金鑫【作者单位】西安邮电学院理学院,陕西西安710121;西安邮电学院通信工程学院,陕西西安710121【正文语种】中文【中图分类】TN911-340 引言所谓谱估计或功率谱估值，就是用已观测到的一定数量的样本数据估计一个平稳随机信号的功率谱。

功率谱估计分为两大类，一类是非参数化方法(又叫经典谱估计法)，另一类是参数化方法(又叫模型法估计)。

谱估计基础及仿真分析0引言 现代信号分析中,对于常见的具有各态历经的平稳随机信号,不可能用清楚的数学关系式来描述,但可以利用给定的N 个样本数据估计一个平稳随机信号的功率谱密度叫做功率谱估计(PSD)。

它是数字信号处理的重要研究内容之一。

功率谱估计可以分为经典功率谱估计(非参数估计)和现代功率谱估计(参数估计)。

功率谱估计在实际工程中有重要应用价值,如在语音信号识别、雷达杂波分析、波达方向估计、地震勘探信号处理、水声信号处理、系统辨识中非线性系统识别、物理光学中透镜干涉、流体力学的内波分析、太阳黑子活动周期研究等许多领域,发挥了重要作用。

一 相关的数学基础 1.1 概率论：1.1.1多维高斯分布高斯分布公式：22()2()x p x μσ--=(1)3σ准则为：(,)σσ-内68.26%；2σ：95.44%；3σ：99.74%。

向量的形式的公式为：11/21/211()exp[()()](2)||2T x xx x N xx p x x C x C μμπ--=---vv v v v (2) 其中[()()]Txx x x C E x x μμ=--v v v v；[]{[[]][[]]}xx ij i i j j C E x E x x E x =--(3)1.1.2 mont-carlo 仿真1ln ,(0,1)y u u λ=-∈为指数分布，y σ=1.2 随即过程平稳随即过程：多维联合概率密度和时间起点无关，狭义的平稳。

数字特征： ()[()]t E x t μ=相关函数：1212(,)[()()]xx r t t E x t x t =(4) 协方差：121122(,){[()()][()()]}xx x x C t t E x t t x t t μμ=--(5) 广义的平稳：1221(,)()()xx xx xx r t t r t t r τ=-=，12(,)()xx xx C t t C τ=(6)各态历经性：时间平均代替集平均。

数字信号处理实验报告姓名： 学号： 日期：2015.12.141. 实验任务信号为两个正弦信号加高斯白噪声，各正弦信号的信噪比均为10dB ，长度为N ，信号频率分别为1f 和2f ，初始相位021==ϕϕ，取2.01=s f f ，s f f 1取不同的数值：0.3，0.25。

s f 为采样率。

（1）分别用 Levinson 递推法和 Burg 法进行功率谱估计，并分析改变数据长度、模型阶数对谱估计结果的影响。

（2）当正弦信号相位、频率、信噪比改变后，上述谱估计的结果有何变化？并作分析说明。

2. 原理分析2.1 现代谱估计中的参数建模根据参数模型来描述随机信号的方法，我们可以知道，如果能确定信号()x n 的信号模型，根据信号观测数据求出模型参数，系统函数用()z H 表示，模型输入白噪声，其方差为2w σ，信号的功率谱用下式求出：()()22iww jw xx e H e P σ=按照这种求功率谱的思路，功率谱估计可分为三个步骤： （1）选择合适的信号模型；（2）根据()n x 有限的观测数据，或者它的有限个自相关函数的估计值，估计模型的参数；（3）计算墨香的输出功率谱。

其中以（1）、（2）两步最为关键。

按照模型的不同，谱估计的方法有许多种，它们共同的特点是对信号观测区以外的数据不假设为0，而先根据信号观测数据估计模型参数，按照求模型输出功率的方法估计信号功率谱，回避了数据观测区以外的数据假设问题。

下面分析AR 谱估计的两种方法：自相关法——列文森（Levenson ）递推法和伯格（Burg ）递推法。

这两种方法均为已知信号观测数据，估计功率谱，两者共同特点是由信号观测数据求模型系数时采用信号预测误差最小的原则。

对于长记录数据，这些方法的估计质量是相似的，但对于短记录数据，不同方法之间存在差别。

2.2 自相关法——列文森（Levenson ）递推法自相关法的出发点是选择AR 模型参数使预测误差功率最小，预测误差功率为()()()21211∑∑∑∞-∞==∞-∞=-+==n pi pi n i n x a n x Nn e Nρ假设信号()x n 的数据区在01n N ≤≤-范围，有P 个预测系数，N 个数据经过冲激响应为()0,1,pi a i p =的滤波器，输出预测误差()e n 的长度为N p +，因此应用下式计算：()()()210121011∑∑∑-+==-+=-+==P N n pi pi P N n i n x a n x Nn e Nρ()e n 的长度长于数据的长度，上式中数据()x n 的两端需补充零点，相当于对无穷长的信号加窗处理，得到长度为N 的数据。