完整word版平稳随机过程的采样和插值.docx
《完整word版平稳随机过程的采样和插值.docx》由会员分享，可在线阅读，更多相关《完整word版平稳随机过程的采样和插值.docx（44页珍藏版）》请在冰豆网上搜索。
完整word版平稳随机过程的采样和插值
随机信号实验
平稳随机过程的采样和插值
一．实验目的
了解确定信号的采样与平稳随机信号的采样之间的关系，掌握信号的采样及分析方法。
二.实验原理
确定信号的采样符合香农定理，那么随机信号的采样是否符合香农定理呢？
答案是定的。
香农定理可以推广到随机信号的采样。
若X（t）为平稳随机过程，且具有零均值，它的功率谱密度
限于（-
，+
）之间。
当满足条件
时，便可将X（t）按它的振幅样本展开为：
上式就是平稳随机过程的采样定理。
式中T为采样周期。
三．实验任务与要求
⑴程序用matlab或c/c++语言编写和仿真。
系统框图如图29、图30所示：
图29抽样系统框图
图30插值系统框图
⑵输入信号x（t）：
x（t）=正弦波信号+n（t），频率为100Hz的正弦波信号，幅值为1v，n（t）为白噪声。
计算输入信号的均值、均方值、方差、相关函数、概率密度、频谱及功率谱密度、相关函数。
⑶低通滤波器设计
低通滤波器技术要求：
通带截止频率1KHz
阻带截止频率2KHz。
过渡带:
1KHz
阻带衰减：
>35DB
通带衰减：
<1DB
采样频率：
≤44.1KHz
计算经低通滤波器后信号的均值、均方值、方差、相关函数、频谱及功率谱密度、相关函数。
⑷对输入信号进行抽样：
采样频率8000Hz。
每间隔4个点和每间隔8个点各抽样一次。
计算抽样信号的均值、均方值、方差、相关函数、频谱及功率谱密度、相关函数。
⑸对采样信号进行插值：
每一个间隔插入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为长度为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_xcorr=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;%频域横坐标
figure（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;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*10000）;%特性分析
figure（6）
plot（w,abs（H））;%画出数字滤波器的波形图
title（'低通滤波器的频谱'）;
xlabel（'f/hz'）;
gridon;
3．x（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采样，波形输出程序为：
sample1=zeros（1,ceil（length（t）/（Fs/fs1）））;
i=1;
fork=1:
length（t）
ifmod（k,（Fs/fs1））==1
sample1（i）=y（k）;
i=i+1;
end
end
figure（12）;
stem（t2,sample1）;%除去相邻采样点间零点的采样信号
xlabel（'t'）,ylabel（'y'）;
title（'采样信号'）;
则波形为：
（1）每间隔4个点抽样一次
4采样波形及均值，均方值，方差
程序代码如下：
sample2=zeros（1,ceil（length（t）/（Fs/fs1）））;
fork=1:
length（t2）
ifmod（k,5）==1
sample2（k）=sample（k）;
else
sample2（k）=0;
end
end
sample3=zeros（1,ceil（length（t）/（Fs/fs1）/5））;
i=1;
fork=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%均方值
4采样信号为：
均值：
-0.0125
方差：
0.3212
均方值：
0.3213
4采样信号的相关函数代码函数如上
4采样信号的概率密度
4采样信号的频谱
4采样信号的功率谱密度
（2）每间隔8个点抽样一次程序和4采样类似，只是间隔8个点
采样信号
均值：
-7.0098e-004
均方值：
0.1099
方差：
0.1099
8采样信号的相关函数
8采样信号的概率密度
8采样信号的频谱
8采样信号的功率谱密度
4．对信号进行插值采样频率为8000Hz
插值思路：
用三次多项式插值法
（1）每一个间隔插入4个值
每间隔插入4个值波形及均值，均方值，方差
程序如下：
t4=zeros（1,（length（t2）-1）*4+length（t2））;%插值信号横坐标
fork=1:
length（t4）
ifmod（k,5）==1
t4（k）=（1/fs1）*fix（k/5）;
end
ifmod（k,5）==2
t4（k）=（1/fs1）*fix（k/5）+（1/fs1）/5;
end
ifmod（k,5）==3
t4（k）=（1/fs1）*fix（k/5）+（2/fs1）/5;
end
ifmod（k,5）==4
t4（k）=（1/fs1）*fix（k/5）+（3/fs1）/5;
end
ifmod（k,5）==0
t4（k）=（1/fs1）*（k/5-1）+（4/fs1）/5;
end
end
inter4=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（'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插值信号频谱'）;
频谱波形为：
每间隔插入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），对采样前后、插值前插值后信号进行比较。
观察在采样频率不变的情况下，信号频谱的变化和频谱的周期延拓性。
问题分析
由随机信号和DSP相关知识知道，对于时间连续但幅度不连续的信号，其DTFT是以采样频率为周期延拓的连续信号，幅度衰减为被采样信号的1/T。
插值时则相反。
因此对采样前后的信号，根据上述结论可得，在采样后信号的频谱是以采样角频率为周期延拓的。
由于信号最初经过40000hz的采样处理,经过滤波器后对输入信号进行了8000HZ的采样，在8000HZ的基础上在进行4间隔采样，8间隔采样。
此处8000HZ处理信号为采样前信号，因此4间隔采样后的信号频谱减小为采样前信号的频谱周期的1/5，其幅度也减小为采样前信号的1/5。
同理对8间隔采样后的信号频谱减小为采样前信号的频谱周期的1/9，其幅度也减小为采样前信号的1/9，对插4点后的信号频谱增为插值前信号的频谱周期的5倍，其幅度也增大为插值前信号的5倍。
对插8点后的信号频谱增为插值前信号的频谱周期的9倍，其幅度也增大为插值前信号的9倍。
实际图像分析：
由图像分析可见其完全符合理论分析。
⑺讨论X（n）的自相关函数、功率谱密度与X（t）的自相关函数、功率谱密度之间的关系。
问题分析
由DSP相关知识知道，对于时间连续但幅度不连续的信号，其DTFT是以采样频率为周期延拓的连续信号，幅度衰减为被采样信号的1/T。
插值时则相反。
并且自相关函数和功率谱密度互为一对傅里叶变换对。
因此对采样前后的信号，根据上述结论可得，4间隔采样后的信号频谱减小为采样前信号的频谱的1/5，由于功率谱密度是频谱平方的时间平均，其幅度减小为采样前信号的1/25。
同理对8间隔采样后的信号频谱减小为采样前信号的1/9，其功率谱密度减小为采样前信号的频谱周期的1/81，对插4点后的信号频谱增为插值前信号的频谱的5倍，其幅度也增大为插值前信号的25倍。
对插8点后的信号频谱增为插值前信号的频谱的9倍，其幅度也增大为插值前信号的81倍。
对于自相关函数，由于其是频谱的傅里叶反变换，对于4采样信号的自相关函数其幅度是采样前的1/5，，对于8采样信号的自相关函数其幅度是采样前的1/9。
同样对于4和8插值信号，其自相关函数幅度分别为采样前的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（'输入信号'）;
%%
disp（'输入信号均值:
'）;
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）;%概率密度
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
%%
%生成数字滤波器
%巴特沃斯模拟滤波器生成数字滤波器
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*10000）;%特性分析
figure（6）
plot（w,abs（H））;%画出数字滤波器的波形图
title（'低通滤波器的频谱'）;
xlabel（'f/hz'）;
gridon;
%%
%信号通过低通滤波器
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_xcorr）;%相关函数
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（'输出信号频谱'）;
%%
%求功率谱密度
P2=y_spectra.*conj（y_spectra）/Ns;%频谱函数的平方求时间平均
figure（11）;
plot（f（1:
300）,abs（P2（1:
300）））;%功率谱密度
xlabel（'f'）,ylabel（'S'）;
title（'输出信号功率谱密度'）;
%%
%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%
fs1=8000;%采样周期取8000Hz
t2=（0:
Fs/fs1:
Ns-1）/Fs;%采样信号横坐标
sample=zeros（1,length（t））;
fork=1:
length（t）
ifmod（k,（Fs/fs1））==1
sample（k）=y（k）;
else
sample（k）=0;
end
end%未除去相邻采样点间零点的采样信号，为计算频谱和功率谱准备
%%
sample1=zeros（1,ceil（length（t）/（Fs/fs1）））;
i=1;
fork=1:
length（t）
ifmod（k,（Fs/fs1））==1
sample1（i）=y（k）;
i=i+1;
end
end
figure（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）））;%频谱
xlabel（'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信号功率谱密度'）;
%%
%间隔四点抽样
sample2=zeros（1,ceil（length（t）/（Fs/fs1）））;
fork=1:
length（t2）
ifmod（k,5）==1
sample2（k）=sample（k）;
else
sample2（k）=0;
end
end
sample3=zeros（1,ceil（length（t）/（Fs/fs1）/5））;
i=1;
fork=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
升级会员 