随机信号分析试验(4)

东北石油大学计算机与信息技术学院李勇勇随机信号分析试验（4）1.编写程序,对均值为0方差为1的平稳高斯随机过程样本进行希尔伯特变换，绘制变换前后的概率密度直方图。注：z=hilbert(x)返回解析信号x+jH(x),所以imag(z)为x的希尔伯特变换clc;clearall;closeall;N=20000;x=random('normal',0,1,1,N);xt=imag(hilbert(x));n=-4:0.1:4;subplot(2,1,1)hist(x,n)subplot(2,1,2)hist(xt,n)2.编写Matlab函数,产生一个采样频率为fs中心频率为f0带宽为的窄带随机信号样本。注：函数Narrowbandsignal(N,f0,deltf,fs,M)产生窄带随机过程样本。N:长度,f0:单边功率谱中心频率,delta:带宽,fs:采样频率,M:产生窄带信号的滤波器阶数。(MN)functionX=Narrowbandsignal(N,f0,deltf,fs,M)N1=N-M;xt=random('normal',0,1,1,N1);f1=f0*2/fs;df1=deltf/fs;ht=fir1(M,[f1-df1,f1+df1]);X=conv(xt,ht);t=(1:N)/fs;plot(t,X)xlabel('时间(s)');ylabel('窄带随机信号样本X');运行结果：x=Narrowbandsignal(1000,2,0.5,10,20)3.编写Matlab函数,求采样频率为fs中心频率为f0的窄带随机过程X(t)的低频过程Ac(t)和As(t)的样本。注：函数Lowfsignal(X,f0,fs)产生Ac和As样本function[Ac,As]=Lowfsignal(X,f0,fs)HX=imag(hilbert(X));[M,N]=size(X);fori=1:1:Mt(i,:)=(0:N-1)/fs;endAc=X.*cos(2*pi*f0*t)+HX.*sin(2*pi*f0*t);As=HX.*cos(2*pi*f0*t)-X.*sin(2*pi*f0*t);subplot(2,1,1);plot(t,Ac(1,:));东北石油大学计算机与信息技术学院李勇勇xlabel('时间(s)');ylabel('窄带低频过程Ac');subplot(2,1,2);plot(t,As(1,:));xlabel('时间(s)');ylabel('窄带低频过程As');运行结果：[Ac,As]=Lowfsignal(rand(10,100),0.3,20)4.编写Matlab程序,对f0=10kHz,带宽500Hz的窄带高斯过程X(t)及低频过程Ac(t),As(t)的功率谱密度进行估计,其中fs=22kHz,采样点数N=10000,滤波器阶数M=200。clc;clearall;closeall;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');Rxw=abs(fft(Rx));Racw=abs(fft(Rac));Rasw=abs(fft(Ras));f=(1:N)/N*fs/2;subplot(3,1,1)plot(f,10*log10(Rxw(1:N)));axis([1,N,-120,20]);subplot(3,1,2)plot(f,10*log10(Racw(1:N)));axis([1,N,-120,20]);subplot(3,1,3)plot(f,10*log10(Rasw(1:N)));axis([1,N,-120,20]);5.编写Matlab函数,建立采样频率fs,中心频率f0的窄带随机过程X(t)的包络A(t),相位和包络的平方样本。注：函数EnvelopPhase(X,f0,fs),用来产生At,Ph,A2.function[At,Ph,A2]=EnvelopPhase(X,f0,fs)HX=imag(hilbert(X));[M,N]=size(X);fori=1:1:Mt(i,:)=(0:N-1)/fs;endAc=X.*cos(2*pi*f0*t)+HX.*sin(1*pi*f0*t);As=HX.*cos(2*pi*f0*t)-X.*sin(2*pi*f0*t);At=(Ac.^2+As.^2).^0.5;Ph=atan(As./Ac);A2=At.^2;subplot(3,1,1);plot(t,At(1,:));xlabel('时间(s)');ylabel('窄带包络At');subplot(3,1,2);plot(t,Ph(1,:));xlabel('时间(s)');ylabel('窄带相位Ph');东北石油大学计算机与信息技术学院李勇勇subplot(3,1,3);plot(t,A2(1,:));xlabel('时间(s)');ylabel('包络平方A2');运行结果：[At,Ph,A2]=EnvelopPhase(rand(10,100),10000,22000)6.编写f0=10kHz,带宽为500Hz,方差为1的窄带高斯过程X(t)的包络A(t),相位,包络平方的产生程序.并对分布情况进行统计,其中fs=22kHz,样本采样点数N=20000,滤波器阶数M=50。clc;clearall;closeall;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);subplot(3,1,1)LA=0:0.05:4.5;hist(At,LA);axis([0,4.5,0,700])subplot(3,1,2)LP=-pi/2:0.05:pi/2;hist(Ph,LP);axis([-pi/2,pi/2,0,500])subplot(3,1,3)LA2=0:0.2:16;hist(A2,LA2);axis([0,16,0,2000])7.编写程序,仿真f0=10kHz,带宽400Hz,归一化方差窄带高斯过程X(t),在三种不同余弦信号情况下,窄带高斯过程加余弦信号的包络A(t),相位,包络平方的分布,其中,fs=22kHz,样本采样点数N=10000,滤波器阶数M=50,三种余弦信号的幅度分别为2,4,8,theta值对应取pi/6,pi/4,pi/3。clc;clearall;closeall;N=10000;f0=10000;deltf=400;fs=22000;M=50;a1=2;a2=4;a3=8;theta1=pi/6;theta2=pi/4;theta3=pi/3;X=Narrowbandsignal(N,f0,deltf,fs,M);X=X/sqrt(var(X));t=(0:N-1)/fs;X1=X+a1*cos(2*pi*f0*t+theta1);X2=X+a2*cos(2*pi*f0*t+theta2);X3=X+a3*cos(2*pi*f0*t+theta3);[At1,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(At1,LA);GA2=hist(At2,LA);GA3=hist(At3,LA);subplot(3,1,1)plot(LA,GA1,'r',LA,GA2,'g',LA,GA3,'b')东北石油大学计算机与信息技术学院李勇勇axis([0,12,0,2000])subplot(3,1,2)LP=-pi/2:0.05:pi/2;GP1=hist((Ph1-theta1),LP);GP2=hist((Ph2-theta2),LP);GP3=hist((Ph3-theta3),LP);plot(LP,GP1,'r',LP,GP2,'g',LP,GP3,'b')axis([-pi/2,pi/2,0,2000])subplot(3,1,3)LA2=0:1:120;GA21=hist(A21,LA2);GA22=hist(A22,LA2);GA23=hist(A23,LA2);GP2=hist((Ph2-theta2),LP);GP3=hist((Ph3-theta3),LP);plot(LA2,GA21,'r',LA2,GA22,'g',LA2,GA23,'b')axis([0,120,0,1200])8.两个随机相位余弦信号,振幅分别为0和2,与均值为零方差为1,中心频率f0=10kHz,带宽500Hz的窄带高斯过程的合成信号,经平方律检波,视频积累8次后,对信号分布进行仿真,其中,fs=22kHz,样本采样点数N=160000,滤波器阶数M=50。注：视频累积8次,即在总序列中,将间隔一定点数(为保证8点不相关)的8个点之和,作为一个点,形成新的序列。clc;clearall;closeall;N=160000;f0=10000;deltf=500;fs=22000;M=50;a=2;theta=pi/3;X=Narrowbandsignal(N,f0,deltf,fs,M);X=X/sqrt(var(X));t=(0:N-1)/fs;X1=X+a*cos(2*pi*f0*t+theta);[At,Ph,A2]=EnvelopPhase(X,f0,fs);[At1,Ph1,A21]=EnvelopPhase(X1,f0,fs);closeall;n=8;m=20;m1=n*m;L=N/m1;G0=zeros(1,L);G1=zeros(1,L);fori=1:1:Lforj=1:1:nG0(i)=G0(i)+A2(i*m1-j*m+1);G1(i)=G1(i)+A21(i*m1-j*m+1);endendly=0:2:120;东北石油大学计算机与信息技术学院李勇勇G0=hist(G0,ly);G1=hist(G1,ly);plot(ly,G0,'r',ly,G1,'b')legend('n=8,a=0','n=8,a=2')