西南交大现代信号处理作业

1、现代信号处理作业1.(5)证明下面定理：任何一个无偏估计子方差的下界叫作Cramer-Rao下界定理：令为一样本向量，是的条件密度，若是的一个无偏估计子，且存在，则式中。其中是的某个不包含的正函数。2.(10)Wiener滤波是信号处理中最常用和基础的波形估计工具之一，对其在自己研究领域的应用情况进行一个简单综述。3.(5)二阶滑动平均过程由定义，式中表示正态分布，其均值为零、方差为。求x(n)的功率谱。4.(20)信号的函数表达式为：，其中，A(t)为一随时间变化的随机过程，dn(t)为经过390-410Hz带通滤波器后的高斯白噪声，n(t)为高斯白噪声，采样频率为1kHz，采样时间为2.0

2、48s。(1) 利用现代信号处理知识进行信号的谱估计；(2) 利用现代信号处理知识进行信号的频率提取；(3) 分别利用Wiener滤波和Kalman滤波进行去噪；(4) 利用Wigner-Ville分布分析信号的时频特征。5.(10)附件中表sheet1 为某地2008年4月28日凌晨12点至2008年5月4日凌晨12点的电力系统负荷数据，采样时间间隔为1小时，利用ARMA方法预测该地5月5日的电力系统负荷，并给出预测误差（5月5日的实际负荷数据如表sheet2）。1、定理：令为一样本向量，是的条件密度，若参数估计是真实参数的一个无偏估计子，且、存在，则的均方误差所能达到的下界（称为Crame

3、r-Rao下界）等于Fisher信息的导数，即： （1-1）不等式中等号成立的充分必要条件是： （1-2） 其中是的某个正函数，与样本无关。证明：由假设条件知，或，因此有 （1-3）对上式两边求关于的偏导，得即有 （1-4）另一方面，由复合函数的求导法，又有 （1-5）由于是的条件概率密度，故 （1-6）将式（1-5）和式（1-6）带入式（1-4），得或改写作 （1-7）由Cauchy-Schwartz不等式知，对于任意两个复函数和，恒有不定式：成立，并且当且仅当，等号成立，将Cauchy-Schwartz不等式应用于式（1-7），则有或等价为 （1-8） 由Cauchy-Schwartz不等

4、式等号成立的条件知：当且仅当，即式（1-2）成立时，不等式（1-8）才取等号。注意到，故有 （1-9）另由公式知 （1-10） 将式（1-9）和式（1-10）代入式（1-8），直接得到不等式（1-1）。根据前面的分析，不等式等号成立的充分必要条件是式（1-2）成立3、解：对取延迟形式：于是有：展开上式得到：即：对上式进行傅里叶变换，则有：从而得到：即4、music算法原理： music(Multiple Signal Classification)算法是针对多元天线阵测向问题提出的。假定M元的均匀线阵，阵元间距为d，信号的工作波长为。空间信号源共有D个，各信号不相关，各阵元的噪声互不相关，噪声

5、和信号也不相关。因此，第m个阵元的输出为 （4-1） 式中，为第个信号源的方向。将式（4-1）写成矩阵形式： （4-2） 式中：、。 求各阵元输出的相关矩阵，有： （4-3） （4-4） 式中：噪声的方差。 对式（4-3）的相关矩阵R作特征分解，其各特征值及其相对应的特征向量分别为：12+1 （4-5）12+1 （4-6）据式式（4-3），可得以下结论：（1）R的最小特征值等于，重数为(MD)，即+1=。据此，空间信号源的个数D可由下式得出：D=M-（R最小特征值的重数）。最小重数为1，因此，M阵元可测向的信号源数目的最大值为。（2）各特征向量相互正交。这些向量为矩阵列空间的基，由于最小特征值

6、为噪声的贡献，因此与最小特征值对应的那些特征向量所张成的子空间也是噪声的贡献，称之噪声子空间，记为。这样的列空间被划分成两个子空间,即信号子空间和噪声子空间： （4-7） （4-8） 由于各特征向量相互正交，故有：。在信号源所在方向上，诸方向向量，均处于信号子空间中，故。构造矩阵： （4-9） 显然有 music算法就是根据式（4-9）来求空间谱，有 （4-10） 谱峰所对应值就是信号源方向的估值。维纳滤波算法原理：维纳(Wiener)是用来解决从噪声中提取信号的一种过滤（或滤波）方法。这种线性滤波问题，可以看做是一种估计问题或一种线性估计问题。一个线性系统，如果它的单位样本响应为，当输入一个

