功率谱和自适应滤波分析.docx

1、功率谱和自适应滤波分析实验一：几种功率谱分析方法的介绍及其性能比较摘 要摘要：谱分析所考虑的问题是，通过非参数或者参数化的技术从有限测量数据中确定时间序列的谱内容(即功率随频率分布的情况)。功率谱分析在很多领域都有广泛的利用，比如语音信号处理、军用雷达信号处理、无线通信等领域。基于零极点信号模型，本实验通过对不同谱分析方法的分析和数值仿真，得到了各种方法的谱分析性能。通过对一多频率信号的数值仿真发现，经典的谱分析方法的分辨率普遍比现代谱分析方法的分辨率低。最后，在本文末尾的附录部分给出了Matlab下的主要仿真代码。关键字：功率谱分析；零极点信号模型；分辨率；1 引言在信号系统中，往往需要研究

2、具有特定统计特性的随机信号。由于随机信号是一类持续时间无限长，具有无限大能量的功率信号，它不满足傅里叶变换条件，而且也不存在解析表达式，因此就不能够应用确定信号的频谱计算方法去分析随机信号的频谱。然而，虽然随机信号的频谱不存在，但其相关函数是可以确定的。如果随机信号是平稳的，那么其相关函数的傅里叶变换就是它的功率谱密度函数，简称功率谱。功率谱反映了单位频带内随机信号的一个样本信号来对该随机过程的功率谱密度函数做出估计。序列功率谱分析是数字信号处理中的主要内容之一，通过对信号序列的一段时间观测从而利用相对应的方法来得到信号对应的频谱成分等信息。一般地，信号模型是零极点模型。谱分析的方法有很多种，

3、根据分别率的高低有低分辨率和高分辨率两类。低分辨率的谱分析方法主要有周期图法，相关图法和改进的相关图法等；高分辨率的谱分析方主要有Capon法、最大熵法、Borgiotti-Laguuas法、Music法、min-norm法等。其中周期图法，相关图法和改进的相关图法、Capon法、最大熵法又叫做非参数法，Borgiotti-Laguuas法、Music法、min-norm法称为参数法。可以看出谱分析的方法有很多种，所以需要针对一些被广泛研究的和代表性的方法来进行分析。因此，本文主要分析和仿真验证了周期图法，相关图法、Capon法、Music法、min-norm法、改进的相关图法原理和性能。本文

4、的基本结构是：第二部分介绍一下各种谱分析方法所基于的信号模型，并对一些参数和默认值进行标准化和预设。第三部分是原理介绍，主要分析各种方法理论知识和关键表达式的推导，最终得到功率谱表达式。各种方法的数值仿真结果及其分析在第四部分给出。接下来是对本文的一个总结。最后给出了主要程序附录。2 信号模型在本文中一般的信号模型可以表示为 其中和分别表示不同的幅度和频率，这里假设，；假设噪声和相位分别服从标准复高斯分布和0到的均匀分布；采样信号，是总采样点数；表示采样频率，这里设为。对于Music法其信号模型为无线通信中多天线接收模型，可以表示成 其中是接收信号；是发送端的接收信号；是噪声；信道增益，可以表

5、示为 其中表示接收天线的个数。假设接收信号是广义平稳的；假设服从零均值方差为1的高斯分布，且和是相互独立的；假设协方差矩阵 是可逆的，同时。3 原理介绍3.1 周期图方法原理Schuster于1899年首先提出周期图法，也称直接法，取平稳随机信号的有限个观察值x(0)，x(1)，x（N-1）对功率谱S（）进行估计 主要性能指标有：3.1.1 估计的均值 W()是窗函数（n）的Fourier变换。当N时，是无偏估计；N是有限值时，是有偏估计，偏差为 3.1.2 估计的方差 其中D0（）是矩形窗d0（n）=1，的Fourier变换，可见不是的一致估计；随着N的增大，谱估计起伏增大，N时，。3.1.

6、3估计的分辨率数据窗为长度为N的矩形窗时，N增大时，分辨率会提高，但会使的起伏加剧，可见方差与分辨率是一对矛盾。周期图法应用比较广泛，主要是由于它与序列的频谱有对应关系，可以采用FFT快速算法来计算。但是，这种方法需要对无限长的平稳序列进行截断，相当于对其加矩形窗，使之成为有限长数据。同时，这也意味着对自相关函数加三角窗，使功率谱与窗函数卷积，从而产生频谱泄漏，容易使弱信号的主瓣被强信号的旁瓣所淹没，造成频谱的模糊和失真，使得谱分辨率较低。经典法功率谱估计是间接法功率谱估计的一个特例，又称为周期图法，它是把随机序列信号的N个观测数据视为一能量有限的序列，接着计算的离散傅里叶变换，得，然后取其幅

