完整word版平稳随机过程的采样和插值.docx

1、完整word版平稳随机过程的采样和插值随机信号实验平稳随机过程的采样和插值 一 实验目的了解确定信号的采样与平稳随机信号的采样之间的关系，掌握信号的采样及分析方法。二 . 实验原理 确定信号的采样符合香农定理，那么随机信号的采样是否符合香农定理呢？答案是定的。香农定理可以推广到随机信号的采样。若X(t)为平稳随机过程，且具有零均值，它的功率谱密度限于（-，+）之间。当满足条件时，便可将X(t)按它的振幅样本展开为：上式就是平稳随机过程的采样定理。式中T为采样周期。三实验任务与要求 程序用matlab或c/c+语言编写和仿真。系统框图如图29、图30所示： 图29 抽样系统框图 图30 插值系统

2、框图 输入信号x(t)：x(t)=正弦波信号+n(t)，频率为100Hz的正弦波信号，幅值为1v，n(t)为白噪声。计算输入信号的均值、均方值、方差、相关函数、概率密度、频谱及功率谱密度、相关函数。 低通滤波器设计低通滤波器技术要求：通带截止频率1KHz阻带截止频率2KHz。过渡带:1KHz阻带衰减：35DB通带衰减：1DB采样频率：44.1KHz 计算经低通滤波器后信号的均值、均方值、方差、相关函数、频谱及功率谱密度、相关函数。 对输入信号进行抽样：采样频率8000Hz。每间隔4个点和每间隔8个点各抽样一次。计算抽样信号的均值、均方值、方差、相关函数、频谱及功率谱密度、相关函数。 对采样信号

3、进行插值：每一个间隔插入4个值和每一个间隔插入8个值。采样频率8000Hz。计算插值信号的均值、均方值、方差、相关函数、频谱及功率谱密度、相关函数。 对采样前后、插值前插值后信号进行比较。观察在采样频率不变的情况下，信号频谱的变化和频谱的周期延拓性。 讨论X(n)的自相关函数、功率谱密度与X(t)的自相关函数、功率谱密度之间的关系。四实验步骤及结果 1输入信号x(t)（1）x(t)=正弦波信号+n(t)，频率为100Hz的正弦波信号，幅值为1v，其中n（t）为高斯白噪声。 Fs=40000;Ns=5000;n=0:Ns-1;t=n/Fs;两点间距取1/40000，共取N=5000个采样点,n为

