随机信号分析实验报告

1、随机信号分析实验报告 学号： 姓名： 日期：目录实验一 随机噪声的产生与性能测试1一实验内容1二实验步骤11.程序：12. 结果：33. 实验结果分析6实验二 随机信号通过线性系统的分析6一实验内容6二实验步骤71.程序72.结果93.实验结果分析：13实验三 窄带随机信号的产生及分析13一实验内容：13二 实验步骤131. 程序132. 结果153.实验结果分析18实验四希尔伯特变换的应用18一实验内容：18二实验步骤：191.程序192. 结果233.实验结果分析25实验一 随机噪声的产生与性能测试一实验内容1产生满足均匀分布、高斯分布、指数分布、瑞利分布、的随机数，长度为N=1024，并

2、计算这些书的均值、方差、自相关函数、概率密度函数、概率分布函数、功率谱密度，画出时域、频域特性曲线；2编程分别确定当五个均匀分布过程和5个指数分布分别叠加时，结果是否是高斯分布； 3采用幅度为2， 频率为25Hz 的正弦信号为原信号，在其中加入均值为2 ， 方差为0.04 的高斯噪声得到混合随机信号，编程求 的均值、相关函数、协方盖函数和方差，并与计算结果进行比较分析。二实验步骤1.程序：（1）内容1实验代码：N=1024;fs=1000;n=0:N-1;% signal=rand(1,N);%均匀% signal=randn(1,N);%高斯% signal=exprnd(2,1,N);%指

3、数% signal=raylrnd(2,1,N);%瑞利signal=chi2rnd(2,1,N);%卡方%signal_mean=mean(signal);%qiu jun zhimeany=mean(signal); figure(1)subplot(3,2,1),plot(0:1023,signal); title('原始信号') subplot(3,2,2),plot(0:1023,meany); title('随机信号均值') vary=var(signal); %方差 subplot(3,2,3),plot(n,vary,'.');

4、title('随机信号方差') cory=xcorr(signal,'unbiased');%自相关函数 subplot(3,2,4),plot(cory); title('随机信号自相关函数') f1,xi=ksdensity(signal); %subplot(3,2,5); plot(xi,f1); title('原信号的概率密度函数'); %p=unidcdf(n,N);k2,n2=ksdensity(signal,'function','cdf');subplot(3,2,6); plot

