随机信号分析试验(4)

随机信号分析试验（4）

1. 编写程序, 对均值为0方差为1的平稳高斯随机过程样本进行希尔伯特变换，绘制变换前后的概率密度直方图。

注：z=hilbert(x)返回解析信号x+jH(x), 所以imag(z)为x的希尔伯特变换

clc;clear all;close all;

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: 产生窄带信号的滤波器阶数。(M<

function X=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);

for i=1:1:M

t(i,:)=(0:N-1)/fs;

end

Ac=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;clear all;close all;

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);

for i=1:1:M

t(i,:)=(0:N-1)/fs;

end

Ac=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;clear all;close all;

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;clear all;close all;

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;clear all; close all;

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);

close all;

n=8;

m=20;

m1=n*m;

L=N/m1;

G0=zeros(1,L);

G1=zeros(1,L);

for i=1:1:L

for j=1:1:n

G0(i)=G0(i)+A2(i*m1-j*m+1);

G1(i)=G1(i)+A21(i*m1-j*m+1);

end

end

ly=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')