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

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

《完整word版平稳随机过程的采样和插值.docx》由会员分享，可在线阅读，更多相关《完整word版平稳随机过程的采样和插值.docx（52页珍藏版）》请在冰豆网上搜索。

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

随机信号实验

平稳随机过程的采样和插值

.实验目的

了解确定信号的采样与平稳随机信号的采样之间的关系，掌握信号的采样及分析方法。

实验原理

确定信号的采样符合香农定理，那么随机信号的采样是否符合香农定理呢？

答案是定的。

香农定理可以推广到随机信号的采样。

若x（t）为平稳随机过程，且具有零均值，它的功率谱密度sx（■）限于（-•飞，+'C）之

间。

当满足条件

X（t）专im'X（nT）sin（tct「n）

n—严ncoct—n兀

上式就是平稳随机过程的采样定理。

式中T为采样周期。

三.实验任务与要求

⑴程序用matlab或C/C++语言编写和仿真。

系统框图如图29、图30所示:

图29抽样系统框图

图30插值系统框图

⑵输入信号x（t）:

x（t）=正弦波信号+n（t）,频率为100Hz的正弦波信号，幅值为1v,n（t）为白噪声。

计算输入信号的均值、均方值、方差、相关函数、概率密度、频谱及功率谱密度、相关函数。

⑶低通滤波器设计

低通滤波器技术要求：

通带截止频率1KHz

阻带截止频率2KHz=

过渡带：

1KHz阻带衰减：

＞35DB

通带衰减：

＜1DB米样频率：

W44.1KHZ

计算经低通滤波器后信号的均值、均方值、方差、相关函数、频谱及功率谱密度、

相关函数。

⑷对输入信号进行抽样：

采样频率8000Hz。

每间隔4个点和每间隔8个点各抽样一

次。

计算抽样信号的均值、均方值、方差、相关函数、频谱及功率谱密度、相关函数。

⑸对采样信号进行插值：

每一个间隔插入4个值和每一个间隔插入8个值。

采样频率

8000Hz。

计算插值信号的均值、均方值、方差、相关函数、频谱及功率谱密度、相关函数。

⑹对采样前后、插值前插值后信号进行比较。

观察在采样频率不变的情况下，信号频

谱的变化和频谱的周期延拓性。

⑺讨论X（n）的自相关函数、功率谱密度与X（t）的自相关函数、功率谱密度之间的关系。

四•实验步骤及结果

1.输入信号x（t）

（1）x（t）=正弦波信号+n（t），频率为为高斯白噪声。

Fs=40000;Ns=5000;n=0:

Ns-1;t=n/Fs;

sine=sin（2*pi*100*t）;

noise=randn（1,length（t））;

x=sine+noise;

figure

（1）;

plot（t,x）;

xlabel（'t'）,ylabel（'x'）;

title（'输入信号’）；

100Hz的正弦波信号，幅值为1v,其中n（t）

两点间距取1/40000，共取N=5000个采样

点,n为长度为N的序列，t为采样时间，作为信号图像横坐标

产生频率为100hz的正弦波

产生高斯白噪声

信号合成为输入信号

产生的输入信号如下:

输入信号

0.14

