经典谱估计算法研究与实现
基于BURG算法的谱估计研究及其MATLAB实现课件
三、经典谱估计
当f1=20, f2=40时得到的仿真结果如图3.5所示。 当f1=100, f2=1000时得到的仿真结果如图3.6 所示。
图3.5
图3.6
三、经典谱估计
当噪声方差 =1和 =10的仿真结果分别如图 3.7所示和图3.8所示。
图3.7
图3.8
三、经典谱估计
仿真结果: 间接法实现时在噪声信号很小的情况 下, 图3.5得到的谱线基本能分辨出两个频率值来, 但是也出现大量假峰。图3.6也能分辨出两个频率值, 假峰减少, 说明间接法在这种情况下得到的效果要
声信号增大到原来的100倍时就无法分辨两个频率 值, 而且通过多次仿真看出, 当噪声信号增大到10倍 以
后就不能分辨出两个频率点的峰值。
三、经典谱估计
■ 3.2 间接法及MATLAB仿真
结果
的估
■ 间接法又称自相关法, 记 为对计, 即
当M较小时, 上式的计算量不是很大, 因此, 此 方法是在FFT问世之前（即周期图被广泛应用之前） 常用的谱估计方法。
四、现代谱估计
当阶数=10时得到的仿真结果如图4.1所示。 当阶数=15时得到的仿真结果如图4.2所示。
图4.1
图4.2
四、现代谱估计
仿真结果: Levinson-Durbin算法得到的谱线波 动性小, 能很好的分辨出两个频率值, 而且没有出 现假峰现象。当增大阶数时得到的结果跟阶数小的 结果不相上下, 并没有增大频率分辨率, 反而增大 了计算次数。
仿真结果比较可得: 四种算法都能分辨出两个频 率值，可以明显看出直接法和间接法得到的仿真结果 出现了大量假峰，且谱的波动性较大，而 LevinsonDurbin算法和BURG算法得到的仿真结果很平滑，没有 出现假峰现象，能清楚的分辨出两个频率点的值，分 辨率远比经典谱估计要好，所以从仿真的结果也可以 看出现代谱估计性能比经典谱估计好。
MATLAB仿真实现经典谱估计(采用周期图法)
grid on;
window = boxcar( length( xn) ) ;%矩形窗
nfft = 512;
[Pxx f]= periodogram( xn,window,nfft,Fs) ;%直接法
subplot(312)
plot( f,10* log10( Pxx) ) ;
plot( f,10*log10( Pxx) ) ;
title('直接法经典谱估计,1024点');
xlabel('频率（Hz）');
ylabel('功率谱密度');
grid on;
六、实验总结
从上图我们可以得到这样的结论：在增加数据长度N时，就会使互不相关的点数增加，提高谱曲线的分辨力，但是加剧谱曲线 的起伏。经典功率谱估计不是一致估计，这是周期图法（直接法）的一个严重的缺点。
数字信号处理课程实验报告
实验指导教师：黄启宏
实验名称
MATLAB仿真实现经典谱估计（采用周期图法）
专业、班级
电子与通信工程
姓名
张帅
实验地点
仿古楼301
实验日期
2013.11.17
一、实验内容
采用周期图法（直接法）实现经典谱估计。
二、实验目的
（1）掌握周期图法（直接法）估计出功率谱的步骤和方法；
（2）在实验的过程中找到影响经典谱估计的因素；
%采用直接法（周期图法）估计功率谱；
clear
Fs = 1000;%采样频率
n = 0:1 /Fs: .3;%产生含有噪声的序列
xn = cos(200*pibplot(311);%输出随机信号xn；
经典功率谱估计方法实现问题的研究
1 随机信号的经典谱估计方法估计功率谱密度的平滑周期图是一种计算简单的经典方法。
它的主要特点是与任何模型参数无关，是一类非参数化方法[4]。
它的主要问题是：由于假定信号的自相关函数在数据观测区以外等于零，因此估计出来的功率谱很难与信号的真实功率谱相匹配。
在一般情况下，周期图的渐进性能无法给出实际功率谱的一个满意的近似，因而是一种低分辨率的谱估计方法。
本章主要介绍了周期图法、相关法谱估计（BT ）、巴特利特(Bartlett)平均周期图的方法和Welch 法这四种方法。
2.1 周期图法周期图法又称直接法。
它是从随机信号x(n)中截取N 长的一段，把它视为能量有限x(n)真实功率谱)(jw x e S 的估计)(jw x e S 的抽样.周期图这一概念早在1899年就提出了，但由于点数Ｎ一般比较大，该方法的计算量过大而在当时无法使用。
只是1965年FFT 出现后，此法才变成谱估计的一个常用方法。
周期图法[5]包含了下列两条假设：1．认为随机序列是广义平稳且各态遍历的，可以用其一个样本x(n)中的一段)(n x N 来估计该随机序列的功率谱。
这当然必然带来误差。
2．由于对)(n x N 采用DFT ，就默认)(n x N 在时域是周期的，以及)(k x N 在频域是周期的。
这种方法把随机序列样本x(n)看成是截得一段)(n x N 的周期延拓，这也就是周期图法这个名字的来历。
与相关法相比，相关法在求相关函数)(m R x 时将)(n x N 以外是数据全都看成零，因此相关法认为除)(n x N 外x(n)是全零序列，这种处理方法显然与周期图法不一样。
但是，当相关法被引入基于FFT 的快速相关后，相关法和周期图法开始融合。
通过比较我们发现：如果相关法中M=N ，不加延迟窗，那么就和补充（N-1）个零的周期图法一样了。
简单地可以这样说：周期图法是M=N 时相关法的特例。
因此相关法和周期图法可结合使用。
2.2 相关法谱估计（BT ）法这种方法以相关函数为媒介来计算功率谱，所以又叫间接法。
chapter_3_经典谱估计
式中P上的 表示估计值。下标p表示 式中 上的 “∧” 号，表示估计值。下标 表示 周期图。很显然， 与真实谱P(k)有 周期图 。 很显然 ， 谱估计p(k) 与真实谱 有 一定的偏差，样本越短偏差可能越大。 一定的偏差，样本越短偏差可能越大。但样本 大小又会受到采样时间和运算时间的限制。 大小又会受到采样时间和运算时间的限制。
五、平均平滑周期图
平均平滑周期图法就是先求信号的分段平滑（加窗） 平均平滑周期图法就是先求信号的分段平滑（加窗） 谱, 然后再求平均谱。 然后再求平均谱。
5 非参数功率谱密度估计方法
• • • • 周期图法 Bartlett法（平均多个周期图, 采用不同数据块） Welch 法 (平均多个周期图, 采用重叠的数据块) Blackman-Tukey 法 (周期图平滑)
四、平滑周期图
平滑周期图分析中的所谓平滑的方法是对样本加各种窗函数使 之平滑化。加窗是用等长的窗函数序列与样本数据序列相乘。 设样本数据序列为x(n)，窗函数为W(n)，则加窗后的序列 xw(n)为: xw(n) = x(n) W(n) (10-8a) 加窗也可在频域完成。根据傅立叶变换定理，两信号乘积的 傅立叶变换等于两信号分别进行傅立叶变换后的结果的卷积 XW(K)可得： XW(K) = X(K)*W(K) (10-8b) 式中K为频域自变量。 下面列举几种生物医学信号处理中常用的窗函数。
高阶累积量(higher-order cumulant)提供 高阶累积量 提供 了非高斯、 了非高斯、非线性信号的高阶相关性的量 对于高斯信号， 度。对于高斯信号，其三阶以上的累积量 皆为0。这就提供了一种可能性： 皆为 。这就提供了一种可能性：采用高阶 统计量分析信号，可以大大提高信噪比。 统计量分析信号，可以大大提高信噪比。 高阶累积量的傅利叶变换叫高阶谱 （higher-order spectrum），如三阶累积 ），如三阶累积 ）， 量的傅利叶变换叫二阶累积量谱或双谱 ）。高阶谱分析可用以检测 （bispectrum）。高阶谱分析可用以检测 ）。 被强高斯噪声淹没的有用信息。 被强高斯噪声淹没的有用信息。
经典功率谱估计
雷达和声呐系统
目标检测
在雷达和声呐系统中，经典功率谱估计常被用于目标检测。通过对接收到的信号进行功率 谱分析，可以判断是否存在目标以及目标的位置和速度等信息。
距离和速度测量
在雷达和声呐系统中，经典功率谱估计还可以用于距离和速度测量。通过对接收到的信号 进行功率谱分析，可以估计出目标与系统之间的距离和相对速度。
信号分类
在雷达和声呐系统中，经典功率谱估计还可以用于信号分类。通过对接收到的信号进行功 率谱分析，可以判断目标的类型，例如区分飞机、船舶或车辆等不同类型目标。
05 经典功率谱估计的改进方 法
基于小波变换的功率谱估计
1
小波变换能够将信号分解成不同频率和时间尺度 的分量，从而更好地揭示信号的内在结构和特征。
然而，这些方法通常需要较长 的数据长度和较为复杂的计算 过程，对于短数据和实时处理 的应用场景具有一定的局限性 。
研究展望
01
随着信号处理技术的发展，经典功率谱估计方法仍有进一步优化的空 间。
02
针对短数据和实时处理的应用场景，研究更为快速、准确的功率谱估 计方法具有重要的实际意义。
03
结合机器学习和人工智能技术，探索基于数据驱动的功率谱估计方法 是一个值得关注的方向。
优点
能够提供较高的频率分辨率和较低的估计误差。
原理
格莱姆-梅尔谱估计利用了信号的模型参数，通过 构造一个模型函数来描述信号的频率响应特性， 并求解该函数的极值问题得到信号的功率谱。
缺点
需要预先设定模型函数的形式和参数，且计算复 杂度较高。
03 经典功率谱估计的优缺点
优点
01
02
03
算法成熟
经典功率谱估计方法经过 多年的研究和发展，已经 相当成熟，具有较高的稳 定性和可靠性。
基于BURG算法的谱估计研究和MATLAB实现毕业论文
毕业设计（论文）题目：基于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上的参数谱估计摘要参数化谱估计模型技术可以提供针对传统的短期快速傅里叶变换提高频率分辨率的方法，克服了窗口采样，时域，输入数据造成的限制。
基于多数据窗的谱估计算法研究及仿真
基于多数据窗的谱估计算法研究及仿真1.1研究思路本课题主要对经典谱估计算法及多数据窗进行研究与分析，在此基础上，将多数据窗与谱：估计算法进行结合，进而利用Matlab进行仿真程序设计，并利用实际汽轮机振动信号的谱分析，主要研究思路如下：1、本课题主要是对经典谱估计算法及多数据窗进行研究，首先对经典谱估计各种算法介绍与分析，算法以傅里叶变换为基础，对平稳随机信号的功率谱分析具有良好效果，目前常用的有周期法、相关图法、改进的周期图法等。
2、其次对多窗功率谱估计算法分析介绍，该算法是将数据工作区外的未知数据假设为零,相当于数据加窗。
有如下研究方法，多窗谱估计算法：多窗谱估计算法也称为MTM算法，是由Thomon提出，建立在经典谱估计算法的基础之上．其算法主要思路是对同一数据序列加多个正交的数据窗分别求直接谱（特征谱），然后加权平均得到谱估计。
3、掌握多窗谱估计和多正弦窗谱估计的算法原理后，用Matlab对其进行编程，并编写图形用户界面程序，形成一套谱估计算法仿真程序。
4、基于实测的汽轮机振动信号对多窗谱估计和多正弦窗谱估计的算法性能进行分析。
5.最后结合实际信号，对算法效果及性能分析，提出其改进算法并加以实现。
1.2相关技术选择1.2.1编程工具MATLAB语言是当今国际上科学界(尤其是自动控制领域)最具影响力、也是最有活力的软件。
它起源于矩阵运算，并已经发展成一种高度集成的计算机语言。
它提供了强大的科学运算、灵活的程序设计流程、高质量的图形可视化与界面设计、便捷的与其他程序和语言接口的功能。
MATLAB语言在各国高校与研究单位起着重大的作用。
MATLAB的含义是矩阵实验室（MATRI某LABORATORY），主要用于方便矩阵的存取，其基本元素是无须定义维数的矩阵。
MATLAB自问世以来,就是以数值计算称雄。
MATLAB进行数值计算的基本单位是复数数组（或称阵列），这使得MATLAB高度“向量化”。
经过十几年的完善和扩充，现已发展成为线性代数课程的标准工具。
经典谱估计(自相关法)
经典谱估计(自相关法)
经典谱估计是一种常用的信号处理方法，其中自相关法是其中一种常见的实现方式。
经典谱估计的主要目的是通过对信号的自相关函数进行分析来估计信号的频谱特性。
自相关函数描述了信号与自身在不同时间点的相关性，通过对自相关函数进行合适的处理，可以得到信号的频谱信息。
自相关法的基本原理是利用信号的自相关函数来估计信号的频谱特性。
自相关函数描述了信号在不同时间点上的相关性，它可以通过计算信号与其自身在不同时间延迟下的乘积来得到。
在实际应用中，可以使用不同的自相关函数估计方法，如周期图谱法、傅里叶变换法等。
在进行自相关法时，需要考虑一些关键因素。
首先是选择合适的信号长度和时间窗口大小，这会影响到自相关函数的准确性和分辨率。
其次是对信号进行预处理，如去除噪声、进行平滑处理等，以提高自相关函数的稳定性和可靠性。
另外，还需要考虑自相关函数的计算方法和参数选择，以确保得到准确的频谱估计结果。
自相关法在实际应用中有着广泛的应用，特别是在信号处理、
通信系统和频谱分析等领域。
它可以用于估计信号的频谱特性，如频率成分、功率谱密度等，对于信号的特征提取和分析具有重要意义。
同时，自相关法也可以用于信号的调制识别、信道估计和系统建模等方面，为工程实践提供了有力的工具和方法。
总的来说，经典谱估计中的自相关法是一种重要的信号处理方法，通过对信号的自相关函数进行分析来估计信号的频谱特性。
在实际应用中，需要综合考虑信号处理的各个环节，合理选择方法和参数，以获得准确可靠的频谱估计结果。
经典谱估计方法
经典谱估计方法嘿，咱今儿就来唠唠经典谱估计方法。
你说这谱估计啊，就像是给一个神秘的信号画像，让咱能看清它的真面目。
想象一下，信号就像一个调皮的小精灵，在我们面前蹦来蹦去，一会儿藏起来，一会儿又冒出来。
而经典谱估计方法呢，就是我们用来抓住这个小精灵的工具啦。
常见的经典谱估计方法有好几种呢，比如周期图法。
这就好像是给小精灵拍了一张快照，能大致看出它的模样。
不过呢，这张快照有时候可不太准哦，会有一些偏差和误差呢。
还有自相关法，这就像是跟小精灵玩一个追踪游戏，通过它留下的痕迹来推测它的样子。
它能让我们对小精灵的了解更深入一些，但也不是十全十美的呀。
你可能会问了，那这些方法都有啥用呢？哎呀，用处可大了去啦！比如说在通信领域，我们得搞清楚信号的特点，才能让信息准确无误地传递呀。
就好像你给朋友发消息，总不能乱码一堆吧。
在音频处理里，经典谱估计方法能帮我们把声音变得更好听，更清晰。
就像给声音化了个妆，让它更漂亮。
在科学研究中，那更是少不了它。
研究各种现象，分析数据，都得靠它来帮忙呢。
可是啊，经典谱估计方法也不是完美无缺的。
它就像一个有点小脾气的朋友，有时候会闹点小情绪，给咱出点难题。
比如说分辨率不高啦，容易受到噪声干扰啦。
那怎么办呢？咱就得想办法哄好它呀，或者找些其他的办法来弥补它的不足。
就像你有个朋友脾气不太好，你得有耐心，还得想办法和他好好相处。
总之呢，经典谱估计方法是我们探索信号世界的重要工具，但我们也得清楚它的优缺点，灵活运用，才能让它发挥最大的作用呀。
别小看了这些方法，它们可是能帮我们解开很多信号的秘密呢！这就是经典谱估计方法，有趣吧？嘿嘿！。
通信工程论文题目推荐
一、通信：1．无线通信信道均衡技术的研究自适应均衡器的研究与仿真设计自适应均衡器及其发展趋势2．信道估计技术研究3．调制解调通信系统中调制技术的研究及其matlab仿真4．信源编码GSM移动通信系统中的语音编码技术研究语音编码及其在移动通信系统中的应用5．信道编码网络编码网技术在协作通信系统中的应用基于网络传输的数字喷泉码研究纠错码及其在移动通信系统中的应用卷积编码在移动通信中的应用及性能分析RS编码在移动通信中的应用及性能分析6．扩频通信系统性能分析与仿真7．LTELTE物理层信道编码技术分析8.物联网物联网核心技术-无线射频识别技术研究无线传感器网络定位技术研究基于无线传感器网络的安全路由协议研究无线传感器网络传输协议研究物联网能量均衡路由算法研究物联网安全与信任机制研究基于信息隐藏的传感器网络安全研究物联网节点任务划分、特征提取与建模方法研究应用无线传输在智能交通上的应用开发二．信号处理1．经典谱估计算法研究与实现2．现代谱估计算法研究与实现（1）基于AR模型的噪声源识别方法（2）空间谱估计中经典算法的研究(3)基于MUSIC的空间谱估计算法的研究与仿真3．自适应滤波器的研究与仿真设计4．多媒体数据压缩方法研究音频、视频5．图像信号处理：图像增强、图像复原基于维纳滤波的图像增强算法研究数字图像增强技术的研究数字图像消噪技术研究（小波变换在图像去噪中的应用研究基于小波变换的图像去噪技术基于Matlab的图像去噪算法的研究与实现）基于DCT变换的图像压缩研究矢量量化技术在图像压缩中的应用图像数字水印技术研究及matlab实现6．压缩标准视频压缩编码标准H.264中的运动估计研究与实现G.729A语音编解码算法分析与实现基于VC6.0的MPEG-4视频编解码分析与实现。
数字信号处理经典谱估计
Nm
r(m)
N
全部样本
所以：
m bi a[rˆ (m)] r(m)
N
含义
(1) m 固定，N , bia[rˆ(m)] 0
渐近无偏估计
(2) N给定，m N, rˆ(m) 接近 r(m) 对固定旳N，此结论 给出了m旳选用原则
(3) E{rˆ(m)} w(m)r(m)
(3) E{rˆ(m)} w(m)r(m)
2.自有关（Blackman-Tukey BT法）法:
Step1
rˆx (m)
1 N
N 1 m n0
x(n)x(n
m)
M
Step2 PˆBT ()
rˆx (m)e jm , M N 1
mM
因为先要估计自有关函数，所以 又称间接法。与此相相应，周期 图法又称直接法。
3.直接法和间接法旳关系：
X N (k) 2
N点离散谱
怎样和 PˆB2TN (k ) 相等？
PˆB2TN (k)
( N 1)
j 2 mk
rˆx (m)e 2N
m( N 1)
k 0,1, 2N 1
PˆBNT
(k
)
2
N
Re[
1
rˆ0
(m)e
j
2 N
mk
]
rˆ(0)
m0
k 0,1, N 1
N点离散谱
（二） M N 1 相当于只用了部分自有关函数
N
经典功率谱估计不是一致估计
解释：
1
2N
P()D0 (
)D0(
2
)d 0
D0 ()
rˆM (m) rˆ(m)v(m), v(m) : M ~ M
功率谱估计的经典方法
∞
=
Ryy (m) =
p =−∞
∑R
k = −∞ ∞
∑ h( k ) R
xx
∞
xx
( m − k ) = Rxx (m) ∗ h( m)
(m − p) Rhh ( p) = Rxx (m) ∗ Rhh (m)
或
= Rxx (m) ∗ h(m) ∗ h(−m) = Rxy (m) ∗ h(−m)
S yy (e jω ) = S xy (e jω ) H (e− jω )
jω
jω
jω
2
离散随机信号通过线性非移变系统
（4）输入随机过程与输出随机过程的互相关序列Rxy(m)
∞ Rxy ( m) = E [x ( n) y ( n + m) ] = E x ( n) ∑ h( k ) x ( n + m − k ) k = −∞
=
k = −∞
∑ h(k ) E[x(n) x(n + m − k )]
ˆ B =α − E [ α ]
无偏估计， 无偏估计 有偏估计，当观测数据为无穷时B = 0，则称其为渐 渐 B = 0时无偏估计 B ≠ 0 有偏估计 进无偏估计。无偏估计和渐进无偏估计又称为是好估计 进无偏估计 好估计。 好估计
均值 均方值
E[xn ] = mxn = ∫ xpxn ( x, n)dx
∞ −∞
E x = ∫ x 2 pxn ( x, n)dx
2 n −∞
[ ]
2
∞
方差
E xn − mxn
[(
) ]= σ
2 xn
=∫
∞ −∞
(x − m )
xn
随机信号的谱估计方法及实现+孟波+20091342083
淮北师范大学2013届学士学位论文随机信号的谱估计方法及实现学院、专业物理与电子信息学院电子信息工程研究方向数字信号处理学生姓名孟波学号20091342083指导教师姓名崔少华指导教师职称讲师2013年4月26日随机信号的谱估计方法及实现孟波淮北师范大学物理与电子信息学院235000摘要对一个确定性信号来进行傅里叶变换是频率域分析研究中的理论基础，但是对于一个随机信号来说，它是不存在傅里叶变换的，因而我们转而去研究它的功率谱。
功率谱估计是数字信号处理中的一个重要组成部分，其理论研究已经非常成熟，但是由于计算量的庞大传统的仿真实现往往令人不愿采用那种方法。
采用功能强大的MATLAB软件对常用的几种线性法功率谱估计进行仿真，通过对于仿真结果进行研究讨论总结出这几种功率谱估计法的各自特点，并对其比较分析，通过对其特点仔细研究，从而在实际工作中做出合理的选择。
本文简要介绍了MATLAB仿真软件的一些基本的功能、使用的特点和特有的优势等等，以及详细介绍了周期图法、平均周期图法、窗函数法等几种经典的功率谱估计方法的原理，并对各种方法进行仿真，通过观察到的结果来分析其优缺点，对这几种方法进行质量评价，并提出了改进方向。
关键词自相关函数;功率谱;MATLAB;随机信号Random signal spectrum estimation method and itsimplementationMeng BoSchool ofPhysics and Electronic Information, Huai Bei Normal University, Anhui Huaibei, 235000 Abstract The Fourier transform of the deterministic signal in the frequency domain analysis is a basical, but for random signals, the Fourier transform does not exist. So we turn to study its power spectrum. Power spectrum estimation is an important part of the digital signal processing.Its theories is very mature now. Because of the hard of the traditional simulation, a lot ofpeople refuse to use it with powerful MATLAB to replace.The results of the simulation are discussed to summarize their respective characteristics of these kinds of power spectrum estimation methods and the comparative analysis. Through the carefully study of its characteristics,It results in the practical work to make the rational choice. This paper briefly introduces the MATLAB simulation software, the use of some of the basic functions of the characteristics , the unique advantage and so on, and the details of the power of several periodic graph method, average periodogram method, window function method and the classicalspectrum estimation method, and the method of simulation, the observed result analysis its advantages and disadvantages, to evaluate the quality of these methods and put forward the improvement direction.Keywords autocorrelation function; power spectrum; MATLAB; random signal目录1 绪论 (5)1.1常用的功率谱研究方法及其发展 (5)1.2功率谱估计方法的产生 (6)1.3功率谱估计应用方向 (6)1.4本课题的主要研究内容 (7)1.5MATLAB应用方向 (7)1.6MATLAB特点及其优势 (8)2 古典法对随机信号的谱估计 (9)2.1周期图法 (9)2.2相关法谱估计（BT） (9)2.3W ELCH法 (11)3 分析比较各种估计方法 (12)3.1直接法 (12)3.2间接法 (13)3.3W ELCH法： (14)3.4对比分析各种估计方法 (16)结论 (17)参考文献 (18)附录 (19)致谢 (21)1 绪论在对信号和系统进行分析研究和处理的时候我们常用的方法主要有两种，一是对信号在时域进行分析处理，二是对信号在频率进行处理研究。
经典功率谱估计方法研究
中文摘要介绍了各种经典功率谱估计方法,不仅从理论上对各种方法的谱估计质量进行了分析比较,而且通过Matlab 实验仿真验证了理论分析的正确性。
着重对使用比较广泛的Welch 法进行了深入的研究,给出了窗函数选择的一般要求,通过仿真分析了不同的窗函数对Welch 法谱估计质量的影响,比较了他们的优缺点。
最后分析了采样点数较少即短数据对Welch 法谱估计质量的影响。
关键词:经典谱估计;估计质量;Welch 法;窗函数;短数据AbstractVarious classical Power Spect rum Density ( PSD) estimation methods are int roduced ,estimation quality of eachmethod is analyzed and compared in both theory and simulation using the sof tware Matlab. Then further study is made inWelch method which is used most widely. General selecting criterion of window function is presented and estimation quality ofWelch method using different window function is compared. Finally ,the impact of fewer data on estimation quality of Welchmethod is analyzed.Keywords:classical PSD estimation ;estimation quality ;Welch method ;window function ;fewer data第1章绪论 (4)1.1 引言 (4)1.2 选择背景与意义 (4)1.3 经典谱估计发展和应用 (4)第2章经典功率谱估计 (5)2.1 引言 (5)2.2 自相关函数法的估计 (10)2.3 周期图作为功率谱的估计 (13)2.4 经典功率谱估计方法的改进 (19)2.4.1 巴特利特(Bartlett)平均周期图的方法 (19)2.4.2 Welch法 (23)第3章 MATLAB仿真 (24)3.1 仿真结果 (24)3.2 仿真结果分析 (24)3.3 不同窗函数的Welch 谱估计 (25)3.4 短数据的Welch 谱估计 (25)3.5 结论 (26)第4章周期图法和Welch法的比较 (27)4.1 周期图法和Welch法 (27)4.1.1周期图法 (27)4.1.2 Welch法 (27)4.2算法流程图、MATLAB程序及谱估计的分析 (27)4.2.1 算法流程 (28)4.2.2 程序 (28)第5章总结 (30)第1章绪论1.1 引言信号的频谱分析是研究信号特性的重要手段之一,对于确定性信号,可以Fourier 变换来考察其频谱性质,而对于广义平稳随机信号,由于它一般既不是周期的,又不满足平方可积,严格来说不能进行Fourier 变换,通常是求其功率谱来进行频谱分析。
经典谱估计算法性能比较
经典谱估计算法性能比较经典谱估计算法是信号处理领域中常用的一类算法，用于从观测到的信号样本中估计信号的频率、振幅、相位等相关参数。
常见的谱估计算法有传统谱估计法、非参数谱估计法和最小二乘谱估计法等。
本文将从算法原理、性能指标和实际应用等方面，对这些经典谱估计算法进行比较和分析。
一、算法原理传统谱估计法是最简单、常用的一类谱估计算法，其基本思想是通过对信号进行线性变换，将频谱估计问题转化为参数估计问题。
常见的传统谱估计算法有周期图法、自相关函数法、特定窗函数法等。
非参数谱估计法则是基于信号样本的统计特性，通过对信号样本进行直接分析来估计信号的频谱。
最常用的非参数谱估计算法有周期图法、Welch法、多普勒谱估计法等。
这类算法通常具有计算量大、辨识能力强的特点。
最小二乘谱估计法是利用线性最小二乘法原理，通过优化目标函数来估计信号的谱。
最小二乘谱估计法的核心是通过最小化残差平方和来获得最佳估计值。
常见的最小二乘谱估计算法有波前源谱估计法、Capon谱估计法等。
二、性能指标1.分辨率：性能指标之一是分辨率，即算法在估计信号频谱时，能否分辨出不同频率成分的能力。
分辨率越高，代表信号的频谱估计结果越精确。
2.偏差：性能指标之二是偏差，即估计结果与真实值之间的差异。
偏差越小，代表算法的估计结果越接近真实值。
3.方差：性能指标之三是方差，即估计结果的波动程度。
方差越小，代表算法的稳定性较好，估计结果相对较稳定。
4.频谱动态范围：性能指标之四是频谱动态范围，即算法在估计信号频谱时，能够估计到的最小和最大频率的能力。
频谱动态范围越宽，代表算法的适用范围越广。
1.分辨率比较：传统谱估计法的分辨率相对较低，非参数谱估计法的分辨率较高，而最小二乘谱估计法的分辨率介于传统谱估计法和非参数谱估计法之间。
2.偏差比较：传统谱估计法的偏差较大，非参数谱估计法的偏差较小，而最小二乘谱估计法的偏差相对较小。
3.方差比较：传统谱估计法的方差较大，非参数谱估计法的方差较小，而最小二乘谱估计法的方差相对较小。
5 谱估计(概述和经典法)
谱 估 计主要内容•引言•经典谱估计•现代谱估计1 引 言✶概述✶估计质量的评价✶功率谱估计的应用✶研究现状•估计质量的评价的偏差(Bias)为零 。
所谓偏差(用B 表示)定义为 无偏估计θ：某个随机变量的真值：它的估计值 ˆθˆθˆˆ[]()B Bias E θθθ∆∆-☠估计1和估计2都属于无偏估计；☠估计2较之估计1方差小；•估计质量的评价均方误差θ：某个随机变量的真值：它的估计值 ˆθ不难证明：22ˆˆ()()MS E e E θθθ⎡⎤⎡⎤==-⎣⎦⎣⎦222ˆE e B θσ⎡⎤=+⎣⎦当N 趋向于无穷大时，谱估计趋向于真实的谱密度。
•估计质量的评价一致估计：ˆ 0ˆ ar 0N Bias N V θθ⎫⎡⎤→∞→⎣⎦⎪⎬⎡⎤→∞→⎪⎣⎦⎭正确的估计应该满足一致估计的条件，此为正确估计的必要条件 反之，若估计方法不满足一致估计的条件，则它一定是不正确的1 引 言•功率谱估计的应用☞在信号处理的许多场所，要求预先知道信号的功率谱密度(或自相关函数)。
☞常常利用功率谱估计来得到线性系统的参数估计。
☞从宽带噪声中检测窄带信号。
•功率谱估计的应用谱估计的分辨率可以粗略地定义为能够分辨出的二个分立的谱分量间的最小频率间隙(距)。
例如：有一个随机信号，它包括二个频率相差1Hz振幅相等的正弦波以及加性白噪声(白色噪声的方差是正弦波功率的10％)。
用三种不同的谱估计方法检测这二个正弦分量的效果。
(a) 经典BT PSD法(b) 最大熵谱估计法(c) Pisavcnko 谐波分解法•研究现状功率谱估计的方法：教材P489 图10.7.1•研究现状☞经典谱估计：固有缺陷：原因：“加窗效应”频率分辨率低原因：加窗截取，认为窗以外的数据为零。
频谱能量向旁瓣泄漏原因：加窗截取，频域产生旁瓣和主瓣宽度不是无限窄的现象。
周期图的缺陷：非一致估计当数据量增至无限多时，周期图的方差并不趋近于零，而是趋近于常数。
矩形序列其傅立叶变换为幅度谱各种窗函数的频谱2 经典谱估计•自相关函数的估计•周期图作为功率谱的估计•平滑后的周期图作为PSD的估计2.2 周期图法进行谱估计求出信号的自相关函数，再求出信号的功率谱密度。
数字信号处理——谱估计实验
实验三 随机信号的功率谱估计方法报告人： #### 报告时间：2013.1.2一、实验目的1．利用自相关函数法和周期图法实现对随机信号的功率谱估计。
2．观察数据长度、自相关序列长度、信噪比、窗函数、平均次数等对谱估计的分辨率、稳定性、主瓣宽度和旁瓣效应的影响。
3．学习使用FFT 提高谱估计的运算速度。
4．体会非参数化功率谱估计方法的优缺点。
二、实验原理与方法假设信号x (n )为平稳随机过程，其自相关序列定义为*()[()()]m E x n x n m φ∆=+ (3-1)其中E 表示取数学期望，*表示共轭运算。
根据定义，x (n )的功率谱密度()P ω与自相关序列()m φ存在下面关系： ()()j m m P m e ωωφ∞-=-∞=∑ (3-2)1()()2j m m P e d πωπφωωπ-=⎰ (3-3)然而，实际中我们很难得到准确的自相关序列()m φ，只能通过随机信号的一段样本序列来估计信号的自相关序列，进而得到信号的功率谱估计。
目前，常用的线性谱估计方法有两种，即自相关函数法和周期图法。
1． 自相关函数法假设我们已知随机信号x(n)的M 长的自相关序列{()m φ}，利用自相关函数法可以得到x(n)的功率谱估计：*11()()()L m i m x i x i m L m φ-Λ==+-∑ 11ˆˆ()()M j m m M P m e ωωφ--=-+=∑ (3-4) 利用窗函数，上式又可表达为ˆˆ()()()R j m M m PW m m e ωωφ∞-=-∞=∑ (3-5) 其中，()R M W m 为矩形窗函数，定义为1()0RM m M W m m M <⎧=⎨≥⎩ (3-6)因此，ˆ()P ω实际上是：真正功率谱()P ω与窗函数()R MW m 傅立叶变换的卷积。
矩形窗函数不仅降低了谱估计的分辨率，而且使谱估计产生了旁瓣。
为了降低旁瓣影响，可以采用具有较小旁瓣的窗函数，如Hamming 窗，它定义为0.540.46cos ()0H M m M m W m M m Mπ⎧<+⎪=⎨≥⎪⎩ (3-7) 这种窗函数可以有效的抑制旁瓣，但是这种方法使主瓣宽度增大，从而降低了谱估计的分辨率，这种主瓣大小和旁瓣干扰之间的矛盾在线性谱估计方法中是无法解决的。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
毕业设计（论文）论文题目：经典谱估计算法研究与实现教学中心：电子科技大学网络教育学院苏州学习中心指导老师：职称：学生姓名：学号：专业：通信工程毕业设计（论文）任务书题目：经典谱估计算法研究与实现任务与要求：探讨我国经典谱估计算法运用和影响下所面临的机遇与挑战以及经典谱估计算法研究与实现,结合所知识,理论联系实际,写出毕业论文。
时间：2014年 1 月25日至 2014年 4 月 14日共 12 周教学中心：电子科技大学网络教育学院苏州学习中心学生姓名：学号：专业：通信工程指导单位或教研室：电子科技大学网络教育学院苏州学习中心指导教师：职称：毕业设计(论文)进度计划表电子科技大学毕业设计(论文)中期检查记录表注:此表同学生毕业设计(论文)一起存档1 随机信号的经典谱估计方法估计功率谱密度的平滑周期图是一种计算简单的经典方法。
它的主要特点是与任何模型参数无关，是一类非参数化方法[4]。
它的主要问题是：由于假定信号的自相关函数在数据观测区以外等于零，因此估计出来的功率谱很难与信号的真实功率谱相匹配。
在一般情况下，周期图的渐进性能无法给出实际功率谱的一个满意的近似，因而是一种低分辨率的谱估计方法。
本章主要介绍了周期图法、相关法谱估计（BT ）、巴特利特(Bartlett)平均周期图的方法和Welch 法这四种方法。
2.1 周期图法周期图法又称直接法。
它是从随机信号x(n)中截取N 长的一段，把它视为能量有限x(n)真实功率谱)(jw x e S 的估计)(jw x e S 的抽样.周期图这一概念早在1899年就提出了，但由于点数Ｎ一般比较大，该方法的计算量过大而在当时无法使用。
只是1965年FFT 出现后，此法才变成谱估计的一个常用方法。
周期图法[5]包含了下列两条假设：1．认为随机序列是广义平稳且各态遍历的，可以用其一个样本x(n)中的一段)(n x N 来估计该随机序列的功率谱。
这当然必然带来误差。
2．由于对)(n x N 采用DFT ，就默认)(n x N 在时域是周期的，以及)(k x N 在频域是周期的。
这种方法把随机序列样本x(n)看成是截得一段)(n x N 的周期延拓，这也就是周期图法这个名字的来历。
与相关法相比，相关法在求相关函数)(m R x 时将)(n x N 以外是数据全都看成零，因此相关法认为除)(n x N 外x(n)是全零序列，这种处理方法显然与周期图法不一样。
但是，当相关法被引入基于FFT 的快速相关后，相关法和周期图法开始融合。
通过比较我们发现：如果相关法中M=N ，不加延迟窗，那么就和补充（N-1）个零的周期图法一样了。
简单地可以这样说：周期图法是M=N 时相关法的特例。
因此相关法和周期图法可结合使用。
2.2 相关法谱估计（BT ）法这种方法以相关函数为媒介来计算功率谱，所以又叫间接法。
它是1958年由Blackman 和Tukey 提出。
这种方法的具体步骤是：第一步：从无限长随机序列x(n)中截取长度N 的有限长序列列)(n x N 第二步：由N 长序列)(n x N 求（2M-1）点的自相关函数)(m R x ∧序列。
即)()(1)(10m n x n xNm R N n N Nx +=∑-=∧(2-1)这里，m=-(M-1)…,-1,0,1…,M-1，M N ，)(m R x 是双边序列，但是由自相关函数的偶对称性式，只要求出m=0,。
,M-1的傅里叶变换，另一半也就知道了。
第三步：由相关函数的傅式变换求功率谱。
即jwm M M m Xjwx e m Re S ----=∧∧∑=)()(1)1((2-2)以上过程中经历了两次截断，一次是将x(n)截成N 长，称为加数据窗，一次是将x(n)截成（2M-1）长，称为加延迟窗。
因此所得的功率谱仅是近似值，也叫谱估计，式中的)(jw x e S 代表估值。
一般取M<<N ，因为只有当M 较小时，序列傅式变换的点数才较小，功率谱的计算量才不至于大到难以实现，而且谱估计质量也较好。
因此，在FFT 问世之前，相关法是最常用的谱估计方法。
当FFT 问世后，情况有所变化。
因为截断后的)(n x N 可视作能量信号，由相关卷积定理可得)(m R x ∧)]()([1m x m x NN N -*=(2-3)这就将相关化为线性卷积，而线性卷积又可以用快速卷积来实现。
我们可对上式两边取（2N-1）点DFT ，则有2121212)(1)()([1)(K X NK x k x N m R N N N x ---∧=*=(2-4)于是将时域卷积变为频域乘积，用快速相关求)(m R x ∧的完整方案如下： 1.对N 长)(n x N 的补充（N-1）个零，成为（2N-1）长的。
2.求（2N-1）点的FFT ，得∑-=----=220121212)()(N N mk N N N W n xK X 。
3.求212)(1K X NN -。
由DFT 性质，)(12n x N -是纯实的，)(12k x N -满足共轭偶对称，而212)(1K X NN -一定是实偶的，且以（2N-1）为周期。
4. 求（2N-1）点的IFFT ：mk N N N k N x W K X N N m R -----=-∧∑-=121)1(212)(1121)((2-5)这里212)(1K X NN -是实偶的，m=-(N-1)...0...N-1。
本来IFFT 求和范围是0至2N-2，由于212)(1K X NN -的实偶性与周期性，求和范围改为-（N-1）至（N-1）不影响计算结果。
同理可将m 的范围改为-（N-1）至（N-1）。
上述的快速相关中，补充零的目的是为了能用圆周卷积代替线性卷积，以便进一步采用快速卷积算法。
快速相关输出是-（N-1）至（N-1）的2N-1点，加)(m W M 窗后截取的是-（M-1）至（M-1）的频段，最后作（2M-1）点FFT ，得)(k S x ∧。
我们注意到：如果数据点数与自相关序列点数相同即M=N ，则（2N-1）点的IFFT 后紧跟一个（2N-1）点的FFT ，利用)(m R x ∧的对称性，FT 运算框的计算式变为=)(K S X ∑---=∧1)1()(N N m Xm R1212---N mkN X W(2-6)由于N=M 并假设窗形状是矩形的，第二次()m W M 的截断就不需要了。
比较式（2-5）和式（2-6）， ](m)R FFT[k S x x ∧=）（,])(1[)( 212K X NIFFT m R N x -∧=正反傅氏变换可以抵消，直接得）（k S x =212)(1K X NN - (2-7)为了实行基2FFT ，也可将（2N-1）点换成2N 点，这样做不影响结果的正确性。
2.3 巴特利特(Bartlett)平均周期图法首先让我们来看一下为什么周期图经过某种平均(或平滑)后会使它的方差当∞→N 时趋于零，达到一致估计的目的。
如果L x x x , , ,21 是不相关的随机变量，每一个具有期望值μ，方差2σ，则可以证明它们的数学平均L x x x x l /)...(21+++=的期望值等于μ，数学平均的方差等于L /2σ，即：[][]μμ=⋅=+++=L Lx x x E L x E L 1121 [][][]222)(][))((x E x E x E x E x Var -=-=[]22212)(1μ-+++=L x x x E L()[]2112222121μ-⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧++++=∑∑=≠=L j L j i i ji L x x E x x x E L [][][]∑∑∑∑=≠==≠=⋅=L j Lji i ijL j ji i ji x E x E x x E 1111222)1(μμμμL L L L -=-=所以[][][]⎥⎦⎤⎢⎣⎡-=-⎥⎦⎤⎢⎣⎡-+=∑∑==2122222122211μμμμL x E L L L x E L x Var L i i L i i[][][][][][]{}22222221212])[(])[(])[(1L L x E x E x E x E x E x E L-++-+-=L L L2221σσ==(2-8)由(2-8)可见，L 个平均的方差比每个随机变量的单独方差小L 倍。
当[]0 →∞→x Var L ，可达到一致谱估计的目的。
因而降低估计量的方差的一种有效方法是将若干个独立估计值进行平均。
把这种方法应用于谱估计通常归功于Bartlett 。
Bartlett 平均周期图的方法是将序列)10( )(-≤≤N n n x 分段求周期图再平均。
设将)(n x 分成L 段，每段有M 个样本，因而LM N =，第i 段样本序列可写成L i M n M iM n x n x i ≤≤-≤≤-+=1 ,10 )()(第i 段的周期图为210)(1)(∑-=-=M n n j ji M e n xMI ωω如果)( ,m M m xx φ>很小，则可假定各段的周期图)(ωiMI 是互相独立的。
对功率谱密度的概念的讨论，谱估计可定义为L 段周期图的平均，即∑==L i iM xx I L P 1)(1)(ˆωω(2-9)于是它的期望值为[][][]∑===L i i M i M xx I E I E L P E 1)()(1)(ˆωωω []⎰--=ππθωθθπωd e W PP E j B xxxx)()(21)(ˆ)(()[()[⎰-⎥⎦⎤⎢⎣⎡--=ππθθωθωθπd M P Mxx 22/sin 2sin )(21 (2-10) 这里L N M /=，因此Bartlett 估计的期望值是真实谱)(ωxx P 与三角窗函数的卷积。
由于三角窗函数不等于δ函数，所以Bartlett 估计也是有偏估计即0≠Bias ，但当∞→N 时，0→Bias 。
由于我们假定各段周期图是相互独立的，所以可按式(2-8)得到下式：[][])(1)(ˆωωMxx I Var LP Var = ⎥⎥⎦⎤⎢⎢⎣⎡⎪⎭⎫ ⎝⎛+≈22sin )sin(1)(1ωωωM M P L xx(2-11)由此可见，随着L 的增加[])(ˆωxx P Var 是下降的，当∞→L 时，[]0)(ˆ→ωxxP Var 。
因此Bartlett 估计是一致估计。
比较式(2-10)的[])(ˆωxxP E 与式(2-1)的[])(ωN I E 可见在二种情况的估计量的期望值都是真值)(ωxx P 与窗口函数)(ωj B e W 的卷积形式，但后者将前者WB 中N 改为M ，N L N M <=/。
因而使)(ωj B e W 主瓣的宽度增大。