7、随机信号，且： （4-1）其中：表示信号，表示噪声，则输出为： （4-2）我们希望通过线性系统后得到的尽量接近于，因此称为的估计值，用表示，即： （4-3）则维纳滤波器的输入输出关系可用下面图1表示。图4-1 维纳滤波器的输入输出关系实际上，式（4-2）所示的卷积形式可以理解为从当前和过去的观察值，来估计信号的当前值。因此，用进行过滤问题实际上是一种统计估计问题。一般地，从当前的和过去的观察值，估计当前的信号值成为过滤或滤波；从过去的观察值，估计当前的或者将来的信号值称为外推或预测；从过去的观察值，估计过去的信号值称为平滑或内插。因此维纳滤波器又常常被称为最佳线性过滤与预测或线性最优估计。这里

8、所谓的最佳与最优是以最小均方误差为准。如果我们分别以与表示信号的真实值与估计值，而用表示他们之间的误差，即： （4-4）显然可能是正值，也可能是负值，并且它是一个随机变量。因此，用它的均方误差来表达误差是合理的，所谓均方误差最小即它的平方的统计期望最小。卡尔曼滤波算法原理：卡尔曼滤波是基于状态空间方法的一套递推滤波算法，在状态空间方法中，引入了状态变量的概念。卡尔曼滤波的模型包括状态空间模型和观测模型。状态模型是反映状态变化规律的模型，通过状态方程来描写相邻时刻的状态转移变化规律；观测模型反映了实际观测量与状态变量之间的关系。Kalman滤波问题就是联合观测信息及状态转移规律来得到系统状态的最

9、优估计。假设动态系统的状态空间模型为 （4-1） （4-2）其中，X(t)为系统在时刻t的状态；Y(t)为对状态的观测值；W(t)为系统噪声，方差阵为Q；V(t)为观测噪声，方差阵为R；为状态转移矩阵；H为观测矩阵；为系统噪声驱动矩阵。卡尔曼滤波的计算流程为：计算状态估计值： （4-3）计算状态一步预测： （4-4）计算新息： （4-5）计算卡尔曼滤波增益： （4-6）计算一步预测均方误差： （4-7）计算一步预测估计均方误差： （4-8）下面给出卡尔曼滤波的系统模型框图：Matlab图形：图1 原始型号时域波形图图2 music方法对信号进行谱估计图3 信号频率提取图4 两种滤波方式时域结果

10、图图5 Wiener滤波去噪信号频域波形图6 Kalman滤波去噪信号频域波形图7 信号Wigner-Ville分布图8 信号伪Wigner-Ville分布Matlab程序：clc; clear allff1=100;ff2=300;ff3=200;fs = 1000; % 采样频率t = 0:1/fs:2.047; % 产生2048个采样时间点N=length(t);At=normrnd(0,1,1,2048); % 产生均值为0，方差为1的正态随机过程dnt1=wgn(1,2048,5) % 产生高斯白噪声%巴特沃斯带通滤波fs1=1000;Wp = 390 410/(fs1/2); %归

11、一化通带Ws = 350 450/(fs1/2);Rp = 3; %通带内的最大衰减Rs = 100; %阻带内的最小衰减n,Wn = buttord(Wp,Ws,Rp,Rs); %计算阶数和截至频率b,a = butter(n,Wn); %计算H(z)freqz(b,a,128,fs1) %画的幅频相频图dnt=filter(b,a,dnt1); %将产生的高斯白噪声进行巴特沃斯带通滤波器滤波 x1= sin(2*pi*ff1*t) +1.5*sin(2*pi*ff2*t)+At.*sin(2*pi*ff3*t)+dnt; % 原始信号x=awgn(x1,5); % 给原始信号加入信噪比为5