5[I[[[l

00.020.040.060.080.10.12

t

（2）计算输入信号的均值、均方值、方差

均值：

x_mean=0.0072

均方值：

x_var=1.4747

方差：

input_fangcha=1.4748

（3）x（t）的相关函数利用MALTLA中的xcorr函数x_xcorr=xcorr（x）;

t1=（-Ns+1:

Ns-1）/Fs;%相关函数图像横坐标

figure

（2）;

plot（t1,x_xcorr）;%相关函数

xlabel（'t'）,ylabel（'R'）;

title（'输入信号相关函数’）；

t

（4）x（t）的概率密度利用ksdensity函数

[F1,y1]=ksdensity（x）;%y1为输入信号所有取值，F1为每个值相应概率

figure（3）;

plot（y1,F1）;%概率密度

xlabel（'y'）,ylabel（'F'）;

title（'输入信号概率密度'）;

0.35

0.3

0.25

0.2

F

0.15

0.1

0.05

0

-6-4-20246

y

（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（'输入信号频谱'）;

1600

1400

1200

1000

Y800

600

400

200

0-

05001000150020002500

f

（6）x（t）的功率谱密度用功率谱的概念求，即频谱函数的平方求时间平均。

P仁x_spectra.*conj（x_spect⑻/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.14

均值：

0.0416

均方值：

0.7248

方差：

0.7265

（2）输出y的相关函数

程序类似求x（t）的相关函数，其波形如下:

（3）y的概率密度

程序类似求x（t）的概率密度，其波形如下：

输出信号概率密度

y

2500

（4）

y的频谱

程序类似求x（t）的频谱，其波形如下：

1600

1400

1200

1000

丫800

600

400

200

0

0500100015002000

f

（5）y的功率谱密度

900

800

700

600

500

400

300

200

100

500

输岀信号功率谱密度

1000

1500

2000

2500

3000

4•对y信号进行抽样采样频率8000Hz

若8000hz采样，波形输出程序为：

sample仁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.5

00.020.040.060.080.1

t

0.12

（1）

每间隔4个点抽样一次

14采样波形及均值，均方值，方差

程序代码如下：

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」unzhi4=mean（sample2）%均值

sample_fangcha4=std（sample2）%方差

sample」unfangzhi4=sample_fangcha4+sample_junzhi4.A2%均方值

4采样信号为：

均值：

-0.0125

方差:

0.3212

均方值：

0.3213

24采样信号的相关函数代码函数如上

③4采样信号的概率密度

y

④4采样信号的频谱

4米样信号频谱

⑤4采样信号的功率谱密度

4采样信号功率谱密度

（2）每间隔8个点抽样一次程序和4采样类似，只是间隔8个点

①采样信号

8米样信号

t

均值：

-7.0098e-004

均方值：

0.1099

方差：

0.1099

2

8采样信号的相关函数

10

8

6

4

2

R

0

-2

-4

-6

-8

-1000-800-600-400-20002004006008001000

t

3

8采样信号的概率密度

0.7

0.6

0.5

0.4

F

0.3

0.2

0.1

-4-3-2-101234

y

48采样信号的频谱

10

9

8

7

6

丫5

4

3

2

1

0

8米样信号频谱

05001000

f

58采样信号的功率谱密度

8采样信号功率谱密度

f

4．对信号进行插值采样频率为8000Hz

插值思路：

用三次多项式插值法

（1）每一个间隔插入4个值

1每间隔插入4个值波形及均值，均方值，方差程序如下：

t4=zeros（1,（length（t2）-1）*4+length（t2））;%插值信号横坐标fork=1:

length（t4）

ifmod（k,5）==1t4（k）=（1/fs1）*fix（k/5）;

end

ifmod（k,5）==2t4（k）=（1/fs1）*fix（k/5）+（1/fs1）/5;

end

ifmod（k,5）==3t4（k）=（1/fs1）*fix（k/5）+（2/fs1）/5;

end

ifmod（k,5）==4t4（k）=（1/fs1）*fix（k/5）+（3/fs1）/5;

end

ifmod（k,5）==0t4（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（'4插值信号均方值:

'）;

inter4_rms=inter4_var+inter4_mean42;%均方值

均值：

-0.0067

均方值：

0.5573

方差：

0.5573

插值信号为:

y

4插值信号

t

②每间隔插入4个值波形的相关函数

4插值信号相关函数

-0.2-0.15-0.1-0.0500.050.10.15

t

3每间隔插入4个值波形的概率密度

y

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插值信号频谱'）;

频谱波形为：

2000

1800

1600

1400

1200

1000

800

600

400

200

0

5每间隔插入4个值波形的功率谱密度

P5=inter4_spectra.*conj（inter4_spectra）/length（t4）;

figure（27）;

plot（f1（1:

300）,P5（1:

300））;%功率谱密度xlabel（'f）,ylabel（'S'）;

title（'4插值信号功率谱密度'）;

900

800

700

600

500

400

300

200

100

（2）每一个间隔插入8个值

①每间隔插入8个值波形及均值，均方值，方差

8插值信号

均值：

-0.0067

均方值：

0.5572

方差：

0.5573

2每间隔插入8个值波形的相关函数

t

3

每间隔插入8个值波形的概率密度

4每间隔插入8个值波形的频谱

8米样信号频谱

10-

~~

1■

9-

一

8-

II'

7-

1

一

11

6

||

1J

!

■\\!

；1

丫5-

-

111/

:

1

IU

li；flJ

J1-.1

4亠

'1|

r1

I

1

1

1"||

I*11=1

II1

3-

1

1

||

1-

11]11

11

11

;1

2-

|1

1

11

11

MHI1

1■

li/

JJ讣

1t

j/1*

\V1,/'?

.1

\9V,if"—

\

i1

E

/1J\

\\\i

0L

r

b

r>1

050010001500

f

5每间隔插入8个值波形的功率谱密度

8插值信号功率谱密度

1600

1400-

1200-

1000-

S800-

600「

400-

200、

0

0

五：

实验结果分析

（1）针对问题（6）,对采样前后、插值前插值后信号进行比较。

观察在采样频率不变的情况

下，信号频谱的变化和频谱的周期延拓性。

1问题分析

由随机信号和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倍。

2实际图像分析：

3

由图像分析可见其完全符合理论分析。

⑺讨论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.A2%均方值

%%

%求相关函数和概率密度以及频谱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.A2;%均方值

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;%频谱函数的平方求时间平均f