cory1=xcorr(X1,'unbiased'); %求自相关 subplot(3,4,8),plot(cory1); title('原信号自相关函数');
meany=mean(X); %求均值 subplot(3,4,10),plot(t,meany); title('混合信号均值');
vary=var(X); %求方差 subplot(3,4,11),plot(t,vary); title('混合信号方差')
cory=xcorr(X,'unbiased'); %求自相关 subplot(3,4,12),plot(cory); title('混合信号自相关函数')
covy=cov(X1,X); %求X1与X的卷积 subplot(3,4,4),plot(covy); title('协方差');
[f1,xi]=ksdensity(X1); %估计样本的概率密度函数 subplot(3,4,5); plot(xi,f1); title('原信号的概率密度分布)');
[f2,xi]=ksdensity(X); %估计样本的概率密度函数 subplot(3,4,9); plot(xi,f2); title('混合信号概率密度分布');
结果
11
三、实验结果分析
随机数产生的理论依据：
生成随机数，需要伪随机数序列这个序列通过设定种子和生成算法来确定，RandStream的构造函数方法来完成此任务使用 RandStream.setDefaultStream函数将确定好的序列对象设置为当前matlab使用的序列。然后通过rand等函数就可以用此序列生成随机数了。 特征量计算的理论依据：
均值是所有值相加求和然后除以这些值的个数，方差是总体各单位标准值与其平均数离差平方。
实验二 随机信号通过线性系统的分析
一、实验内容
1．已知输入信号为X(t)?sin2?f1t?sin2?f2t?sin2?f3t?N(t)，其中
f1?nkHz(n为学号)，f2?2nkHz，f3?3nkHz，N(t)为高斯白噪声；系统为一FIR低通滤波器，通带截止频率为
f1，f通带最大衰减为40dB，
阻带截止频率为2，
阻带最小衰减为1dB。试分析信号通过滤波器前后的时域信号、频谱、自相关、功率谱密度、期望、方差等性能参数的差异；
22R(m)???(m),??5；线性系 2．已知平稳随机过程X(n)的相关函数为：
统的单位冲击响应为
h(k)?rk,k?0,r?1?12
1n?1（n为学号）。试分别用时域分析
法和频域分析法分析信号通过滤波器前后的时域信号、频谱、自相关、功率谱密度、期望、方差等性能参数的差异；求解输入输出的互相关函数、互功率谱密度。
二、实验步骤 1.程序
N=1024;
Ts=1/200000; n=1:1:N;
w=randn(1,N); a=40; %学号 f1=a*1000; f2=a*1000*2; f3=a*1000*3;
x=sin(2*pi*f1*n*Ts)+sin(2*pi*f2*n*Ts)+sin(2*pi*f3*n*Ts)+w(n); x_mean=mean(x); x_var=var(x); R=xcorr(x);
[S,W]=periodogram(x); m=-(N-1):1:(N-1); %f=n/(N*Ts); figure(1)
subplot(3,1,1);plot(n,x);title('原始信号'); subplot(3,1,2);plot(m,R);title('原始自相关'); subplot(3,1,3);plot(W,S);title('原始功率谱');
fp=f1;fs=f2;ap=1;as=40; %数字滤波器指标 T=Ts,fsa=1/T; %采样频率与间隔
wp=2*pi*fp/fsa;ws=2*pi*fs/fsa; %转换为数字角频率
Wp=2/T*tan(wp/2);Ws=2/T*tan(ws/2); %由数字角频率转换为模拟角频率
[N1,Wc]=buttord(Wp,Ws,ap,as,'s'); %获取模拟滤波器的阶数和3dB截止频率 [Z,P,K]=buttap(N1); %归一化模拟滤波器模型的零极点形式参数 [B,A]=zp2tf(Z,P,K); %归一化模拟滤波器传递函数的系数
[Bl,Al]=lp2lp(B,A,Wc); %把模拟滤波器原型转换成截至频率为Wc的低通滤波器 [b,a]=bilinear(Bl,Al,fsa);%用双线性变换法实现模拟滤波器到数字滤波器的转换 [H,w]=freqz(b,a); %获取频率响应 figure(2)
plot(w*fsa/(2*pi),abs(H));grid; %绘制频率响应曲线 xlabel('频率(Hz)');ylabel('频率响应幅度'); title('FIR低通滤波器频率响应');
y=filter(b,a,x); y_mean=mean(y); y_var=var(y); Ry=xcorr(y);
13
[Sy,Wy]=periodogram(y); m=-(N-1):1:(N-1); %f=n/(N*Ts); figure(3)
subplot(3,1,1);plot(n,y);title('滤波后信号'); subplot(3,1,2);plot(m,Ry);title('滤波后自相关'); subplot(3,1,3);plot(Wy,Sy);title('滤波后功率谱');
结果
14
2.程序
clc;
R_x=zeros(1,81);R_x(41)=1; %产生全零序列，并在序列中间设置一个冲激 S_x=fftshift(abs(fft(R_x))); %进行快速傅里叶变换 No=40; %学号 r=1-1/(No+1); h0=zeros(1,40); i=1:41;
h1=r.^i;
h=[h0,h1]; %把h0，h1并在一起
H=fftshift(abs(fft(h)));%进行快速傅里叶变换 m_x=0; sigma_x=R_x(41); P_x=R_x(41); figure(1);
subplot(221),stem(R_x),title('Rx'); subplot(222),stem(S_x),title('Sx'); subplot(223),stem(h),title('h'); subplot(224),stem(H),title('H');
%%%%%%%%%%%%%%%%%% %时域法求解
R_xy=conv(R_x,h);R_xy=R_xy(41:121);
R_yx=conv(R_x,fliplr(h));R_yx=R_yx(41:121);
15