5、(n2,k2);title('原信号的概率分布函数'); figure1s,w=periodogram(signal);plot(w,s);title('功率谱密度函数')figure(2)xn=signal;%信号cxn=xcorr(xn,'unbiased'); %计算序列的自相关函数CXk=fft(cxn,N);%对自相关函数进行傅里叶变换Pxx=abs(CXk);index=0:round(N/2-1);k=index*fs/N;%频谱序列subplot(312),plot(k,Pxx(index+1);%绘制功率谱曲线title(

6、9;功率谱密度函数');（2）内容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(2,2,1);hist(signal);title('叠加均匀分布随机数直方图');subplot(2,2,2);plot(n1,k1);t

7、itle('叠加均匀分布的概率密度');signal_6=exprnd(2,1,N);signal_7=exprnd(2,1,N);signal_8=exprnd(2,1,N);signal_9=exprnd(2,1,N);signal_10=exprnd(2,1,N);signal2=signal_6+signal_7+signal_8+signal_9+signal_10;k2,n2=ksdensity(signal);subplot(2,2,3);hist(signal2);title('叠加指数分布随机数直方图');subplot(2,2,4);plot

8、(n2,k2);title('叠加指数分布的概率密度');（3）内容3实验代码：clc; clear all; t=0:320; x=2*sin(2*pi*t*25); %y=trapz(t,z); %积分%z=int(x,t); y=-2/(2*pi*25)*cos(2*pi*t*25);z=imnoise(y,'gaussian',2,4); subplot(3,2,1),plot(t,z); axis(0 320 -3 3);title('随机信号序列') meany=mean(z); figure(1)subplot(3,2,3),plo

9、t(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('随机信号协方差'

10、)2. 结果：（1）程序1运行结果均匀分布：高斯分布：指数分布：瑞利分布：卡方分布：（2） 程序2运行结果（3）程序3运行结果 实验二 随机信号通过线性系统的分析一实验内容1已知输入信号为，其中(为学号)，为高斯白噪声；系统为一FIR低通滤波器，通带截止频率为，阻带截止频率为，通带最大衰减为40dB，阻带最小衰减为1dB。试分析信号通过滤波器前后的时域信号、频谱、自相关、功率谱密度、期望、方差等性能参数的差异；2已知平稳随机过程X(n)的相关函数为：；线性系统的单位冲击响应为（为学号）。试分别用时域分析法和频域分析法分析信号通过滤波器前后的时域信号、频谱、自相关、功率谱密度、期望、方差等性能参

11、数的差异；求解输入输出的互相关函数、互功率谱密度。二实验步骤1.程序（1）内容1实验代码N=1024;Ts=1/50000;n=1:1:N;w=randn(1,N);a=32;f1=a*500;f2=a*500*2;f3=a*500*3;x=sin(2*pi*f1*n*Ts)+sin(2*pi*f2*n*Ts)+sin(2*pi*f3*n*Ts)+w(n);R=xcorr(x);S,W=periodogram(x);%功率谱密度m=-(N-1):1:(N-1);figure(1)subplot(3,2,1);plot(n,x);title('滤波前时域')subplot(3,2

12、,2);plot(m,R);title('滤波前自相关')subplot(3,2,3);plot(W,S);title('滤波前功率谱密度')meany=mean(x);subplot(3,2,4),plot(0:1023,meany); title('滤波前均值') vary=var(x); %方差 subplot(3,2,5),plot(n,vary,'.'); title('滤波前方差') f=n/(N*Ts);FFT=fftshift(abs(fft(x);subplot(3,2,6),plot(f,FFT

13、,'.'); 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); %归一

14、化模拟滤波器传递函数的系数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('频率响应幅度');y=filter(b,a,x);%滤波器处理Ry=xcorr(y);Sy,Wy=periodogram(y);m=-(N-1):1:(N-1

15、);figure(3)subplot(3,2,1);plot(n,y);title('滤波后时域')subplot(3,2,2);plot(m,Ry);title('滤波后自相关')subplot(3,2,3);plot(Wy,Sy);title('滤波后功率谱密度')meany=mean(y);subplot(3,2,4),plot(0:1023,meany); title('滤波后均值') vary=var(y); %方差 subplot(3,2,5),plot(n,vary,'.'); title('

16、滤波后方差') f=n/(N*Ts);FFTY=fftshift(abs(fft(y);subplot(3,2,6),plot(f,FFTY,'.'); title('滤波后频谱') （2）内容1实验代码clc;R_x=zeros(1,81);R_x(41)=1; %S_x=fftshift(abs(fft(R_x); %No=41; %r=1-1/(No+1);h0=zeros(1,40);i=1:41;h1=r.i; h=h0,h1; %H=fftshift(abs(fft(h);%m_x=0; %sigma_x=R_x(41); P_x=R_x(4

17、1);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);figur

18、e(2),subplot(221),stem(R_x);title('Rx'); subplot(222),stem(R_xy);title('Rxy'); %subplot(223),stem(R_yx);title('Ryx');subplot(224),stem(R_y);title('Ry'); %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

19、),stem(S_x);title('Sx');subplot(222),stem(S_xy);title('Sxy'); %subplot(223),stem(S_yx);title('Syx');subplot(224),stem(S_y);title('Sy'); %频域分析法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);

20、title('S0xy'); %subplot(223),stem(S0_yx);title('S0yx');subplot(224),stem(S0_y);title('S0y'); %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('R

21、0xy'); %subplot(223),stem(R0_yx);title('R0yx');subplot(224),stem(R0_y);title('R0y'); %2.结果（1）（2）实验三 窄带随机信号的产生及分析一实验内容：1. 产生一输入信号，其中，(为学号)，与一样，为高斯白噪声；2按下图系统，设计一个低通滤波器，使得通过系统后的输出为窄带信号，设计一包络检波器，提取出，设计一相位检波器，提取出2 实验步骤1. 程序clcfs=16000; %设定采样频率N=1300;n=0:N-1; %取的样本点数t=n/fs; %获得以1/16000

22、为时间间隔采样序列noisy=wgn(1,N,0); %1*N强度为0的高斯噪声a=cos(32*2*pi*1000*t*0.2); %x=(1+a).*cos(32*2*pi*1000*t*0.2+1+a)+noisy/10; %输入信号x（t） %以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;subplot(2,1,2); n=0:N

23、-1; plot(n,x); title('x(t)的时域波形'); %滤波器设计f1=2800;f3=5200;%通带截止频率上下限 fsl=2500;fsh=5500;%阻带截止频率上下限 rp=0.1;rs=30;%通带边衰减DB值和阻带边衰减DB值 wp1=2*pi*f1/fs; wp3=2*pi*f3/fs; wsl=2*pi*fsl/fs; wsh=2*pi*fsh/fs; wp=wp1 wp3; ws=wsl,wsh; n,wn=cheb1ord(ws/pi,wp/pi,rp,rs);%切比雪夫1型带通滤波器 bz1,az1=cheby1(n,rp,wp/pi);

24、%切比雪夫1型带通滤波器 h,w=freqz(bz1,az1,256,fs); Wt_dai=filter(bz1,az1,x);figure(2); subplot(2,1,1); plot(w,abs(h);%绘制带通滤波器的幅频曲线 xlabel('频率/HZ'); ylabel('h/db'); title('带通滤波器的幅频曲线'); grid on; subplot(2,1,2); plot(w,angle(h);%绘制带通滤波器的相频曲线 xlabel('频率/HZ'); ylabel('w/db')

25、; title('带通滤波器的相频曲线'); %figure(3)n=0:N-1;subplot(211),plot(n,Wt_dai);axis(0 200 -3 3);title('W(T)时域波形');len=length(Wt_dai);ff=(0:len-1)*16000/len;subplot(212),plot(ff,abs(Wt_dai);title('W(T)频域波形'); %包络检波后的统计特性ac=abs(hilbert(Wt_dai); %包络检波器%画Ac(t)的波形figure(4); subplot(2,1,1);

26、n=0:N-1; plot(n,ac); axis(0 300 -3 3);xlabel('采样点');ylabel('幅值');title('Ac(t)的时域波形');grid on;%进行fft变换并做频谱图yc=fft(ac,length(ac); %longc=length(ac); %求傅里叶变换后的序列长度labelx=(0:longc-1)*16000/longc;magnl=abs(ac); %subplot(2,1,2); plot(labelx,magnl); %axis (0 1000 -3 3); xlabel('

27、频率(Hz)'); ylabel('幅值'); title('Ac(t)频谱图'); grid on; %求得As(t)的统计特性x2=Wt_dai.*cos(32*2*pi*1000*t*0.2); %滤波器设计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(5);H,

28、W=freqz(b,a); %求得滤波器传输函数的幅频特性%Wt_di=filter(b,a,x);subplot(2,2,1); plot(W*fs/(2*pi),abs(H); %在02pi区间内作幅度谱title('低通滤波器幅度谱'); grid on;subplot(2,2,2); plot(W*fs/(2*pi),angle(H); %在02pi区间内作相位谱title('低通滤波器相位谱'); grid on; as=filter(b,a,x2); %相位检波%作As(t)的时域波形subplot(2,2,3); n=0:N-1;plot(n,as

29、); % axis(0 400 -0.5 2); xlabel('采样点');ylabel('幅值');title('As(t)的时域波形');grid on;%对As(t)进行FFT变换并做频谱图ys=fft(as,length(as); % longs=length(ys); % labelx=(0:longs-1)*16000/longs;magn2=abs(ys); % subplot(2,2,4); plot(labelx,magn2); % axis tight;xlabel('频率(Hz)');ylabel('

30、;幅值');title('As(t)的频谱图');grid on;2. 结果实验四希尔伯特变换的应用一实验内容：1. 产生一输入信号，其中，(为学号)，与一样，为高斯白噪声；设计一个低通滤波器，使得通过系统后的输出为窄带信号；2利用希尔伯特变换实现单边带的调幅，如下图所示 二实验步骤：1.程序%1）参数设定clc clear allfs=15000; %采样频率ts=1/fs; %采样周期t=0:ts:0.01; %时间序列df=0.2; %采样分辨率M=2048; %频率点数n=1:M-1;fc = 4000; %载波频率Lt=length(t); %时间序列长度%2

31、）产生高斯白噪声n(t)并进行频谱分析nt = wgn(1,length(t),0.1);n_1=nt/max(abs(nt); %噪声figure(1);subplot(321);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(322);plot(f,mag0);title('高斯白噪声频谱分析'

32、;);xlabel('f/KHz');ylabel('幅度/v');axis(0 10 0 20);grid on;%3）产生基带信号s(t)并进行频谱分析a=1+cos(32*1000*2*pi*t);%学号改动st=a.*cos(32*1000*2*pi*t+a);%学号改动subplot(323);plot(t,st);title('初始信号'); xlabel('t/s');ylabel('幅度/v');grid on;y1=fft(st,M);mag1=(abs(y1);f=n*fs/(1000*M);s

33、ubplot(324);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(325);plot(t,xt);title('调制信号xt=st+nt（初始信号+噪声）'); xlabel('t/s');ylabel('幅度/v');grid on;y3=fft(xt,M);mag3=(abs(y3)

34、;f=n*fs/(1000*M);subplot(326);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

35、(B,A,xt); %经过低通滤波器的a点信号figure(2)subplot(321);plot(W,abs(H);title('a点低通滤波器信号'); xlabel('t/s');ylabel('幅度/v');grid on;meany=mean(at); subplot(322),plot(0:1023,meany); title('a点均值') vary=var(at); %方差 subplot(323),plot(n,vary,'.'); title('a点方差') cory=xcorr

36、(at,'unbiased');%自相关函数 subplot(324),plot(cory); title('a点自相关函数')s,w=periodogram(at);subplot(325),plot(w,s);title('a点功率谱密度函数')%6）信号经过希尔伯特变换产生SSB调制L=2*min(at); R=2*max(abs(at);figure(3);subplot(3,2,1);plot(t,at);title('经过滤波器后的调制信号')xlabel('t/s');ylabel('幅度/v

37、');grid on;%pause %按任意键可看到调制信号的曲线c1=cos(2*pi*fc*t);c2=sin(2*pi*fc*t);subplot(3,2,3);u1=at(1:Lt).*c1(1:Lt)+imag(hilbert(at(1:Lt).*c2(1:Lt);plot(t,u1);title('下边带调制信号');xlabel('t/s');ylabel('幅度/v');grid on;axis(0 0.01 -R R);y2=fft(u1,M);mag2=(abs(y2);f=n*fs/(1000*M);subplot(3,2,4);plot(f,mag2);title('下边带频域信号');xlabel('f/KHz');ylabel('幅度/v');grid on;axis(0 8 0 100);u2=at(1:Lt).*c1(1:Lt)-imag(hilbert(at(1:Lt).*c2(1:Lt);subplot(3,2,5);plot(t,u2);title('上边带调制信号');xlabel('t/s');ylabel('幅度/v'