基于MATLAB的随机信号及其参数建模分析.docx

基于MATLAB的随机信号及其参数建模分析.docx

《基于MATLAB的随机信号及其参数建模分析.docx》由会员分享，可在线阅读，更多相关《基于MATLAB的随机信号及其参数建模分析.docx（29页珍藏版）》请在冰豆网上搜索。

基于MATLAB的随机信号及其参数建模分析

随机信号基础

Binornd函数示例：

>>R=binornd（10,0.5）

R=

5

>>R=binornd（10,0.5,1,6）

R=

357482

>>R=binornd（10,0.5,[1,10]）

R=

5835753652

>>R=binornd（10,0.5,[2,3]）

R=

497

565

>>n=10:

10:

60;

>>r1=binornd（n,1./n）

r1=

000211

>>r2=binornd（n,1./n,[16]）

r2=

130120

Normrnd函数示例：

>>n1=normrnd（1:

6,1./（1:

6））

n1=

0.56741.16723.04184.07194.77076.1985

>>n2=normrnd（0,1,[1,5]）

n2=

1.1892-0.03760.32730.1746-0.1867

>>n3=normrnd（[123;456],0.1,2,3）

n3=

1.07262.21833.0114

3.94124.98646.1067

>>R=normrnd（10,0.5,[23]）

R=

10.02969.58389.3319

9.952210.147210.3572

产生15（3行5列）个均值为2、标准差为5的正态分布随机数：

>>y=random（'norm',2,0.5,3,5）

y=

2.81182.62702.28562.40782.3343

1.65411.20311.80012.35602.5954

2.42901.27952.34502.64511.3988

计算正态分布N（0，1）的随机变量X在点0.6578的密度函数值：

>>pdf（'norm',0.6578,0,1）

ans=

0.3213

计算自由度为8卡方分布，在点2.18处的密度函数值：

>>pdf（'chi2',2.18,8）

ans=

0.0363

绘制卡方分布密度函数在自由度分别为2、6、16时的图像：

>>x=0:

0.1:

30;

>>y1=chi2pdf（x,2）;

>>plot（x,y1,'-.'）;

>>holdon;

>>y2=chi2pdf（x,6）;

>>plot（x,y2,'*'）;

>>y3=chi2pdf（x,16）;

>>plot（x,y3）;

>>axis（[03000.2]）;

二项分布示例：

>>x=0:

10;

>>y=binopdf（x,10,0.5）;

>>plot（x,y,'o'）;

卡方分布示例：

>>x=0:

0.2:

15;

>>y=chi2pdf（x,4）;

>>plot（x,y,'+'）;

非中心卡方分布图示例：

>>x=（0:

0.1:

10）';

>>y1=ncx2pdf（x,4,2）;

>>y=chi2pdf（x,4）;

>>plot（x,y,':

',x,y1,'-'）;

指数分布示例：

>>x=0:

0.1:

10;

>>y=exppdf（x,2）;

>>plot（x,y,'--'）;

F分布示例：

>>x=0:

0.1:

10;

>>y=fpdf（x,5,3）;

>>plot（x,y,'-'）;

非中心F分布：

>>x=（0.01:

0.1:

10.01）';

>>y1=ncfpdf（x,5,20,10）;

>>y=fpdf（x,5,20）;

>>plot（x,y,':

',x,y1,'-'）;

『分布：

>>x=gaminv（（0.005:

0.01:

0.995）,100,10）;

>>y=gampdf（x,100,10）;

>>y1=normpdf（x,1000,100）;

>>plot（x,y,'--',x,y1,':

'）;

正态分布示例：

>>x=-3:

0.2:

3;

>>y=normpdf（x,0.1）;

>>plot（x,y）;

负二项分布示例：

>>x=（0:

10）;

>>y=nbinpdf（x,3,0.5）;

>>plot（x,y,'*'）;

对数正态分布示例：

>>x=（10:

1000:

125010）';

>>y=lognpdf（x,log（20000）,1.0）;

>>plot（x,y,'-.'）;

瑞利分布：

>>x=[0:

0.01:

2];

>>y=raylpdf（x,0.5）;

>>plot（x,y,'-.'）;

泊松分布示例：

>>x=0:

15;

>>y=poisspdf（x,5）;

>>plot（x,y）;

威布尔分布示例：

>>t=0:

0.1:

3;

>>y=weibpdf（t,2,2）;

>>plot（y）

T分布示例：

>>x=-5:

0.1:

5;

>>y=tpdf（x,5）;

>>z=normpdf（x,0,1）;

>>plot（x,y,':

',x,z,'-'）;

随机信号的数字特征

计算方差N=100000的正态分布高斯随机信号的平均值、平方值、均方差、方差：

>>N=100000;

randn（'state',0）;

