运用matlab进行随机信号的功率谱密度估计仿真【急求】
怎么都没人啊自己顶一下。。 查看原帖>>
clear all;
close all;
warning off all;
Fs = 1000;
nfft=1024;
idx = 0:round(nfft/2-1);
k = idx*Fs/nfft;
t = 0:1/Fs:1;
x1 = rand(1,1001);
[cor1 lag1] = xcorr(x1,'unbiased');
figure(1);
subplot(211),plot(lag1/Fs,cor1),title('(0,1)均匀分布自相关函数');
Xk1 = fft(cor1,nfft);
Px1 = abs(Xk1);
subplot(212),plot(k,10*log10(Px1(idx+1))),title('(0,1)均匀分布功率谱密度');
x2 = normrnd(2,5,1,1001);
[cor2 lag2] = xcorr(x2,'unbiased');
figure(2);
subplot(211),plot(lag2/Fs,cor2),title('(2,5)正态分布自相关函数');
Xk2 = fft(cor2,nfft);
Px2 = abs(Xk2);
subplot(212),plot(k,10*log10(Px2(idx+1))),title('(2,5)正态分布功率谱密度');
x3 = cos(600*pi*t)+cos(640*pi*t)+randn(1,1001);
[cor3 lag3] = xcorr(x3,'unbiased');
figure(3);
subplot(211),plot(lag3/Fs,cor3),title('随机信号自相关函数');
Xk3 = fft(cor3,nfft);
Px3 = abs(Xk3);
subplot(212),plot(k,10*log10(Px3(idx+1))),title('随机信号功率谱密度');
fai = 2*pi*rand(1,1001);
x4 = 2*cos(1000*pi*t+fai);
[cor4 lag4] = xcorr(x4,'unbiased');
figure(4);
subplot(211),plot(lag4/Fs,cor4),title('随机相位信号自相关函数');
Xk4 = fft(cor4,nfft);
Px4 = abs(Xk4);
subplot(212),plot(k,10*log10(Px4(idx+1))),title('随机相位信号功率谱密度');
这份是我在网络上找的：
（1）.周期图法： 思想：周期图法是为了得到功率谱估值，先取信号序列的离散傅里叶变换，然后取其幅频特性的平方并除以序列长度 N。由于序列 x(n)的离散傅里叶变换 X(k) 具有周期性，因而这种功率谱也具有周期性，常称为周期图。周期图是信号功率谱的一个有偏估值；而且，当信号序列的长度增大到无穷时，估值的方差不趋于零。因此，随着所取的信号序列长度的不同，所得到的周期图也不同，这种现象称为随机起伏。由于随机起伏大，使用周期图不能得到比较稳定的估值。
程序：
%首先，生成输入信号的程序为：
clear ; fs=20000; n=0:1/ fs: 0. 1; N=lengt h( n); W=2000*pi;%因方波频率 F=1000HZ 所以角频率 W=2000pi X1n=square( W*n) ;%方波信号 X2n=randn( 1, N);%白噪声信号 xn=X1n+0. 2*X2n; subplot (3, 1, 1) plot (n, xn) ; xlabel( ' n' ) ylabel('输入信号') %其次，开始用周期图法进行估计； clear all fs=20000; n=0:1/fs:0.1; N=length(n); W=2000*pi; x1n=square(W*n); x2n=randn(1,N); xn=x1n+0.2*x2n; 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,10*log10(Pxx)),%转成DB 单位
（2）.自相关函数法： 思想： 随机信号 x(n)的相关函数是在时间域内描述随机过程的重要特征。自相关函数是随机信号在不同时刻的值之间的依赖性的量度，是一个很有用的统计平均量，其定义为自相关函数 (1) 式中E【·】表示数学期望，*表示共轭值,m 为时间滞后数。在随机信号处理中，自相关函数可以用来检测淹没在随机噪声干扰中的信号，随机信号的自功率谱等于它的自相关函数的傅里叶变换。因此，通过自相关估计可求得信号的功率谱。利用计算机计算自相关估值有两种方法。一种是直接方法，先计算出随机信号和它的滞后序列的乘积，再取其平均值即得相关函数的估计值。另一种是间接方法，先用快速变换算法计算随机序列的功率谱密度，再作反变换计算出相关函数。
程序： n=0:1/fs:1; N=length(n); W=2000*pi; x1n=square(W*n); x2n=randn(1,N); xn=x1n+0.2*x2n; figure(2) subplot(3,1,1) plot(n,xn);%输入信号 axis([0 0.01 -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(-li*pi/500).^(m'*k));%对R 求傅里叶变换,li 就是j magX=abs(X); subplot(3,1,3) plot(M,10*log10(magX)); xlabel('功率谱的改进直接法估计')
（3）.自协方差法： 思想：在实际中,有时用随机信号在两个不同时刻t 1 ,t 2 的取值X(t 1 )和X(t 2 )之间的二阶混合中心矩来描述随机信号 X(t)在任意两个时刻取值起伏变化的相依程度, 即自协方差。自协方差函数与自相关函数描述的特性是一致的，所以其原理与相关函数法近似。
clear all; fs=20000; n=0:1/fs:0.1; P=2000*pi; y=square(P*n); xn=y+0.2*randn(size(n));%绘制信号波形 figure(3) subplot(2,1,1) plot(n,xn) xlabel('时间（s）') ylabel('幅度') title('y+randn(size(n))') ymax_xn=max(xn)+0.2; ymin_xn=min(xn)-0.2; axis([0 0.3 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=10*log10(xpsd+0.000001); 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=2000*pi; x1n=square(W*n); x2n=randn(1,N); xn=x1n+0.2*x2n; figure(5) subplot(3,1,1) plot(n,xn); asis([0 0.01 -5 5]) Nfft=256;%分段长度256 [Pxx,f]=pmem(xn,14,Nfft,fs);%调用最大熵函数pmem，滤波器阶数14 subplot(2,1,2), plot(f,10*log10(Pxx)), title(' 最大熵法，滤波器14'),xlabel('频率'),ylabel('功率谱db');
（5）.最大似然法： 思想：最大似然法原理是让信号通过一个滤波器，选择滤波器的参数使所关心的频率的正弦波信号能够不失真地通过，同时，使所有其他频率的正弦波通过这个滤波器后输出的均方值最小。在这个条件下，信号经过这个滤波器后输出的均方值就作为其最大似然法功率谱估值。可以证明,如果信号x 是由一个确定性信号S 加上一个高斯白噪声n 所组成，则上述滤波器的输出是信号S 的最大似然估值。如果n 不是高斯噪声，则上述滤波器的输出是信号S 的最小方差的线性的无偏估值，而且它能得到比使用固定的窗口函数的周期图法更高的分辨率。
程序： fs=20000; n=0:1/fs:1; N=length(n); W=2000*pi; x1n=square(W*n); x2n=randn(1,N); xn=x1n+0.2*x2n; figure(4) subplot(3,1,1) plot(n,xn); axis([0 0.01 -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:0.002*pi:2*pi]; B=[(exp(-li)).^(wm'*p)];%li 就是虚数单位j A=B; %最小方差功率谱估计 z=A*inv(R)*A'; Z=diag(z'); pmv=1./Z; subplot(2,1,2) plot(wm/pi,pmv); title('基于最大似然的功率谱估计') ylabel('功率谱幅度（db）') xlabel('角度频率w/pi')
自己分下段就能运行了。
楼主牛B
Matlab如何进行换行?
一、明确答案 在Matlab中换行，可以通过以下几种方式实现：1. 使用回车键进行换行。2. 在编写代码时，使用三个连续的单引号开始多行字符串，并在每行的末尾自动换行。3. 使用MATLAB的编程环境特有的换行符或快捷键Ctrl+Enter或Command+Return进行换行。二、详细解释 1. 使用回车键换行 在Matlab的命令...
怎么用matlab设置x和y坐标轴的值?
可以用matlab提供的坐标轴设置函数axis进行设置（对当前坐标轴进行设置），具体用法如下：1、plot(x, y); \/\/ 画图后用axis函数设置坐标轴的范围。2、axis([xmin xmax ymin ymax]); % 设置坐标轴在指定的区间。3、xmin、xmax 表示设置横坐标的最小最大值。4、ymin、ymax 表示设置纵坐标的最小...
如何有效地利用Matlab进行纯数学研究?
1.学习和掌握Matlab的基本语法和函数：这是使用Matlab进行数学研究的基础。你需要了解如何使用Matlab编写脚本，如何调用内置函数，以及如何创建自定义函数。2.利用Matlab的数值计算功能：Matlab提供了丰富的数值计算工具，如矩阵运算、微积分、线性代数等。你可以利用这些工具进行复杂的数学计算和分析。3.使用Matl...
matlab如何做线性拟合
首先，最常用的线性拟合可以通过内置的多项式拟合功能实现。在命令窗口中，输入自变量x和因变量y，使用polyfit函数进行操作。例如，若要进行二次多项式拟合，只需键入p = polyfit(x,y,2)。若需拟合更高阶的多项式，只需调整括号中的数字即可。另一种方法是利用MATLAB自带的曲线拟合工具。在APPS选项中，找...
怎么用matlab仿真啊?
如何用MATLAB进行仿真？MATLAB是一款强大的数学计算与仿真软件，广泛应用于各个领域。在MATLAB中进行仿真，通常涉及以下几个步骤：一、明确仿真目的 首先，需要明确仿真的目的。这涉及到确定要模拟的现实系统或现象，以及希望通过仿真得到哪些数据或结果。明确目的有助于设定仿真模型的参数和范围。二、建立仿真...
matlab怎么运行代码
在命令窗口中，您可以输入任何Matlab代码，并按“Enter”键运行它。例如，以下代码将打印出“Hello, world！”：disp(\\'Hello, world!\\')可以在Matlab的命令窗口中执行单个命令，也可以编写一个脚本来执行多个命令。要运行脚本，请使用“run”命令，后跟文件名，如下所示：run myscript.m 编辑器运行...
新手如何利用matlab软件进行简单的矩阵运算
首先我们先打开matlab这款软件。提示一下：不同版本的计算方式都是一样的。在Fx的后面输入：A=[1 2 3;4 5 6;7 8 9]这些字符。其中分号是表示换行，对于这个矩阵，也就是有三行三列。我们按一下回车键，看会出现什么结果。这样一个A矩阵就输入到matlab当中了 接下来是利用同样的方法再输入一个...
如何使用MATLAB进行编程以及绘图?
用MATLAB进行编程以及绘图很简单，具体操作步骤如下：1、找到matlab软件，并打开它。2、我们需要新建一个界面并在其中书写程序。首先我们在左上角找到File这一选项（点击），找到new这一选项鼠标移到该位置并找到script选项（点击）。就会弹出编写程序的界面了。3、我们进入我们创建的matlab的全新界面后，...
MATLAB中如何进行换行?
按下Enter键即可实现换行。这种换行不涉及MATLAB编程语言的特性，仅仅是文本编辑的基本操作。总的来说，在MATLAB中换行主要依赖于回车键和分号的使用，以及特定场合下三个引号的使用来实现多行字符串的处理。编程时可以根据需要和代码的上下文选择合适的方式进行换行，以提高代码的可读性和效率。
matlab怎么求反函数
有时候我们在使用matlab进行数学运算的时候，想求反函数，怎么求呢，下面来分享一下方法 第一步我们首先需要知道在matlab中求反函数用到的是finverse函数，在命令行窗口中输入“helpfinverse”，可以看到函数的使用方法，第二步g=finverse(f)格式，f符号函数表达式，变量x，求得的反函数g是满足g(f(x))...
昌善科米： 关于谱能量,有这样一种解释,你可以试着去算一算 信号可以分成能量信号与功率信号,非周期能量信号具有能量谱密度,是傅立叶变换的平方,功率信号具有功率谱密度,其与自相关函数是一对傅立叶变换对,等于傅立叶变换的平方/区间长...
叶县15865031093： 用matlab求电流的功率谱!? - ？
昌善科米： 设信号为x,则求功率谱的MATLAB命令是:Pw=fft(conv(x,fliplr(x))), 即利用功率谱等于自相关系数的傅里叶变换来计算.
叶县15865031093： matlab中如何生成一个随机信号 - ？
昌善科米： 原发布者:李世民的111Matlab中随机信号的产生在matlab编程中,我们所能用到的用于产生随机信号的函数有三:Rand,randn,randi下面我们详细的了解一下这三个函数.1.Rand功能是生产均匀分布的伪随机数,并且所生成的伪随机数分布在...
叶县15865031093： 请问matlab simulink中怎么测量信号功率啊,有直接测试的模块吗?我找了好久都没找到啊~跪求高人指点？
昌善科米： 印象中是没有这样的模块的.估计对于信号处理,功率并不重要.
叶县15865031093： 明天老师要收的作业,求大神救命,随机信号处理,用MATLAB编程 - ？
昌善科米： a=[1 -1.6 1.28]; b=[0.5 0.1]; n=64; pulse=[1,zeros(1,n-1)]; %单位脉冲序列 %1)单位冲激响应 h = filter(b,a,pulse); %或h1 = impz(b,a); figure(1); stem(h); title('单位冲激响应(filter函数)'); %2)输入响应 x=zeros(1,n-1); x(1)=1; x(2)=2; x(3)=3; ...
叶县15865031093： 求用matlab表示下面信号的功率谱的程序:x(n)=10sin(2π*0.2n+π/3)+ - ？
昌善科米： Matlab代码示例:clear; Fs=1000; %采样频率 n=0:1/Fs:1;%产生含有噪声的序列 xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n)); window=boxcar(length(xn)); %矩形窗 nfft=1024; [Pxx,f]=periodogram(xn,window,nfft,Fs); %直接法 plot(f,10*log10(Pxx)); 带入你自己的信号就行了
叶县15865031093： 在matlab中混入噪声功率为4W的随机噪声的代码是什么 - ？
昌善科米： M=imread('dl011.jpg') %读取MATLAB中的名为cameraman的图像 subplot(3,3,1) imshow(M) %显示原始图像 title('original') P1=imnoise(M,'gaussian',0.02) %加入高斯躁声 subplot(3,3,2) imshow(P1) %加入高斯躁声后显示图像 title('gaussian noise'...
叶县15865031093： 怎样用matlab产生宽带随机信号 - ？
昌善科米： Matlab 伪随机信号的产生(M序列)这里主要介绍两种方法,一是使用m文件编辑,二是用simulink产生 一、用m文件编辑 Np=63;%循环周期 delta_T = 1;%时钟节拍 a=1;%幅度 %初始化M序列 M(1)=1; M(2)=0; M(3)=0; M(4)=1; M(5)=1; M(6)...
叶县15865031093： 功率的求法 - ？
昌善科米： 力学中的功:P=W/t P=F*v P=ma/v 电学中的功:P=UI P=RI^2 P=U^2/R
叶县15865031093： 如何用MATLAB绘制功率谱密度图形 - ？
昌善科米： 功率谱:信号先自相关再作FFT 频谱:信号直接作FFT. 幅度谱的平方(二次量纲)又叫能量谱,它描述了信号能量的频域分布; 功率信号的功率谱描述了信号功率随频率的分布特点; 相频体现了图像的总体信息(低)和细节信息(高);
