运用matlab进行随机信号的功率谱密度估计仿真【急求】

好眼熟啊。。。通院08级的随机信号实验题目？

这份是我在网络上找的：

（1）周期图法： 思想：周期图法是为了得到功率谱估值，先取信号序列的离散傅里叶变换，然后取其幅频特性的平方并除以序列长度 N。由于序列 x(n)的离散傅里叶变换 X(k) 具有周期性，因而这种功率谱也具有周期性，常称为周期图。周期图是信号功率谱的一个有偏估值；而且，当信号序列的长度增大到无穷时，估值的方差不趋于零。因此，随着所取的信号序列长度的不同，所得到的周期图也不同，这种现象称为随机起伏。由于随机起伏大，使用周期图不能得到比较稳定的估值。

程序：

%首先，生成输入信号的程序为：

clear ; fs=20000; n=0:1/ fs: 0 1; N=lengt h( n); W=2000pi;%因方波频率 F=1000HZ 所以角频率 W=2000pi X1n=square( Wn) ;%方波信号 X2n=randn( 1, N);%白噪声信号 xn=X1n+0 2X2n; subplot (3, 1, 1) plot (n, xn) ; xlabel( ' n' ) ylabel('输入信号') %其次，开始用周期图法进行估计； clear all fs=20000; n=0:1/fs:01; N=length(n); W=2000pi; x1n=square(Wn); x2n=randn(1,N); xn=x1n+02x2n; subplot(2,1,1) plot(n,xn); Nfft=256;N=256;%傅里叶变换的采样点数256 Pxx=abs(fft(xn,Nfft)^2)/N; f=(0:length(Pxx)-1)fs/length(Pxx); subplot(2,1,2), plot(f,10log10(Pxx)),%转成DB 单位

（2）自相关函数法： 思想： 随机信号 x(n)的相关函数是在时间域内描述随机过程的重要特征。自相关函数是随机信号在不同时刻的值之间的依赖性的量度，是一个很有用的统计平均量，其定义为自相关函数 (1) 式中E·表示数学期望，表示共轭值,m 为时间滞后数。在随机信号处理中，自相关函数可以用来检测淹没在随机噪声干扰中的信号，随机信号的自功率谱等于它的自相关函数的傅里叶变换。因此，通过自相关估计可求得信号的功率谱。利用计算机计算自相关估值有两种方法。一种是直接方法，先计算出随机信号和它的滞后序列的乘积，再取其平均值即得相关函数的估计值。另一种是间接方法，先用快速变换算法计算随机序列的功率谱密度，再作反变换计算出相关函数。