y=randn（1,N）;

>>disp（'平均值:

'）;

平均值:

>>yMean=mean（y）

yMean=

0.0032

>>disp（'平方值:

'）;

平方值:

>>y2p=y*y'/N

y2p=

1.0073

>>disp（'均方差'）;

均方差

>>ystd=std（y,1）

ystd=

1.0037

>>disp（'方差'）;

方差

>>yd=ystd.*ystd

yd=

1.0073

随机信号的自相关与协方差

编写程序实现白噪声、正弦信号和带有白噪声的正弦信号的自相关函数：

>>fs=1000;

>>t=0:

1/fs:

（1-1/fs）;

>>maxlag=100;

>>x=randn（1,1000）;

>>[c,maxlags]=xcorr（x,maxlag）;

>>figure

（1）;

>>subplot（2,2,1）;plot（t,x）;

>>xlabel（'t'）;ylabel（'x（t）'）;

>>title（'白噪声'）;

>>subplot（2,2,2）;plot（maxlags/fs,c）;

>>xlabel（'x'）;ylabel（'Rx（t）'）;

>>title（'白噪声自相关'）;

>>y=sin（2*pi*20*t）+0.8*randn（1,1000）;

>>[c,lags]=xcorr（y,maxlag）;

>>subplot（2,2,3）;plot（t,y）;

>>xlabel（'x'）;ylabel（'y（t）'）;

>>title（'原始信号'）;

>>subplot（2,2,4）;plot（lags/fs,c）;

>>xlabel（'t'）;ylabel（'Ry（t）'）;

>>title（'自相关'）;

比较一下某函数的互相关函数和互协方差函数：

>>t=1:

20;x=t.^2;

>>y=（t+2）.^2;

>>R=xcorr（x,y）;

>>c=xcov（x,y）;

>>n=1:

length（c）;

>>subplot（2,1,1）;stem（n,c）;

>>title（'互协方差'）;

>>subplot（2,1,2）;stem（n,R）;

>>title（'互相关'）;

功率谱估计

1、最大熵法

用最大熵法估计含有噪声和50Hz、120Hz周期信号的功率谱密度，并与Welch法得到的结果进行比较：

>>Fs=1000;

>>N=1024;Nfft=256;

>>n=0:

N-1;t=n/Fs;

>>window=hanning（256）;

>>randn（'state',0）;

>>xn=sin（2*pi*50*t）+2*sin（2*pi*120*t）+randn（1,N）;

>>[Pxx1,f]=pmem（xn,14,Nfft,Fs）;

>>subplot（2,1,1）;plot（f,10*log10（Pxx1））;

>>xlabel（'频率/Hz'）;ylabel（'功率谱/dB'）;

>>title（'最大熵法（MTM）Order=14'）;

>>gridon;

>>noverlap=128;

>>dflag='none';

>>subplot（2,1,2）;

>>psd（xn,Nfft,Fs,window,noverlap,dflag）;

>>title（'Welch方法估计的功率谱'）;

>>gridon;

2、多窗口法

假定NW=4和NW=2，使用多窗口法估计含有噪声和50Hz、120Hz周期信号的功率谱密度：

>>Fs=1000;

N=1024;Nfft=256;

n=0:

N-1;t=n/Fs;

randn（'state',0）;

xn=sin（2*pi*50*t）+2*sin（2*pi*120*t）+randn（1,N）;

[Pxx1,f]=pmtm（xn,4,Nfft,Fs）;

>>subplot（2,1,1）;plot（f,10*log10（Pxx1））;

xlabel（'频率/Hz'）;ylabel（'功率谱/dB'）;

>>title（'多窗口法（MTM）NW=4'）;

>>gridon;

>>[Pxx2,f]=pmtm（xn,2,Nfft,Fs）;

>>subplot（2,1,2）;plot（f,10*log10（Pxx2））;

>>xlabel（'频率/Hz'）;ylabel（'功率谱/dB'）;

title（'多窗口法（MTM）NW=2'）;

>>gridon;

3、信号分类法

使用多信号分类，采用7个窗口，估计含有噪声和50Hz、120Hz周期信号的功率谱密度：

>>Fs=1000;

>>N=1024;Nfft=256;

>>n=0:

N-1;t=n/Fs;

>>randn（'state',0）;

>>xn=sin（2*pi*50*t）+2*sin（2*pi*120*t）+randn（1,N）;

>>pmusic（xn,[7,1.1],Nfft,Fs,32,16）;

>>xlabel（'频率/Hz'）;ylabel（'功率谱/dB'）;

>>title（'通过MUSIC法估计的伪谱'）;

>>gridon;

倒谱分析

设原信号是一个45Hz的正弦波，在传输过程中遇到障碍产生回声，回声振幅误差为原信号的0.5，并与原信号有0.2s的延迟。