12、dB的高斯白噪声 % 信号时域图figureplot(t,x);title(观测信号时域波形);xlabel(时间/s); ylabel(幅值);axis tight; figuresubplot(211);plot(t,x);title(观测信号时域波形);xlabel(时间/s);ylabel(幅度);axis tight;subplot(212);plot(t,x1);title(原始信号时域波形);xlabel(时间/s);ylabel(幅度);axis tight; % 信号频域图Y = fft(x,N); %做FFT变换Ayy=(abs(Y)*2/N; %为了与真实振幅对应，需要将

13、变换后结果乘以2除以N。Ayy(1)=Ayy(1)/2; %直流量不同与其他f=(1:N-1)*fs/N; %换算成实际的频率值 ,等同于F=(0:N-1)*Fs/N;figure;plot(f(1:N/2),Ayy(1:N/2); title(N=1024 信号幅度-频率曲线图);xlabel(频率/Hz);ylabel(幅值);axis tight; % music法信号频谱图nfft=1024;figurepmusic(x,7,1.1,nfft,fs,32,16); title(music法对信号进行谱估计); xlabel(f/kHz); ylabel(功率/db);axis tigh

14、t; % wiener滤波去噪M=500; %维纳滤波器长度Rxx=xcorr(x(1:1024); %计算观测信号的自相关函数A=Rxx(1024:1523);Rx1x=xcorr(x1(1:1024),x(1:1024); %计算观测信号与期望信号的互相关函数B=Rx1x(1024:1523); y=wienerfilter(x,A,B,M); %用维纳滤波进行去噪 figuresubplot(2,1,1)plot(t,x);title(观测信号时域波形);xlabel(时间);ylabel(幅度);axis tightsubplot(2,1,2)plot(t,y);title(经过维纳滤

15、波器滤波后的信号时域波形);xlabel(时间);ylabel(幅度);axis tight Y = fft(y,N); %做FFT变换Ayy=(abs(Y)*2/N; %为了与真实振幅对应，需要将变换后结果乘以2除以N。Ayy(1)=Ayy(1)/2; %直流量不同与其他f=(1:N-1)*fs/N; %换算成实际的频率值 ,等同于F=(0:N-1)*Fs/N;figureplot(f(1:N/2),Ayy(1:N/2); title(经过维纳滤波后信号的频域图形);xlabel(频率/Hz);ylabel(幅值);grid on; % Kalman滤波去噪a=1; %a为状态转移阵w=0.

16、1*randn(1,N); %w为过程噪声V=0.02*randn(1,N);%V为测量噪声q1=std(V); Rvv=q1.2; q2=std(x); Rxx=q2.2; q3=std(w); Rww=q3.2; %Rvv、Rww分别为过程噪声和测量噪声的协方差c=1; %c为测量矩阵Y=c*x+V; %测量方差p(1)=10; %初始最优化估计协方差s(1)=100; %初始最优化估计for t1=2:Np1(t1)=a.2*p(t1-1)+Rww; %进一一步估计的协方差，从t-1时刻最优化估计s的协方差得到t-1时刻到t时刻一步估计的协方差b(t1)=c*p1(t1)/(c.2*p1

17、(t1)+Rvv); %b为卡尔曼增益，表示为状态误差的协方差与量测误差的协方差之比s(t1)=a*s(t1-1)+b(t1)*(Y(t1)-a*c*s(t1-1); %称之为新息，是观测值与一步估计得到的观测值之差，此式由上一时刻状态的最优化估计s(t-1)得到当前时刻的最优化估计s(t)p(t1)=p1(t1)-c*b(t1)*p1(t1); %此式由一步估计的协方差得到此时刻最优化估计的协方差end Y = fft(s,N); %做FFT变换Ayy=(abs(Y)*2/N; %为了与真实振幅对应，需要将变换后结果乘以2除以N。Ayy(1)=Ayy(1)/2; %直流量不同与其他f=(1:

18、N-1)*fs/N; %换算成实际的频率值 ,等同于F=(0:N-1)*Fs/N;figureplot(f(1:N/2),Ayy(1:N/2); title(经过卡尔曼滤波后信号的频域图形);xlabel(频率/Hz);ylabel(幅值);grid on; % tfrwv分布figure;tfr,t,f=tfrwv(x);contour(t/fs,f*fs,tfr);title(维格纳-威利分布);xlabel(时间 t);ylabel(频率 f);figure;contour(t/fs,f*fs,abs(tfr);title(伪维格纳-威利分布);xlabel(时间 t);ylabel(频

19、率 f)function y=wienerfilter(x,Rxx,Rxd,N)%进行维纳滤波%x是输入信号，Rxx是输入信号的自相关向量%Rxx是输入信号和理想信号的的互相关向量，N是维纳滤波器的长度%输出y是输入信号通过维纳滤波器进行维纳滤波后的输出h=yulewalker(Rxx,Rxd,N);%求解维纳滤波器系数t=conv(x,h); %进行滤波Lh=length(h); %得到滤波器的长度Lx=length(x); %得到输入信号的长度y=t(double(uint16(Lh/2):Lx+double(uint16(Lh/2)-1);%输出序列y的长度和输入序列x的长度相同 %以下

20、是维纳滤波器系数的求解function h=yulewalker(A,B,M) %求解Yule-Walker方程%A是接收信号的自相关向量为 Rxx(0),Rxx(1),.,Rxx(M-1)%B是接收信号和没有噪声干扰信号的互相关向量为 Rxd(0),Rxd(1),.,Rxd(M-1)%M是滤波器的长度%h保存滤波器的系数%例如A=6,5,4,3,2,1;B=100,90,120,50,80,200;M=6;%求解出h=26.4286 -20.0000 50.0000 -50.0000 -45.0000 81.4286T1=zeros(1,M);%T1存放中间方程的解向量T2=zeros(1,

21、M);%T2存放中间方程的解向量T1(1)=B(1)/A(1);T2(1)=A(2)/A(1);X=zeros(1,M);for i=2:M-1temp1=0;temp2=0; for j=1:i-1 temp1=temp1+A(i-j+1)*T1(j); temp2=temp2+A(i-j+1)*T2(j); end X(i)=(B(i)-temp1)/(A(1)-temp2); for j=1:i-1 X(j)=T1(j)-X(i)*T2(j); end for j=1:i T1(j)=X(j); endtemp1=0;temp2=0; for j=1:i-1 temp1=temp1+A(

22、j+1)*T2(j); temp2=temp2+A(j+1)*T2(i-j); end X(1)=(A(i+1)-temp1)/(A(1)-temp2); for j=2:i X(j)=T2(j-1)-X(1)*T2(i-j+1); end for j=1:i T2(j)=X(j); endendtemp1=0;temp2=0;for j=1:M-1 temp1=temp1+A(M-j+1)*T1(j); temp2=temp2+A(M-j+1)*T2(j);endX(M)=(B(M)-temp1)/(A(1)-temp2);for j=1:M-1 X(j)=T1(j)-X(M)*T2(j);

23、endh=X;function tfr,t,f = tfrwv(x,t,N,trace);%TFRWV Wigner-Ville time-frequency distribution.% TFR,T,F=TFRWV(X,T,N,TRACE) computes the Wigner-Ville distribution% of a discrete-time signal X, % or the cross Wigner-Ville representation between two signals. % % X : signal if auto-WV, or X1,X2 if cross-

24、WV.% T : time instant(s) (default : 1:length(X).% N : number of frequency bins (default : length(X).% TRACE : if nonzero, the progression of the algorithm is shown% (default : 0).% TFR : time-frequency representation. When called without % output arguments, TFRWV runs TFRQVIEW.% F : vector of normal

25、ized frequencies.% Example :% sig=fmlin(128,0.1,0.4); tfrwv(sig);% % See also all the time-frequency representations listed in% the file CONTENTS (TFR*) % F. Auger, May-August 1994, July 1995.% Copyright (c) 1996 by CNRS (France).% - CONFIDENTIAL PROGRAM - % This program can not be used without the

26、authorization of its% author(s). For any comment or bug report, please send e-mail to % if (nargin = 0), error(At least one parameter required);end;xrow,xcol = size(x); if (nargin = 1), t=1:xrow; N=xrow ; trace=0;elseif (nargin = 2), N=xrow ; trace=0;elseif (nargin = 3), trace = 0;en

27、d; if (N2), error(X must have one or two columns);elseif (trow=1), error(T must only have one row); elseif (2nextpow2(N)=N), fprintf(For a faster computation, N should be a power of twon);end; tfr= zeros (N,tcol); if trace, disp(Wigner-Ville distribution); end;for icol=1:tcol, ti= t(icol); taumax=mi

28、n(ti-1,xrow-ti,round(N/2)-1); tau=-taumax:taumax; indices= rem(N+tau,N)+1; tfr(indices,icol) = x(ti+tau,1) .* conj(x(ti-tau,xcol); tau=round(N/2); if (ti=tau+1), tfr(tau+1,icol) = 0.5 * (x(ti+tau,1) * conj(x(ti-tau,xcol) + . x(ti-tau,1) * conj(x(ti+tau,xcol) ; end; if trace, disprog(icol,tcol,10); e

29、nd;end; tfr= fft(tfr); if (xcol=1), tfr=real(tfr); end ; if (nargout=0), tfrqview(tfr,x,t,tfrwv);elseif (nargout=3), f=(0.5*(0:N-1)/N);end;5、ARMA模型算法： ARMA模型时间序列分析法简称为时序分析法，是一种利用参数模型对有序随机振动响应数据进行处理，从而进行模态参数识别的方法。参数模型包括AR自回归模型、MA滑动平均模型和ARMA自回归滑动平均模型。 N个自由度的线性系统激励与响应之间的关系可用高阶微分方程来描述，在离散时间域内，该微分方程变成由一系

30、列不同时刻的时间序列表示的差分方程，即ARMA时序模型方程： （5-1） 式（5-1）表示响应数据序列与历史值的关系，其中等式的左边称为自回归差分多项式，即AR模型，右边称为滑动平均差分多项式，即MA模型。N为自回归模型和滑动均值模型的阶次，、分别表示待识别的自回归系数和滑动均值系数，表示白噪声激励。当k0时，设。 由于ARMA过程具有唯一的平稳解为 （5-2）式中：为脉冲响应函数。 的相关函数为 （5-3） 是白噪声，故 （5-4）式中：为白噪声方差。 将此结果代人式（5-3），即可得 （5-5） 因为线性系统的脉冲响应函数，是脉冲信号，激励该系统时的输出响应，故由ARMA过程定义的表达式为

31、 （5-6） 利用式（5-5）和式（5-6），可以得出： （5-7） 对于一个ARMA过程，当是大于其阶次2N时，参数0。故当l2N时，式（5-7）恒等于零，于是有 （5-8）或写成 （5-9） 设相关函数的长度为L，并令M2N。对应不同的l值，由代人以上公式可得一组方程： （5-10）或缩写为矩阵形式: （5-11） 式（11）为推广的Yule-walker方程。一般情况下，由于L比2N大得多，采用伪逆法可求得方程组的最小二乘解，即 （5-12）由此求得自回归系数。滑动平均模型系数可通过以下非线性方程组来求解： （5-13）其中 （5-14）式中：为响应序列的自协方差函数。 当求得自回归系数

32、和滑动均值系数后，可以通过ARMA模型传递函数的表达式计算系统的模态参数，ARMA模型的传递函数为 （5-15）ARMA模型Matlab图形：图1 电力系统负荷数据曲线和平稳化处理后数据曲线图2 电力系统负荷数据自相关和偏相关函数曲线图3 使用ARMA模型对电力系统负荷进行预测曲线图4 预测误差曲线ARMA模型Matlab程序：clc;clear;x1=xlsread(C:Documents and SettingsAdministrator桌面现代信号处理作业负荷数据.xls,Sheet1);x1=x1(:,2);x2=xlsread(C:Documents and SettingsAdmi

33、nistrator桌面现代信号处理作业负荷数据.xls,Sheet2);x2=x2(:,2);x=x1;x2; % 电力系统负荷数据，前168个点为历史数据，后24个点为待预测数据N1=length(x1);N=length(x);t=1:N;figure;plot(x);title(电力系统负荷数据);xlabel(时间);ylabel(电力系统负荷);axis tight; %对原始数据进行去趋势处理，即零均值化、平稳化处理z1=diff(x); %差分Y=z1-mean(z1); %零均值H,PValue,TestStat,CriticalValue = adftest(Y); %检验是

34、否为时间平稳序列 figure;plot(Y,r)title(平稳化处理后电力系统负荷数据);xlabel(时间),ylabel(功率); %计算自相关函数、偏相关函数figuresubplot(2,1,1)autocorr(Y); %计算置信度为95%的acf，并画出其自相关函数曲线；subplot(2,1,2) parcorr(Y); %计算置信度为95%的pacf，并画出其偏自相关函数曲线； %自相关系数和偏相关系数均拖尾性，所以选择arma模型对电力系统负荷数据进行预测 %估计arma模型参数，以AIC标准来定阶z=iddata(Y);test = ; for p = 1:10 %自回归对应PACF,给定滞后长度上限p和q for q = 1: