【随机信号分析2】随机信号通过线性系统仿真和窄带随机信号仿真

1.讨论随机信号通过线性系统之后的自相关函数和功率谱密度。
2.窄带随机信号的产生与窄带高斯信号包络和相位的概率分布。

1、设计一个带宽为 0.3 的 FIR 带通数字滤波器，其低端截止频率为 0.25，并画出滤波器的 幅频特性和相频特性。

实验程序：
h=fir1(71,[0.25 0.55]);
[HH,WW]=freqz(h,1,512);
subplot(211);
plot(WW/pi,20*log10(abs(HH)));
title(“幅度谱”)
grid on;
axis([0 1 -100 0]);
subplot(212);
plot(WW/pi,unwrap(angle(HH))*180/pi);
title(“相位谱”)
grid on;

2、平均功率为 1 的高斯白噪声信号通过一个带通系统，其系统的两个截止频率分别为 2kHz 和 6kHz，求输出的自相关函数和功率谱密度。

实验程序：
N=500;
xt=random(‘norm’,0,1,1,N);
ht=fir1(101,[0.2,0.6]);
HW=fft(ht,2N);
Rxx=xcorr(xt,‘biased’); Sxx=abs(fft(xt,2N).^2)/(2*N); HW2=abs(HW).^2;
Syy=Sxx.HW2;
Ryy=fftshift(ifft(Syy));
w=(1:N)/N;
t=(-N:N-1)/N(N/20000);
subplot(4,1,1);
plot(w,abs(Sxx(1:N)));
title(“输入信号功率谱密度”)
subplot(4,1,2);
plot(w,abs(HW2(1:N)));
title(“系统功率传输函数”)
subplot(4,1,3);
plot(w,abs(Syy(1:N)));
title(“输出信号功率谱密度”)
subplot(4,1,4);
plot(t,Ryy);
title(“输出信号自相关函数”)

3、

（1）产生中心频率为 f0为10kHz、带宽为∆f为500Hz 的窄带高斯随机过程 X (t)的样本，
（2）产生其同相分量AC(t)和正交分量AS(t)的样本，
（3）对X(t)、AC(t)和AS(t)进行功率谱估计。

实验程序：
function X=Narrowbandsignal(N,f0,deltf,fs,M)
N1=N-M;
xt=random(‘norm’,0,1,[1,N1]);
f1=f0*2/fs;
df1=deltf/fs;
ht=fir1(M,[f1-df1 f1+df1]);
X=conv(xt,ht);
return

function [Ac As]=Lowfsignal(X,f0,fs)
HX=imag(hilbert(X));
[M N]=size(X);
t=0:1/fs:((N-1)/fs);
Ac=X.cos(2pif0t)+HX.sin(2pif0t);
As=HX.cos(2pif0t)-X.sin(2pif0t);
return

N=10000;f0=10000;deltf=500;fs=22000;M=200;
X=Narrowbandsignal(N,f0,deltf,fs,M);
[Ac As]=Lowfsignal(X,f0,fs);
Rx=xcorr(X,‘biased’);
Rac=xcorr(Ac,‘biased’);
Ras=xcorr(As,‘biased’);
Racw=abs(fft(Ras));
Rasw=abs(fft(Ras));
Rxw=abs(fft(Rx));
N1=2N-1;
f=fs/N1:fs/N1:fs/2;
subplot(3,1,1); plot(f,10log10(Rxw(1:(N1-1)/2)+eps));
title(“X(t)的功率谱密度”);
subplot(3,1,2); plot(f,10log10(Racw(1:(N1-1)/2)+eps));
title(“Ac(t)的功率谱密度”);
subplot(3,1,3);
plot(f,10log10(Rasw(1:(N1-1)/2)+eps)); title(“As(t)的功率谱密度”);

4、分别产生中心频率为f0为10kHz、带宽为∆f为500Hz、方差为1的窄带高斯随机过程X (t)，及其包络A(t)和相位Φ(t)、包络平方A2(t)的样本，并对它们的概率分布进行估计。

实验程序：
function X=Narrowbandsignal(N,f0,deltf,fs,M)
N1=N-M;
xt=random(‘norm’,0,1,[1,N1]);
fl=f0*2/fs; dfl=deltf/fs;
ht = fir1(M,[fl-dfl fl+dfl]);
X=conv(xt,ht);
return