在某测点测到的信号是原信号和回声信号的叠加：

>>Fs=100;

>>t=0:

1/Fs:

1.27;

>>s1=sin（2*pi*45*t）;

>>s2=s1+0.5*[zeros（1,20）s1（1:

108）];

>>c=cceps（s2）;

>>plot（t,c）;

>>gridon;

时域建模

使用一个3阶前向预测器来估计模型系数，并用AR输出最后4096点数据：

>>randn（'state',0）;

>>noise=randn（50000,1）;

>>x=filter（1,[11/21/31/4],noise）;

>>x=x（45904:

50000）;

>>a=lpc（x,3）;

>>est_x=filter（[0-a（2:

end）],1,x）;

>>e=x-est_x;

>>[acs,lags]=xcorr（e,'coeff'）;

>>plot（1:

97,x（4001:

4097）,1:

97,est_x（4001:

4097）,'--'）;

>>gridon;

>>plot（lags,acs）;

建立一个4阶巴特沃思原始滤波器，并从脉冲响应巴特沃思过滤器中恢复出它的系数：

>>[b,a]=butter（4,0.2）

b=

0.00480.01930.02890.01930.0048

a=

1.0000-2.36952.3140-1.05470.1874

>>h=filter（b,a,[1zeros（1,25）]）;

>>[bb,aa]=prony（h,4,4）

bb=

0.00480.01930.02890.01930.0048

aa=

1.0000-2.36952.3140-1.05470.1874

建立一个6阶的巴特沃思脉冲响应滤波器，并利用stmcb函数对该滤波器进行预测：

>>[b,a]=butter（6,0.2）;

>>h=filter（b,a,[1zeros（1,100）]）;

>>freqz（b,a,128）;

>>[bb,aa]=stmcb（h,4,4）;

>>freqz（bb,aa,128）;

>>x=[1111];

>>y=filter（1,[11],x）;

>>[b,a]=stmcb（y,x,0,1）

b=

1.0000

a=

1.00001.0000

比较函数lpc、prony及stmcb的效果：

>>randn（'state',0）;

>>aa=[10.10.10.10.1];

>>bb=1;

>>x=impz（bb,aa,10）+randn（10,1）/10;

>>a1=lpc（x,3）;

>>[b2,a2]=prony（x,3,3）;

>>[b3,a3]=stmcb（x,3,3）;

>>[x-impz（1,a1,10）x-impz（b2,a2,10）x-impz（b3,a3,10）]

ans=

-0.04330-0.0000

-0.024000.0234

-0.00400-0.0778

-0.0448-0.00000.0498

-0.2130-0.1303-0.0742

0.15450.16550.1270

0.14260.14040.1055

0.00680.01880.0465

0.03290.03760.0530

0.01080.0081-0.0162

>>sum（ans.^2）

ans=

0.09530.06590.0471

频域建模

模拟滤波器的频率建模

从一个简单的数字滤波器的频率响应中恢复出原始滤波器的系数：

>>a=[123214];b=[12323];

>>[h,w]=freqs（b,a,64）;

>>[bb,aa]=invfreqs（h,w,4,5）

bb=

1.00002.00003.00002.00003.0000

aa=

1.00002.00003.00002.00001.00004.0000

>>[bbb,aaa]=invfreqs（h,w,4,5,[],30）

bbb=

0.68162.10152.66940.9113-0.1218

aaa=

1.00003.46767.40606.21022.54130.0001

数字滤波器的频率建模

>>[b,a]=butter（4,0.4）;

>>[h,w]=freqz（b,a,64）;

>>wt=ones（size（w））;

>>[bbb,aaa]=invfreqz（h,w,3,3,wt,30）

bbb=

0.04640.18290.25720.1549

aaa=

1.0000-0.86640.6630-0.1614

数字滤波器建模与模拟滤波器建模的比较：

>>[b,a]=butter（4,0.4）;

>>[h,w]=freqz（b,a,64）;

>>[bb,aa]=invfreqz（h,w,3,3）;

>>wt=ones（size（w））;

>>[bbb,aaa]=invfreqz（h,w,3,3,wt,30）;

>>[H1,W1]=freqz（b,a）;

>>[H2,W2]=freqz（bb,aa）;

>>[H3,W3]=freqz（bbb,aaa）;

>>plot（W1/pi,abs（H1）,W2/pi,abs（H2）,':

',W3/pi,abs（H3）,'-'）;

>>gridon;

>>fvtool（b,a,bb,aa,bbb,aaa）;

>>sum（abs（h-freqz（bb,aa,w））.^2）

ans=

0.0200

>>sum（abs（h-freqz（bbb,aaa,w））.^2）

ans=

0.0096