7、值的平方，并除以N，作为对功率谱的估计。 3.2 自相关函数法原理（Blackman-Tukey 法）假设已知随机信号的N个观测样本，则其自相关序列可以用下式进行估计 当仅使用长度为2M-1的自相关序列时，对其进行傅立叶变换即可得到功率谱估计如下 上式可以通过对自相关序列加窗表示如下 其中M为加窗长度，为矩形窗函数，定义如下 因此，在一定程度上可以看作是“真正的功率谱P(w)”与窗函数傅立叶变换的卷积。 矩形窗函数不仅降低了谱估计的分辨率，而且使谱估计产生了旁瓣，旁瓣效应使那些处于旁瓣附近功率较小的频率分量被淹没。为了降低旁瓣的影响，可以采用具有较小旁瓣的窗函数，如Hamming窗，定义为：

8、这种窗函数可以有效的抑制旁瓣，但此时主瓣宽度增大，从而降低了谱估计的分辨率。这种主瓣和旁瓣的矛盾在非参数化功率谱估计方法中是无法解决的。3.3 Capon 方法原理Capon滤波器组谱估计方法就是设计一种有限冲击响应（FIR）滤波器，使得滤波器输入信号中的某个频率成分完全通过的同时，使滤波器的输出功率最小，进而将x(n)的谱估计转换成求高斯白噪声下谐波信号的幅度或功率问题。假设x(n);n=0,1,2N-1为平稳随机过程的N点采样数据，h（n）,n=0,1,2,N-1为N阶FIR滤波器的冲激响应函数，则过程x（n）经过FIR滤波器的输出y（k）为 h(n),n=0,1n-1为n阶FIR滤波器的

9、冲击响应函数， 对频率为w的谐波分量的响应由滤波器的频率响应来确定，即 Capon滤波器组设计目的是使观测数据x（n）无失真的通过以频率w为中心的带通滤波器，其思想在数学上可表示为下列优化问题 用拉格朗日数乘法来求最优解 对式(18)中的求偏导得到 将代入约束条件中，我们可以得到 所以， 因此，从Capon方法得到功率谱估计 最后,估计的谱密度为 3.4 Music 方法原理Music算法是一种基于矩阵特征空间分解的方法。从几何角度讲，信号处理的观测空间可以分解为信号子空间和噪声子空间，显然这两个空间是正交的。信号子空间由阵列接收到的数据协方差矩阵中与信号对应的特征向量组成，噪声子空间由协方差

10、矩阵中所有最小特征值（噪声方差）对应的特征向量组成。Music算法的核心就是利用这两个互补空间之间的正交特性来估计空间信号的方位。噪声子空间的所有向量被用来构造谱，所有空间方位谱中的峰值位置对应信号的来波方位。下面重点来分析Music算法的原理，根据式(2)和(3)对接收信号计算协方差矩阵R 其中，是满秩的对角矩阵。如果对协方差矩阵进行特征值分解，可以看到在信号子空间由 同理噪声子空间有 所以信号子空间对应的特征向量和噪声子空间对应特征向量可以写成 因此 由上式可以看出。根据算法的核心思想估计的频率点可以从下式中求出 其中，。最后Music算法的谱估计可以写成 其中 3.5 Min-norm

11、法和Music法一样，Min-norm法也是一种高分辨率的谱分析方法，并且拥有较强的鲁棒性。在Music法基础上，这种方法可以描述成以下优化问题 我们知道 所以，式(32)可以重写为 因此，可以利用拉格朗日函数法 利用式(32)的第一个限制条件，可以解出拉格朗日乘子 所以，这种方法的最优解为 回顾以下Min-norm法的功率谱表达式 其中向量有以下特点：1，它的第一个元素为1和最小范数；2，属于噪声样本空间。最后，将式(37)得到的解代入式(38)可以得到 由于分子是一个常数，对最后的谱密度形状不会产生影响，一次略去分子式(39)可以重新写成 4 仿真分析本节给出的是本实验基于MATLAB的数

12、值仿真结果。信号所在的频率点分别设为0.1Hz,0.15Hz和0.4Hz,采样点数为100。以下是Capon法、Blackman-Tukey法、Music法、min-norm法，periodogram法，correlogram法在不同信噪比条件下的功率谱图形。图1-1 给出的是采样点数为100，信噪比为100dB时的功率谱图形，从图中可以看出periodogram法，Blackman-Tukey法以及min-norm法在非信号频率点上功率仍有较大波动，并且Blackman-Tukey法通过加窗并没有改善谱性能。Capon法利用对其它非信号频率点功率的强抑制性能，使得其功率谱曲线特别光滑。图2-