4、长度为N的序列，t为采样时间，作为信号图像横坐标sine=sin(2*pi*100*t); 产生频率为100hz的正弦波noise=randn(1,length(t); 产生高斯白噪声 x=sine+noise; 信号合成为输入信号figure(1);plot(t,x);xlabel(t),ylabel(x);title(输入信号);产生的输入信号如下：（2）计算输入信号的均值、均方值、方差 均值: x_mean = 0.0072均方值: x_var = 1.4747 方差 ： input_fangcha =1.4748 （3） x(t)的相关函数 利用MALTLAB中的xcorr函数 x_x

5、corr=xcorr(x);t1=(-Ns+1:Ns-1)/Fs;%相关函数图像横坐标figure(2);plot(t1,x_xcorr);%相关函数xlabel(t),ylabel(R);title(输入信号相关函数); （4） x(t)的概率密度 利用ksdensity函数 F1,y1=ksdensity(x);%y1为输入信号所有取值，F1为每个值相应概率figure(3);plot(y1,F1);%概率密度xlabel(y),ylabel(F);title(输入信号概率密度); （5） x(t)的频谱 利用fft函数 x_spectra=fft(x);f=Fs*n/Ns;%频域横坐标f

6、igure(4);plot(f(1:300),abs(x_spectra(1:300);xlabel(f),ylabel(Y);title(输入信号频谱);（6） x(t)的功率谱密度 用功率谱的概念求，即频谱函数的平方求时间平均。P1=x_spectra.*conj(x_spectra)/Ns;%频谱函数的平方求时间平均figure(5);plot(f(1:300),abs(P1(1:300);%功率谱密度xlabel(f),ylabel(S);title(输入信号功率谱密度) 2滤波器的设计 设计思路是：用巴特沃斯模拟滤波器生成数字滤波器，并画出滤波器的频谱。 程序如下：fp=1000;f

7、s=2000;rp=1;rs=35;Fs1=40000;wp=2*pi*fp/Fs1;ws=2*pi*fs/Fs1;wap=tan(wp/2);was=tan(ws/2);Fs1=Fs1/Fs1;N,Wn=buttord(wap,was,rp,rs,s);z,p,k=buttap(N);bp,ap=zp2tf(z,p,k);%得到传输函数bs,as=lp2lp(bp,ap,wap);%低通到低通，频谱变换bz,az=bilinear(bs,as,Fs1/2); %将模拟滤波器传递函数转为数字滤波器传递函数H,w=freqz(bz,az,256,Fs1*10000);%特性分析figure(6)

8、plot(w,abs(H);%画出数字滤波器的波形图title(低通滤波器的频谱);xlabel(f/hz);grid on;3x(t)通过滤波器的信号y （1）输出y的波形及均值，均方值，方差 y=filter(bz,az,x)均值: 0.0416 均方值: 0.7248方差 ： 0.7265 （2）输出y的相关函数 程序类似求x(t)的相关函数，其波形如下： （3） y的概率密度 程序类似求x(t)的概率密度，其波形如下： （4） y的频谱 程序类似求x(t)的频谱，其波形如下：（5）y的功率谱密度 4对y信号进行抽样 采样频率8000Hz 若8000hz采样，波形输出程序为：sample

9、1=zeros(1,ceil(length(t)/(Fs/fs1);i=1;for k=1:length(t) if mod(k,(Fs/fs1)=1 sample1(i)=y(k); i=i+1; endendfigure(12);stem(t2,sample1);%除去相邻采样点间零点的采样信号xlabel(t),ylabel(y);title(采样信号);则波形为： （1） 每间隔4个点抽样一次 4采样波形及均值，均方值，方差程序代码如下：sample2=zeros(1,ceil(length(t)/(Fs/fs1);for k=1:length(t2) if mod(k,5)=1 sa

10、mple2(k)=sample(k); else sample2(k)=0; endendsample3=zeros(1,ceil(length(t)/(Fs/fs1)/5);i=1;for k=1:5:length(t2) sample3(i)=sample1(k); i=i+1;end figure(13);stem(1:length(t2)/5)/300,sample3);%除去相邻采样点间零点的采样信号xlabel(t),ylabel(y);title(4采样信号);sample_junzhi4=mean(sample2)%均值sample_fangcha4=std(sample2)%

11、方差sample_junfangzhi4=sample_fangcha4+sample_junzhi4.2%均方值4采样信号为：均值：-0.0125方差：0.3212均方值：0.32134采样信号的相关函数 代码函数如上 4采样信号的概率密度4采样信号的频谱4采样信号的功率谱密度（2） 每间隔8个点抽样一次 程序和4采样类似，只是间隔8个点采样信号 均值：-7.0098e-004均方值：0.1099方差：0.1099 8采样信号的相关函数 8采样信号的概率密度 8采样信号的频谱 8采样信号的功率谱密度 4对信号进行插值 采样频率为8000Hz插值思路：用三次多项式插值法 （1）每一个间隔插入4

12、个值 每间隔插入4个值波形及均值，均方值，方差 程序如下： t4=zeros(1,(length(t2)-1)*4+length(t2);%插值信号横坐标for k=1:length(t4) if mod(k,5)=1 t4(k)=(1/fs1)*fix(k/5); end if mod(k,5)=2 t4(k)=(1/fs1)*fix(k/5)+(1/fs1)/5; end if mod(k,5)=3 t4(k)=(1/fs1)*fix(k/5)+(2/fs1)/5; end if mod(k,5)=4 t4(k)=(1/fs1)*fix(k/5)+(3/fs1)/5; end if mod

13、(k,5)=0 t4(k)=(1/fs1)*(k/5-1)+(4/fs1)/5; end endinter4=interp1(t2,sample1,t4,cubic);figure(23)stem(t4,inter4);%插值信号xlabel(t),ylabel(y);title(4插值信号);figure(80)plot(t4,inter4);%插值信号xlabel(t),ylabel(y);title(4插值信号);disp(4插值信号均值:);inter4_mean=mean(inter4);%均值disp(4插值信号方差:);inter4_var=var(inter4);%方差disp

14、(4插值信号均方值:);inter4_rms=inter4_var+inter4_mean.2;%均方值均值：-0.0067均方值：0.5573方差：0.5573插值信号为： 每间隔插入4个值波形的相关函数 每间隔插入4个值波形的概率密度每间隔插入4个值波形的频谱f1=5*fs1*(0:length(t4)-1)/length(t4);%频域横坐标inter4_spectra=fft(inter4);figure(26);plot(f1(1:300),abs(inter4_spectra(1:300);%频谱xlabel(f),ylabel(Y);title(4插值信号频谱);频谱波形为： 每

15、间隔插入4个值波形的功率谱密度P5=inter4_spectra.*conj(inter4_spectra)/length(t4);figure(27);plot(f1(1:300),P5(1:300);%功率谱密度xlabel(f),ylabel(S);title(4插值信号功率谱密度);（2）每一个间隔插入8个值 每间隔插入8个值波形及均值，均方值，方差 均值：-0.0067均方值：0.5572方差：0.5573 每间隔插入8个值波形的相关函数每间隔插入8个值波形的概率密度每间隔插入8个值波形的频谱每间隔插入8个值波形的功率谱密度五：实验结果分析(1) 针对问题（6）， 对采样前后、插值前

16、插值后信号进行比较。观察在采样频率不变的情况下，信号频谱的变化和频谱的周期延拓性。问题分析由随机信号和DSP相关知识知道，对于时间连续但幅度不连续的信号，其DTFT是以采样频率为周期延拓的连续信号，幅度衰减为被采样信号的1/T。插值时则相反。 因此对采样前后的信号，根据上述结论可得，在采样后信号的频谱是以采样角频率为周期延拓的。由于信号最初经过40000hz的采样处理,经过滤波器后对输入信号进行了8000HZ的采样，在8000HZ的基础上在进行4间隔采样，8间隔采样。此处8000HZ处理信号为采样前信号，因此4间隔采样后的信号频谱减小为采样前信号的频谱周期的1/5，其幅度也减小为采样前信号的1

17、/5。同理对8间隔采样后的信号频谱减小为采样前信号的频谱周期的1/9，其幅度也减小为采样前信号的1/9，对插4点后的信号频谱增为插值前信号的频谱周期的5倍，其幅度也增大为插值前信号的5倍。对插8点后的信号频谱增为插值前信号的频谱周期的9倍，其幅度也增大为插值前信号的9倍。实际图像分析：由图像分析可见其完全符合理论分析。 讨论X(n)的自相关函数、功率谱密度与X(t)的自相关函数、功率谱密度之间的关系。问题分析由DSP相关知识知道，对于时间连续但幅度不连续的信号，其DTFT是以采样频率为周期延拓的连续信号，幅度衰减为被采样信号的1/T。插值时则相反。并且自相关函数和功率谱密度互为一对傅里叶变换对

18、。 因此对采样前后的信号，根据上述结论可得， 4间隔采样后的信号频谱减小为采样前信号的频谱的1/5，由于功率谱密度是频谱平方的时间平均，其幅度减小为采样前信号的1/25。同理对8间隔采样后的信号频谱减小为采样前信号的1/9，其功率谱密度减小为采样前信号的频谱周期的1/81，对插4点后的信号频谱增为插值前信号的频谱的5倍，其幅度也增大为插值前信号的25倍。对插8点后的信号频谱增为插值前信号的频谱的9倍，其幅度也增大为插值前信号的81倍。 对于自相关函数，由于其是频谱的傅里叶反变换，对于4采样信号的自相关函数其幅度是采样前的1/5，对于8采样信号的自相关函数其幅度是采样前的1/9。同样对于4和8插

19、值信号，其自相关函数幅度分别为采样前的5倍和9倍实际图像分析：由图像分析可见其完全符合理论分析。六源程序%生成输入信号Fs=40000;Ns=4096;n=0:Ns-1;t=n/Fs;%用离散序列模拟的正弦信号两点间距取1/40000，共取N=5000个采样点,n为长度为N的序列，t为采样时间，作为信号图像横坐标sine=sin(2*pi*100*t);%要求的正弦信号noise=randn(1,length(t);%高斯白噪声x=sine+noise; %合成输入信号figure(1);plot(t,x);%输入信号图形xlabel(t),ylabel(x);title(输入信号);%dis

20、p(输入信号均值:);x_mean=mean(x)%均值disp(输入信号方差:);x_var=var(x)%方差disp(输入信号均方值:)input_fangcha=x_var+x_mean.2%均方值%求相关函数和概率密度以及频谱x_xcorr=xcorr(x);t1=(-Ns+1:Ns-1)/Fs;%相关函数图像横坐标figure(2);plot(t1,x_xcorr);%相关函数xlabel(t),ylabel(R);title(输入信号相关函数);F1,y1=ksdensity(x);%y1为输入信号所有取值，F1为每个值相应概率figure(3);plot(y1,F1);%概率密

21、度xlabel(y),ylabel(F);title(输入信号概率密度);x_spectra=fft(x);f=Fs*n/Ns;%频域横坐标figure(4);plot(f(1:300),abs(x_spectra(1:300);%频谱xlabel(f),ylabel(Y);title(输入信号频谱);P1=x_spectra.*conj(x_spectra)/Ns;%频谱函数的平方求时间平均figure(5);plot(f(1:300),abs(P1(1:300);%功率谱密度xlabel(f),ylabel(S);title(输入信号功率谱密度);%hh%生成数字滤波器%巴特沃斯模拟滤波器

22、生成数字滤波器fp=1000;fs=2000;rp=1;rs=35;Fs1=40000;wp=2*pi*fp/Fs1;ws=2*pi*fs/Fs1;wap=tan(wp/2);was=tan(ws/2);Fs1=Fs1/Fs1;N,Wn=buttord(wap,was,rp,rs,s);z,p,k=buttap(N);bp,ap=zp2tf(z,p,k);%得到传输函数bs,as=lp2lp(bp,ap,wap);%低通到低通，频谱变换bz,az=bilinear(bs,as,Fs1/2); %将模拟滤波器传递函数转换为数字滤波器传递函数H,w=freqz(bz,az,256,Fs1*1000

23、0);%特性分析figure(6)plot(w,abs(H);%画出数字滤波器的波形图title(低通滤波器的频谱);xlabel(f/hz);grid on;%信号通过低通滤波器y=filter(bz,az,x);y_junzhi=mean(y);%均值y_fangcha=std(y);%方差y_junfangzhi=y_fangcha+y_junzhi.2;%均方值figure(7)plot(t,y);title(通过低通滤波器的波形);%求相关函数y_xcorr=xcorr(y);t1=(-Ns+1:Ns-1)/40000;%相关函数图像横坐标figure(8);plot(t1,y_xc

24、orr);%相关函数xlabel(t),ylabel(R);title(输出信号相关函数);%求概率密度F2,y2=ksdensity(y);%y1为输入信号所有取值，F1为每个值相应概率figure(9);plot(y2,F2);%概率密度xlabel(y),ylabel(F);title(输出信号概率密度);%求频谱y_spectra=fft(y);f=40000*n/Ns;%频域横坐标figure(10);plot(f(1:300),abs(y_spectra(1:300);%频谱%plot(abs(y_spectra);xlabel(f),ylabel(Y);title(输出信号频谱)

25、;%求功率谱密度P2=y_spectra.*conj(y_spectra)/Ns;%频谱函数的平方求时间平均figure(11);plot(f(1:300),abs(P2(1:300);%功率谱密度xlabel(f),ylabel(S);title(输出信号功率谱密度);%fs1=8000;%采样周期取8000Hzt2=(0:Fs/fs1:Ns-1)/Fs;%采样信号横坐标sample=zeros(1,length(t);for k=1:length(t) if mod(k,(Fs/fs1)=1 sample(k)=y(k); else sample(k)=0; endend%未除去相邻采样点

26、间零点的采样信号，为计算频谱和功率谱准备%sample1=zeros(1,ceil(length(t)/(Fs/fs1);i=1;for k=1:length(t) if mod(k,(Fs/fs1)=1 sample1(i)=y(k); i=i+1; endendfigure(12);stem(t2,sample1);%除去相邻采样点间零点的采样信号xlabel(t),ylabel(y);title(采样信号);%k_spectra=fft(sample);f=Fs*n/Ns;%频域横坐标figure(1221);plot(f(1:800),abs(k_spectra(1:800);%频谱x

27、label(f),ylabel(Y);title(8000信号频谱);%k_xcorr=xcorr(sample);t1=(-Ns+1:Ns-1)/Fs;%相关函数图像横坐标figure(11111);plot(t1,k_xcorr);%相关函数xlabel(t),ylabel(R);title(8000信号相关函数);%q1=k_spectra.*conj(k_spectra)/Ns;%频谱函数的平方求时间平均figure(22222);plot(f(1:1000),abs(q1(1:1000);%功率谱密度xlabel(f),ylabel(S);title(8000信号功率谱密度);%间隔

28、四点抽样sample2=zeros(1,ceil(length(t)/(Fs/fs1);for k=1:length(t2) if mod(k,5)=1 sample2(k)=sample(k); else sample2(k)=0; endendsample3=zeros(1,ceil(length(t)/(Fs/fs1)/5);i=1;for k=1:5:length(t2) sample3(i)=sample1(k); i=i+1;end figure(13);stem(1:length(t2)/5)/300,sample3);%除去相邻采样点间零点的采样信号xlabel(t),ylabel(y);title(4采样信号);%sample_junzhi4=mean(sample2)%均值sample_fangcha4=std(sample2)%方差sample_junfangzhi4=sample_fangcha4+sample_junzhi4.2%均方值%sample_xcorr4=x