function [At Ph A2]=EnvelopPhase(X,f0,fs)
HX=imag(hilbert(X));
[M N]=size(X);
t=0:1/fs:((N-1)/fs);
Ac=X.cos(2pif0t)+HX.sin(2pif0t);
As=HX.cos(2pif0t)-X.sin(2pif0t);
Ph=atan(As./Ac);
A2=Ac.*Ac+As.*As;
At=sqrt(A2);

N=20000;f0=10000;
deltf=500;
fs=22000;
M=50;
X=Narrowbandsignal(N,f0,deltf,fs, M) ;
X=X/sqrt(var(X) ) ;
[At,Ph ,A2]=EnvelopPhase(X, f0, fs) ;
LA=0:0.05:4.5;
hist(At, LA) ;
title(“包络A(t)样本值的分布直方图”);
LP=-pi/2:0.05:pi/2;
figure;
hist(Ph, LP) ;
title(“相位Φ(t)样本值的分布直方图”);
LA2=0:0.2:16;
figure;
hist(A2, LA2) ;
title(“包络平方A2(t)样本值的分布直方图”);

5、仿真中心频率为f0为10kHz、带宽为∆f为400Hz、方差为1的窄带高斯随机过程X (t)，在下面三种余弦信号情况下：

（1）余弦信号的幅度为2，相位为π/6；
（2）余弦信号的幅度为 4，相位为π/4；
（1）余弦信号的幅度为 8，相位为π/3；
分析窄带高斯随机过程加余弦信号的合成信号的包络、相位、包络平方的概率分布。

实验程序：
function [At,Ph,A2]=EnvelopPhase(X,f0,fs)
HX=imag(hilbert(X));
[M,N]=size(X);
t=0:1/fs:((N-1)/fs);
Ac=X.cos(2pif0t)+HX.sin(2pif0t);
As=HX.cos(2pif0t)-X.sin(2pif0t);
Ph=atan(As./Ac);
A2=Ac.*Ac+As.*As;
At=sqrt(A2);

function X=Narrowbandsignal(N,f0,delt,fs,M);
N1=N-M;
xt=random(‘norm’,0,1,[1,N1]);
f1=f0*2/fs;
df=delt/fs;
ht=fir1(M,[f1-df f1+df]);
X=conv(xt,ht)
return

N=10000;f0=10000;delt=400;fs=22000;M=50;
al=2;a2=4;a3=8;
sitl=pi/6;sit2=pi/4;sit3=pi/3
X=Narrowbandsignal(N,f0,delt,fs,M);
X=X/sqrt(var(X));
t=0:1/fs:((N-1)/fs); X1=X+alcos(2pif0t+sitl);
X2=X+a2cos(2pif0t+sit2);
X3=X+a3cos(2pif0t+sit3);
[Atl,Ph1,A21]=EnvelopPhase(X1, f0, fs);
[At2,Ph2,A22]=EnvelopPhase(X2, f0, fs);
[At3,Ph3,A23]=EnvelopPhase(X3, f0, fs);
LA=0:0.4:12;
GA1=hist(Atl, LA);
GA2=hist(At2, LA);
GA3=hist(At3, LA);
plot(LA,GA1,‘:’,LA,GA2,‘-’,LA,GA3,‘–’);
title(“包络的分布直方图”);
figure;
LP=-pi/2:0.05:pi/2;
GP1=hist((Ph1-sitl),LP);
GP2=hist((Ph2-sit2),LP);
GP3=hist((Ph3-sit3),LP);
plot(LP,GP1,‘:’,LP,GP2,‘-’,LP,GP3,‘–’);
title(“相位的分布直方图”);
figure;
LA2=0:1:120;
GA21=hist(A21,LA2);
GA22=hist(A22,LA2);
GA23=hist(A23,LA2);
plot(LA2,GA21,‘:’, LA2,GA22,‘-’,LA2,GA23,‘–’);
title(“包络平方值的分布直方图”);

选做题：

2、模拟产生一个窄带随机过程，首先产生两个互相独立的随机过程AC(t)和AS(t)，并将用两个正交载波cos(2πf0t)和sin(2πf0t)进行调制，如下图所示，然后进行抽样得到窄带随机过程的抽样。

实验程序：
N=1000;f0=1000/pi;
T=1;
x1=randn(N,1);
x2=randn(N,1);
b=[1 -0.9]
Ac=filter(b,1,x1);
As=filter(b,1,x2);
t=0:T;(N-1)/T;
S=Ac.cos(2pif0t)- As.sin(2pif0t);
c = xcorr(S);
pxx = periodogram(S)
figure(1);
plot©;
title(“自相关”);
figure(2);
plot(pxx)
title(“功率谱”);