程序： n=0:1/fs:1; N=length(n); W=2000pi; x1n=square(Wn); x2n=randn(1,N); xn=x1n+02x2n; figure(2) subplot(3,1,1) plot(n,xn);%输入信号 axis([0 001 -5 5]) m =-100:100 [r,lag]=xcorr(xn,100,'biased')%求XN 的自相关函数R，biased 为有偏估计,lag 为R 的序列号 subplot(3,1,2) hndl=stem(m,r);%绘制离散图，分布点从-100—+100 set(hndl,'Marker','') set(hndl,'MarkerSize',2); ylabel('自相关函数R(m)') %利用间接法计算功率谱 k=0:1000;%取1000 个点 w=(pi/500)k; M=k/500; X=r(exp(-lipi/500)^(m'k));%对R 求傅里叶变换,li 就是j magX=abs(X); subplot(3,1,3) plot(M,10log10(magX)); xlabel('功率谱的改进直接法估计')

（3）自协方差法： 思想：在实际中,有时用随机信号在两个不同时刻t 1 ,t 2 的取值X(t 1 )和X(t 2 )之间的二阶混合中心矩来描述随机信号 X(t)在任意两个时刻取值起伏变化的相依程度, 即自协方差。自协方差函数与自相关函数描述的特性是一致的，所以其原理与相关函数法近似。

clear all; fs=20000; n=0:1/fs:01; P=2000pi; y=square(Pn); xn=y+02randn(size(n));%绘制信号波形 figure(3) subplot(2,1,1) plot(n,xn) xlabel('时间（s）') ylabel('幅度') title('y+randn(size(n))') ymax_xn=max(xn)+02; ymin_xn=min(xn)-02; axis([0 03 ymin_xn ymax_xn]) %使用协方差法估计序列功率谱 p=floor(length(xn)/3)+1; nfft=1024; [xpsd,f]=pcov(xn,p,nfft,fs,'half'); %绘制功率谱估计 pmax=max(xpsd); xpsd=xpsd/pmax; xpsd=10log10(xpsd+0000001); subplot(2,1,2) plot(f,xpsd) title('基于协方差的功率谱估计') ylabel('功率谱估计（db）') xlabel('频率') grid on; ymin=min(xpsd)-2; ymax=max(xpsd)+2; axis([0 fs/2 ymin ymax])

（4）最大熵法 思想：上网查的最大熵法原理是取一组时间序列，使其自相关函数与一组已知数据的自相关函数相同，同时使已知自相关函数以外的部分的随机性最强，以所取时间序列的谱作为已知数据的谱估值。它等效于根据使随机过程的熵为最大的原则,利用 N 个已知的自相关函数值来外推其他未知的自相关函数值所得到的功率谱。最大熵法功率谱估值是一种可获得高分辨率的非线性谱估值方法，特别适用于数据长度较短的情况。最大熵法谱估值对未知数据的假定 ，一个平稳的随机序列，可以用周期图法对其功率谱进行估值。这种估值方法隐含着假定未知数据是已知数据的周期性重复。现有的线性谱估计方法是假定未知数据的自相关函数值为零，这种人为假定带来的误差较大。最大熵法是利用已知的自相关函数值来外推未知的自相关函数值，去除了对未知数据的人为假定，从而使谱估计的结果更为合理。熵在信息论中是信息的度量，事件越不确定，其信息量越大，熵也越大。对于上述问题来说，对随机过程的未知的自相关函数值，除了从已知的自相关函数值得到有关它的信息以外，没有其他的先验知识。因而，在外推时，不希望加以其他任何新的限制，亦即使之“最不确定”。换言之，就是使随机过程的熵最大。

程序： fs=20000; n=0:1/fs:1; N=length(n); W=2000pi; x1n=square(Wn); x2n=randn(1,N); xn=x1n+02x2n; figure(5) subplot(3,1,1) plot(n,xn); asis([0 001 -5 5]) Nfft=256;%分段长度256 [Pxx,f]=pmem(xn,14,Nfft,fs);%调用最大熵函数pmem，滤波器阶数14 subplot(2,1,2), plot(f,10log10(Pxx)), title(' 最大熵法，滤波器14'),xlabel('频率'),ylabel('功率谱db');

（5）最大似然法： 思想：最大似然法原理是让信号通过一个滤波器，选择滤波器的参数使所关心的频率的正弦波信号能够不失真地通过，同时，使所有其他频率的正弦波通过这个滤波器后输出的均方值最小。在这个条件下，信号经过这个滤波器后输出的均方值就作为其最大似然法功率谱估值。可以证明,如果信号x 是由一个确定性信号S 加上一个高斯白噪声n 所组成，则上述滤波器的输出是信号S 的最大似然估值。如果n 不是高斯噪声，则上述滤波器的输出是信号S 的最小方差的线性的无偏估值，而且它能得到比使用固定的窗口函数的周期图法更高的分辨率。

程序： fs=20000; n=0:1/fs:1; N=length(n); W=2000pi; x1n=square(Wn); x2n=randn(1,N); xn=x1n+02x2n; figure(4) subplot(3,1,1) plot(n,xn); axis([0 001 -5 5]) %估计自相关函数 m=-500:500; [r,lag]=xcorr(xn,500,'biased'); R=[r(501) r(502) r(503) r(504); r(500) r(501) r(502) r(503); r(499) r(500) r(501) r(502); r(498) r(499) r(500) r(501)]; [V,D]=eig(R); V3=[V(1,3),V(2,3),V(3,3),V(4,3)]'; V3=[V(1,4),V(2,4),V(3,4),V(4,4)]'; p=0:3; wm=[0:0002pi:2pi]; B=[(exp(-li))^(wm'p)];%li 就是虚数单位j A=B; %最小方差功率谱估计 z=Ainv(R)A'; Z=diag(z'); pmv=1/Z; subplot(2,1,2) plot(wm/pi,pmv); title('基于最大似然的功率谱估计') ylabel('功率谱幅度（db）') xlabel('角度频率w/pi')

自己分下段就能运行了。

你运行了这个程序，发现报错了。

在Matlab里面，这个shang只是一个函数。这不是这样运行的。

在其他M文件里面是可以调用这个函数的，

比如写个：

先定义个x 。然后运行；x是输入

shang(x)

之后，你就得到了结果。输出的w

这就是matlab里面函数的概念。

你可以试试这个程序，这样：

先定义

a=magic（4）；

然后把第二行到结尾都复制进去，再试。就应该有结果了。

最大化就实现figure与屏幕大小一致，所以可以先获取屏幕大小，之后将figure的位置属性设置为获取到的屏幕大小。

具体实现方法可以参考如下程序段：

% figure 窗口最大化，坐标轴也随着窗口变大而相应变大

scrsz = get(0,'ScreenSize'); % 是为了获得屏幕大小，Screensize是一个4元素向量[left,bottom, width, height]

set(gcf,'Position',scrsz); % 用获得的screensize向量设置figure的position属性，实现最大化的目的

L(i,j)=mm(i)mm2(j); 这一条意味着两个变量相互独立，最终答案肯定为0，

m=[783 132 432 432 72 245 339 679 143 20 267 469 230 1598 128 2156 139 139 702 1296 ];

m = mapminmax(m,0,1);

m2=[7861 9698 8186 8291 8135 8395 7876 8572 8398 8447 8473 7614 12469 12469 12469 62347 12469 12469 1875 63002 ];

m2 = mapminmax(m2,0,1);

Y=hist3([m',m2'])/100;

mm=sum(Y,1);

mm2=sum(Y,2);

mm2=mm2';

%

% mm=hist(m)/10;

% mm

mm(mm==0) = [ ];

% mm

Hm = -sum(mmlog2(mm));

Hm

mm2(mm2==0) = [ ];

mm2

Hm2 = -sum(mm2log2(mm2));

Hm2

Y=reshape(Y,1,100);

Y(Y==0)=[];

Hab=-sum(sum(Ylog2(Y)));

disp('Hab信源联合熵为：');

Hab

%计算a和b的互信息

mi =( Hm+Hm2)-Hab;

mi