13、2 是采样点数为100，信噪比为100dB时的功率谱图形，很显然随着信噪比的降低，传统的谱分析方法periodogram法，Blackman-Tukey法和correlogram法的分别率下降显著。到信噪比下降到0dB时，图2-3给出了各种方法的功率谱图形。其中Capon法、Music法、min-norm法的性能还比较稳定，仍能稳健的分辨出各频率分量。图1-1信噪比为100dB，采样点为100时，各方法得到的功率谱图1-2 信噪比为10dB，采样点为100时，各方法得到的功率谱图1-3信噪比为0dB，采样点为100时，各方法得到的功率谱5 总结本文分析了几种重要功率谱分析方法的原理和数学推导，

14、并用仿真结果验证分析的正确性。发现相比传统的功率谱分析方法现代谱估计算法拥有较好的分别率，也就是说基于样本协方差矩阵得到的功率谱相比于直接利用采样点或相关序列得到的功率谱，其谱拥有曲线光滑、分辨率高、噪声功率低等特点。当然，这些算法还不能满足实际中的应用，比如Capon法对角度非常敏感，当角度出现一点偏差时其效果会变得很差。所以，一些改进的算法也逐渐被提出，以便解决实际应用中出现的问题。实验二：基于LMS算法的自适应滤波1 系统模型考虑一种AR模型，其中 w(n)WGN(0,)。当信噪比较小时如何从中提取出我们感兴趣的信号。自适应滤波是常采用的一种方法，其基本实现框图如图2-1所示图2-1 自

15、适应滤波一般模型2 LMS 自适应滤波器原理自适应滤波器由参数可调的数字滤波器和自适应算法两部分组成，参数可调的数字滤波器可以是FIR滤波器或者IIR滤波器，也可以是格型数字滤波器。自适应滤波器的基本原理图如图1所示。在闭环自适应滤波器中，输入信号x(n)通过参数可调数字滤波器后产生输出信号y(n)，将其与参考信号d(n)进行比较，形成误差信号e(n)。e(n)通过某种自适应算法对滤波器参数进行调整，最终使e(n)的均方值最小。因此，实际上自适应滤波器是一种能够自动调整本身参数的特殊维纳滤波器，在设计的时候不需要事先知道关于输入信号和噪声统计特性的知识，它能够在自己的工作过程中逐渐“了解”或估

16、计出所需的统计特性，并以此为依据自动调整自己的参数，以达到最佳滤波效果。一旦输入信号的统计特性发生变化，他又能够跟踪这些变化，自动调整参数，使滤波器性能重新达到最佳。当输入过程的统计特性未知时，自适应滤波器调整自己参数的过程称为“学习”过程；当输入过程的统计特性变化时，自适应滤波器调整自己参数的过程称为“跟踪”过程。图2-1所示的自适应滤波器的原理图有两个输入：x(n)和d(n)，两个输出：y(n)和e(n)，均为时间序列。其中x(n)可以是单输入信号，也可以是多输入信号。在不同的应用背景下这些信号代表不同内容。输入矢量为 ()加权矢量（即滤波器参数矢量）为 ()滤波器的输出为 ()y(n)相

17、对于滤波器期望输出d(n)的误差为 ()根据最小均方误差准则，最佳的滤波器参量应使得均方误差为最小，在W(n)是常数矢量的情况下，在时刻n的均方误差表达式为 ()其中，是期望响应d(n)的方差，是输入矢量和期望响应d(n)的互相关矢量，是输入矢量X(n)的自相关矩阵。上式对W(n)求偏导，并令导数等于零，可得到正则方程 ()当为满秩时，正则方程有唯一解 ()这个解称为维纳解。当时，均方误差性能函数的最小值（即最小均方误差）等于 ()以上直接对求逆来得到的方法运算量很大，实际中常用到梯度法。3 仿真结果如图2-2所示，表示两个系数随着迭代次数的增加，逐渐逼近设定值的学习曲线。a中的两个系数，随着

18、迭代次数的增加，分别收敛到原信号的两个设定值系数的大小。b中是不同的 对模型参数收敛的影响，可以看出 的收敛性要优于。图2-3给出的是基于自适应滤波的线性预测仿真图，从中可以看出随着迭代次数的增加其预测的准确度逐渐增加。图2-2 LMS算法的收敛图形图2-3 基于自适应滤波的线性预测4 总结自适应滤波在信号处理中有着广泛的应用，如系统辨识、自适应均衡、噪声（干扰）消除以及线性预测。本文对自适应滤波进行了简单探讨，分析了其模型和数学原理，最后给出具体实现方法。并对其在线性预测方面的应用进行了仿真。5 附录5.1程序清单一：谱估计初始版本：2016/3/9clc;close all; clear

19、all;f1 = 0.1; f2 = 0.15;f3 = 0.4; %frequecy unit:(khz)SNRdB = -3;N0 = 10(-SNRdB/20);fs = 0.5; %sample rate khzTs = 1/fs; M = 10; %capon filter lengthN = 100; %samplesD = 1/(N*Ts); %minimun resolution of periodogramn = 0:N-1; %samples numbernoise = (randn(1,N)+1i*randn(1,N)/sqrt(2);x = exp(1i*2*pi*f1

20、*n*Ts)+exp(1i*2*pi*f2*n*Ts)+exp(1i*2*pi*f3*n*Ts)+N0*noise;x1 = exp(1i*2*pi*f1*n)+exp(1i*2*pi*f2*n)+exp(1i*2*pi*f3*n)+N0*noise;%=periodogram=P_peri = 10*log10(abs(fft(x,N).2/Ts);figure(1);subplot(131);box on;plot(D*n,P_peri);set(gca,FontSize,20);title(periodogram ,FontSize,20);xlabel(Hz);ylabel();xli

21、m(0,0.5);%=correlogram=corx = xcorr(x,unbiased);P_cor=10*log10(abs(fft(corx,N);subplot(132);box on;set(gca,FontSize,20);plot(D*n,P_cor);xlabel(Hz);title(correlogram ,FontSize,20);xlim(0,0.5);%=Blackman-Tukey=K = N/4;wd = hamming(K);corx = xcorr(x,unbiased);corx = corx(1:K);xx = wd.*corx;P_bt=10*log1

22、0(abs(fft(xx,N);subplot(133);box on;plot(D*n,P_bt);set(gca,FontSize,20);title(Blackman-Tukey ,FontSize,20);xlabel(Hz);xlim(0,0.5);%=capon=P_capon,f = capon(x1,M,N);figure(2);subplot(131);box on;plot(f,P_capon);set(gca,FontSize,20);xlabel(Hz);ylabel();title(capon ,FontSize,20);set(gca,FontSize,20);%=

23、music=fest,P_music,f = music(x1,3,N,N);subplot(132);box on;set(gca,FontSize,20);plot(f,P_music);str = music ,num2str(fest);display(sprintf(str);title(music ,FontSize,20);xlim(0,0.5);xlabel(Hz);%=min-norm=fest1,P_min = minnorm(x1,3,N,N);subplot(133);box on;plot(f,P_min);str = minnorm ,num2str(fest1);

24、display(sprintf(str);set(gca,FontSize,20);title(minnorm ,FontSize,20);xlim(0,0.5);xlabel(Hz);更新版2016/10/19clc;close all; clear all;%the correlogram method is outpreform preiodogram method when the samples%are large enough.f = 0.1 0.3 0.4; %frequecy unit:(khz)Len_f = length(f);SNRdB = 10; %signal-to-

25、noise ratioSNR_Lin = 10(-SNRdB/20);fs = 1; %sample rate khzTs = 1/fs; M = 10; %filter length for capon method in spectral analysisN = 100; % samples (using small data, so the methods preformance of preiodogram and correlogram is equal )K = N-2; %Blackman-Tukey parameterD_peri = 1/(N*Ts); % the minim

26、un resolution of periodogramD_cor = 1/(2*N-1)*Ts); % the minimun resolution of correlogramD_BT = 1/(2*K-1)*Ts); % the minimun resolution of correlogramn = 0:N-1; %samples linespacenoise = (randn(1,N)+1i*randn(1,N)/sqrt(2); x = exp(1i*2*pi*f(1)*n*Ts)+exp(1i*2*pi*f(2)*n*Ts)+exp(1i*2*pi*f(3)*n*Ts)+SNR_

27、Lin*noise; %generate testing data=periodogram=P_peri = (abs(fft(x,N).2)/N;x_peri, ix_peri = sort(P_peri);Est_peri = ix_peri(end-Len_f+1:end)*D_peri;figure(1);subplot(131);box on;plot(D_peri*n,P_peri);set(gca,FontSize,20);grid ontitle(periodogram ,f=, num2str(Est_peri),FontSize,10);xlabel(Hz);ylabel(

28、);axis tight;% =correlogram=corx = xcorr(x,unbiased);P_cor = abs(fft(corx);x_cor, ix_cor = sort(P_cor);Est_cor = ix_cor(end-Len_f+1:end)*D_cor;subplot(132);box on;set(gca,FontSize,20);plot(D_cor*(0:2*N-2),P_cor);xlabel(Hz); title(correlogram ,f=, num2str(Est_cor),FontSize,10);axis tight;grid on% =Blackman-Tukey=wd = hamming(2*K-1);x_