随机信号分析实验报告.doc

随机信号分析实验报告 班级：13050141 姓名： 学号： 日期：2016年5月24日实验一 随机噪声的产生与性能测试一、实验内容 1产生满足均匀分布、高斯分布、指数分布、瑞利分布的随机数，长度为N=1024，并计算这些数的均值、方差、自相关函数、概率密度函数、概率分布函数、功率谱密度，画出时域、频域特性曲线；2编程分别确定当五个均匀分布过程和5个指数分布分别叠加时，结果是否是高斯分布； 3采用幅度为2， 频率为25Hz 的正弦信号为原信号，在其中加入均值为2 ， 方差为0.04 的高斯噪声得到混合随机信号，编程求 的均值、相关函数、协方差函数和方差，并与计算结果进行比较分析。二、实验步骤1.程序N=1024;fs=1000;n=0:N-1;signal=chi2rnd(2,1,N); %rand(1,N)均匀分布 ,randn(1,N)高斯分布,exprnd(2,1,N)指数分布,raylrnd(2,1,N)瑞利分布,chi2rnd(2,1,N)卡方分布signal_mean=mean(signal);signal_var=var(signal);signal_corr=xcorr(signal,signal,unbiased);signal_density=unifpdf(signal,0,1);signal_power=fft(signal_corr);%s,w=periodogram(signal);k1,n1=ksdensity(signal);k2,n2=ksdensity(signal,function,cdf);figure;hist(signal);title(频数直方图);figure;plot(signal);title(均匀分布随机信号曲线);f=n*fs/N; %频率序列figure;plot(abs(signal_power);title(功率幅频);figure;plot(angle(signal_power);title(功率相频);figure;plot(1:2047,signal_corr);title(自相关函数);figure;plot(n1,k1);title(概率密度);figure;plot(n2,k2);title(分布函数);结果（1）均匀分布 （2）高斯分布 （3）指数分布 （4）瑞利分布 （5）卡方分布 2.程序N=1024;signal_1=rand(1,N);signal_2=rand(1,N);signal_3=rand(1,N);signal_4=rand(1,N);signal_5=rand(1,N);signal=signal_1+signal_2+signal_3+signal_4+signal_5;k1,n1=ksdensity(signal);figure(1)subplot(1,2,1);hist(signal);title(叠加均匀分布随机数直方图);subplot(1,2,2);plot(n1,k1);title(叠加均匀分布的概率密度); 结果指数分布叠加 均匀分布叠加 结果：五个均匀分布过程和五个指数分布分别叠加时，结果是高斯分布。3.程序clear all;clc; t=0:320; x=2*sin(2*pi*t*25); z=imnoise(x,gaussian,2,0.04); % z=x+x1; y=trapz(t,z); %y=int(z,x,0,t); subplot(3,2,1),plot(z); title(随机信号序列) meany=mean(z); figure(1)subplot(3,2,3),plot(t,meany,.); title(随机信号均值) vary=var(y); %方差 subplot(3,2,4),plot(t,vary,.); title(随机信号方差) cory=xcorr(z,unbiased);%自相关函数 subplot(3,2,2),plot(cory); title(随机信号自相关函数) covv=cov(y); subplot(3,2,5),plot(t,covv,.);title(随机信号协方差) t=0:0.0005:0.045;X1=sin(2*pi*25*t);%正弦 figure(2)subplot(3,4,1);plot(t,X1);grid title(正弦函数序列); X2=randn(1,length(t); %产生标准正态分布的随机数%X2=normrnd(2,0.04); %产生正态分布随机数subplot(3,4,2); plot(t,X2); title(高斯噪声序列); X=X1+X2; %混合随机信号X(t)subplot(3,4,3); plot(t,X);gridtitle(混合随机信号); meany1=mean(X1); %求均值subplot(3,4,6),plot(t,meany1);title(原信号均值); vary1=var(X1); %求方差subplot(3,4,7),plot(t,vary1);title(原信号方差); 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(混合信号概率密度分布);结果三、实验结果分析随机数产生的理论依据：生成随机数，需要伪随机数序列这个序列通过设定种子和生成算法来确定，RandStream的构造函数方法来完成此任务使用 RandStream.setDefaultStream函数将确定好的序列对象设置为当前matlab使用的序列。然后通过rand等函数就可以用此序列生成随机数了。特征量计算的理论依据：均值是所有值相加求和然后除以这些值的个数，方差是总体各单位标准值与其平均数离差平方。实验二 随机信号通过线性系统的分析一、实验内容 1已知输入信号为，其中(为学号)，为高斯白噪声；系统为一FIR低通滤波器，通带截止频率为，阻带截止频率为，通带最大衰减为40dB，阻带最小衰减为1dB。试分析信号通过滤波器前后的时域信号、频谱、自相关、功率谱密度、期望、方差等性能参数的差异； 2已知平稳随机过程X(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);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(滤波后功率谱);结果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);R_y=conv(R_yx,h);R_y=R_y(41:121);figure(2);subplot(221),stem(R_x);title(Rx); subplot(222),stem(R_xy);title(Rxy); %产生XY的互相关函数subplot(223),stem(R_yx);title(Ryx);subplot(224),stem(R_y);title(Ry); %产生Y的自相关函数S_xy=abs(fft(R_xy);S_xy=fftshift(S_xy);S_yx=fftshift(abs(fft(R_yx);S_y=fftshift(abs(fft(R_y);figure(3);subplot(221),stem(S_x);title(Sx);subplot(222),stem(S_xy);title(Sxy); %产生XY的互谱密度subplot(223),stem(S_yx);title(Syx);subplot(224),stem(S_y);title(Sy); %产生Y的功率谱密度%频域分析法S0_xy=S_x.*H;S0_yx=S_x.*fliplr(H);S0_y=S0_yx.*H;figure(4);subplot(221),stem(S_x);title(Sx);subplot(222),stem(S0_xy);title(S0xy); %产生XY的互谱密度subplot(223),stem(S0_yx);title(S0yx);subplot(224),stem(S0_y);title(S0y); %产生Y的功率谱密度R0_xy=fftshift(abs(ifft(S0_xy);R0_yx=fftshift(abs(ifft(S0_yx);R0_y=fftshift(abs(ifft(S0_y);figure(5);subplot(221),stem(R_x);title(Rx); subplot(222),stem(R0_xy);title(R0xy); %产生XY的互相关函数subplot(223),stem(R0_yx);title(R0yx);subplot(224),stem(R0_y);title(R0y); %产生Y的自相关函数结果三、实验结果分析信号进入系统前的特征是波形的周期明显，具有稳定的周期性，自相关函数的波形窄，存在三个数字角频率，信号进入系统后波形的周期性减弱，甚至看不到周期，自相关函数的波形变宽，只剩下一个数字角频率。实验三 窄带随机信号的产生及分析一、实验内容 1. 产生一输入信号，其中，(为学号)，与一样，为高斯白噪声； 2按下图系统，设计一个低通滤波器，使得通过系统后的输出为窄带信号，设计一包络检波器，提取出，设计一相位检波器，提取出二、实验步骤程序clcfs=16000; %设定采样频率N=1300;n=0:N-1; %取的样本点数t=n/fs; %获得以1/16000为时间间隔采样序列noisy=wgn(1,N,0); %高斯白噪声w1=2*40*pi*1000;%包络的频率（学号40）w0=w1;%a=1+cos(w1*t); %包络函数x=a.*cos(w0*t+a)+noisy/10; %输入信号%以t为横坐标画出x(t)的时域图型figure(1); subplot(2,1,1); plot(n,x); axis(0 140 -3 3);xlabel(采样点);ylabel(X(t)/V);title(窄带随机信号波形);grid on;%求X(t)的统计特性 并画出来disp(X(t)的均值为); Ex=mean(x); disp(Ex);%输出均值hold on; plot(n,Ex,r.); disp(X(t)的方差为);Dx=var(x); disp(Dx);%输出方差hold on; plot(n,Dx,g.);%画出X(t)的概率分布函数each=linspace(min(x),max(x),14); %创建起始值为min(x),终止值为max(x),长度为14的向量nr=hist(x,each); %产生x的分布情况nr=nr/length(x); %归一化subplot(2,1,2); p=polyfit(each,nr,9); %多项式拟合,多项式次数为9bar(each,nr); %绘制直方图hold on; plot(each,nr,g);eachi=-2:0.1:2;nri=polyval(p,eachi);plot(eachi,nri,r)axis tight;title(X(t)概率密度分布);xlabel(X(t);ylabel(P(x);grid on;%对X(t)进行频谱分析Fx=fft(x,N); %对x(t)进行fft变换，在016000区间内得到2N-1个频率值magn=abs(Fx); %求x(t)幅值xangle=angle(Fx); %求X(t)相位labelang=(0:length(x)-1)*16000/length(x); %在016000区间内求横坐标刻度figure(2); plot(labelang,magn*10); %在016000区间内做频谱和相位图axis(0 16000 -0.5 600); xlabel(频率/Hz);ylabel(幅值);title(X(t)频谱图);grid on;%求X(t)的自相关函数c,lags=xcorr(x,coeff); %求自相关函数figure(3); subplot(2,1,1); plot(lags/fs,c); %绘制自相关函数axis tight; xlabel(T);ylabel(Rx(T);title(X(t)的自相关函数);grid on;%求X(t)的功率谱密度long=length(c); Sx=fft(c,long); labelx=(0:long-1)*2*pi;plot_magn=10*log10(abs(Sx);subplot(2,1,2); plot(labelx,plot_magn); %绘制功率谱密度axis tight;xlabel(w);ylabel(Sx(w);title(X(t)的功率谱密度);grid on;%窄带系统检测z1=2.*cos(2*pi*4000*t); z2=-2.*sin(2*pi*4000*t); Ac=z1.*x; %Ac分量As=z2.*x; %As分量y=Ac.*cos(2*pi*4000*t)-As.*sin(2*pi*4000*t);%滤波器设计f_p=1000;f_s=1600;R_p=1;R_s=35; %设定滤波器参数; 通、阻带截止频率,通、阻带衰减Ws=2*f_s/fs;Wp=2*f_p/fs; %频率归一化n,Wn=buttord(Wp,Ws,R_p,R_s); %采用巴特沃思滤波器b,a=butter(n,Wn); %求得滤波器传输函数的多项式系数figure(4);H,W=freqz(b,a); %求得滤波器传输函数的幅频特性subplot(2,1,1); plot(W*fs/(2*pi),abs(H); %在02pi区间内作幅度谱title(低通滤波器幅度谱); grid on;subplot(2,1,2); plot(W*fs/(2*pi),angle(H); %在02pi区间内作相位谱title(低通滤波器相位谱); grid on;%求Ac(t)滤波后的统计特性mc=filter(b,a,Ac); %上支路通过滤波器 Ac(t)disp(Ac(t)的均值);Eh=mean(mc) %求均值disp(Ac(t)的均方值是);E2h=mc*mc/N %求均方值disp(Ac(t)的方差);Dh=var(mc) %求方差%画Ac(t)的 波形figure(6); subplot(2,1,1); n=0:N-1; plot(n,mc); axis(0 300 -1 1);xlabel(采样点);ylabel(幅值);title(Ac(t)的时域波形);grid on;%画Ac(t)的 图yc=fft(mc,length(mc); %作傅里叶变换longc=length(yc); %求傅里叶变换后的序列长度labelx=(0:longc-1)*16000/longc;magnl=abs(yc); %取幅值subplot(2,1,2); plot(labelx,magnl); %绘制幅频曲线axis tight; xlabel(频率(Hz); ylabel(幅值); title(Ac(t)频谱图); grid on;%求Ac(t)的 函数c1,lags1=xcorr(mc,coeff); %求自相关函数figure(7); subplot(2,1,1); plot(lags1/fs,c1); %绘制自相关函数曲线xlabel(T);ylabel(Rx(T);axis tight;title(Ac(t)的自相关函数);grid on;%求Ac(t)的 功率谱Sac=fft(c1,length(c1); % 求频率响应magnc=abs(Sac); % 求幅频响应long=length(Sac); % 频谱序列长度labelc=(0:long-1)*16000/long;subplot(2,1,2); plot(labelc,10*log10(magnc); % 绘制双边功率谱xlabel(频率(Hz);ylabel(功率谱(dbW);axis tight;title(Ac(t)的双边功率谱);grid on;%求得As(t)的统计特性ms=filter(b,a,As); %对下支路信号进行滤波得As(t)disp(As(t)的均值); Eh=mean(ms) % 求均值disp(As(t)的均方值是); E2h=ms*ms/N % 求均方值disp(As(t)的方差); Dh=var(ms) % 求方差%作As(t)的时域波形figure(8);subplot(2,1,1); n=0:N-1;plot(n,ms); % 绘制时域波形axis(0 300 -0.5 2); xlabel(采样点);ylabel(幅值);title(As(t)的时域波形);grid on;%对As(t)进行FFT变换并做频谱图ys=fft(ms,length(ms); % 进行傅里叶变换longs=length(ys); % 求傅里叶变换后的序列长度labelx=(0:longs-1)*16000/longs;magn2=abs(ys); % 取幅值subplot(2,1,2); plot(labelx,magn2); % 绘制幅频曲线axis tight;xlabel(频率(Hz);ylabel(幅值);title(As(t)的频谱图);grid on;%求As(t)的自相关函数c2,lags2=xcorr(ms,coeff); %求出As(t)的自相关序列figure(9);subplot(2,1,1);plot(lags2/fs,c2); %绘制自相关函数曲线xlabel(T);ylabel(Rx(T);axis tight;title(As(t)的自相关函数);grid on;%求As(t)的双边功率谱Sas=fft(c2,length(c2); %对As(t)的自相关函数进行傅里叶变换magnc=abs(Sac); %求As(t)的双边功率谱幅值long=length(Sas); %求傅里叶变换后的序列长度labels=(0:long-1)*16000/long;subplot(2,1,2); plot(labelc,10*log10(magnc); %画As(t)的自相关函数频谱xlabel(频率(Hz);ylabel(功率谱(dbW);axis tight;title(As(t)的双边功率谱); % 求y(t)的统计特性disp(输出信号Y(t)的均值);Eh=mean(y) % 求均值disp(输出信号Y(t)的均方值);E2h=y*y/N %求均方值disp(输出信号Y(t)的方差);Dh=var(y) % 求方差%作输出信号Y(t)的时域波形figure(10); subplot(2,1,1);n=0:N-1;plot(n,y); axis(0 150 -2 2);xlabel(采样点);ylabel(幅值);title(Y(t)的时域波形);grid on;%进行FFT变换并做频谱图yy=fft(y,length(y); %进行傅里叶变换longy=length(yy); %求傅里叶变换后的序列长度labelx=(0:longy-1)*16000/longy;magn3=abs(yy); %取幅值subplot(2,1,2); plot(labelx,magn3); %绘制幅频曲线axis tight;xlabel(频率(Hz);ylabel(幅值);title(Y(t)的频谱图);grid on;%求输出信号Y(t)的自相关函数c3,lags3=xcorr(y,coeff); %求自相关函数figure(11); subplot(2,1,1); plot(lags3/fs,c3); %绘制自相关函数曲线xlabel(T);ylabel(Rx(T);axis tight;title(Y(t)的自相关函数);grid on;%求输出信号Y(t)的双边功率谱Sy=fft(c3,length(c3); %进行傅里叶变换magny=abs(Sy); %取幅值long=length(Sy);labely=(0:long-1)*16000/long; subplot(2,1,2); plot(labely,10*log10(magny); %*画Y(t)的功率谱密度xlabel(频率(Hz);ylabel(功率谱(dbW);axis tight;title(Y(t)的双边功率谱);grid on;结果3、 实验结果分析由于高斯白噪声Xt是标准正态的，所以均值趋近于零，而t是由Xt通过一个线性系统（低通滤波器）得到的，所以输出均值不变，仍为零，从程序运行结果可以看出，Xt，at，t均值都趋近于零。由自相关函数图形可看出，中心点上相关程度最高，在其他地方，自相关函数接近于零。高斯白噪声经过低通滤波器后：通过二者的自相关图像对比可知，限带白噪声的自相关图像是一个振荡图像，与先前的分析相符合，通过二者的功率谱密度图像对比可知，限带白噪声的功率谱具有选择性，只有低频成分通过。实验四希尔伯特变换的应用一、实验内容 1. 产生一输入信号，其中，(为学号)，与一样，为高斯白噪声；设计一个低通滤波器，使得通过系统后的输出为窄带信号； 2利用希尔伯特变换实现单边带的调幅，如下图所示二、实验步骤程序%1）参数设定clc clear allfs=15000; %采样频率ts=1/fs; %采样周期t=0:ts:0.01; %时间序列df=0.2; %采样分辨率M=2048; %频率点数fc = 4000; %载波频率Lt=length(t); %时间序列长度%L=2*min(at); %R=2*max(abs(at); %2）产生高斯白噪声n(t)并进行频谱分析nt = wgn(1,length(t),0.1);n_1=nt/max(abs(nt); %噪声figure(1);subplot(211);plot(t,n_1);title(高斯白噪声nt信号);xlabel(t/s);ylabel(幅度/v);grid on; n=0:M-1; %t=n/fs; %时间序列y0=fft(n_1,M);mag0=(abs(y0);f=n*fs/(1000*M);subplot(212);plot(f,mag0);title(高斯白噪声频谱分析);xlabel(f/KHz);ylabel(幅度/v);axis(0 10 0 20);grid on; %3）产生基带信号s(t)并进行频谱分析%st=sin(1000*2*pi*t);w1 = 2*40*pi*1000; %学号40w0 = w1;at = 1+cos(w1*t);wt = cos(w0*t+at);st = wt.*at;subplot(211);plot(t,st);title(初始信号); xlabel(t/s);ylabel(幅度/v);grid on; y1=fft(st,M);mag1=(abs(y1);f=n*fs/(1000*M);subplot(212);plot(f,mag1);title(初始信号频谱分析);xlabel(f/KHz);ylabel(幅度/v);grid on;axis(0 10 0 100); %4）调制信号（s(t)+n(t)）进行频谱分析xt=st+n_1;subplot(211);plot(t,xt);title(调制信号xt=st+nt（初始信号+噪声）); xlabel(t/s);ylabel(幅度/v);grid on; y3=fft(xt,M);mag3=(abs(y3);f=n*fs/(1000*M);subplot(212);plot(f,mag3);title(调制信号频谱分析);xlabel(f/KHz);ylabel(幅度/v);axis(0 10 0 100);grid on; %5）调制信号通过滤波器哟后a点的信号分析wp=2*2200/fs; %通带边界频率ws=2*2800/fs; %阻带边界频率Rp=1; %通带最大衰减度As=30; %阻带最小衰减度V,wc=buttord(wp,ws,Rp,As);B,A=butter(V,wc);H,W=freqz(B,A);at=filter(B,A,xt); %经过低通