清华大学测试信号分析与处理实验一和实验二报告

测试信号分析与处理实验报告（实验一和实验二）

姓名：

学号：

日期：

实验一：利用FFT 作快速相关估计

1. 实验目的：

a. 掌握信号处理的一般方法，了解相关估计在信号分析与处理中的 作用。

b. 熟悉FFT 算法程序；熟练掌握用FFT 作快速相关估计的算法。

c. 了解快速相关估计的谱分布的情况。

2.实验内容：

a. 读入实验数据Qjt.dat 。

b. 编写一利用FFT 作相关估计的程序。

c. 将计算结果表示成图形的形式，给出相关谱的分布情况图。

3.算法讨论及分析：

遍历离散随机序列的自相关函数的理论计算式：

∑-=∞→+=1

)()(1lim )(N n N xx m n x n x N m R

可以用一个有限长度为N 的样本来估计自相关函数：

无偏估计：

)1(,......,1,0 )()(||1)(1

||0

-±±=+-=∑--=∧

N m m n x n x m N m R m N n xx 可以在数学上证明，当|m|值接近N 时，估计方差很大。因此实际往往采用下面的有偏估计。 有偏估计：

)1(,......,1,0 )()(1

)(1||0

-±±=+=

∑--=∧

N m m n x n x N

m R m N n xx

可以证明，当N 变大时，估计偏差渐渐趋于0，估计方差也趋于0.因此这是一个一致估计。

为了提高估计质量，需要增大N ，而增大N 将带来计算量的剧增。因此在算法上需要用到FFT 快速算法，步骤如下：

Step1：对原序列补N 个零，得到新序列)(2n x N Step2：作)()]([FFT 22ωj N N e X n x ⇒或)(2k X N

Step3：取)(2k X N 的共轭)(*

2k X N ，得

22|)(|1

k X N

N Step4:作)(|)(|1IFFT 22m R k X N o N ∧

⇒⎥⎦

⎤⎢⎣⎡ Step5：调整)(m R o ∧

与)(m R xx ∧

的错位

4.实验结果分析

4.1关于原信号是否可作为平稳信号的分析：

以dL=2^12=4096为单位长度，截取原信号前5个dL长的数据，分别求均值，结果如下：

Ex1 Ex2 Ex3 Ex4 Ex5 Ex

1.2708 0.0430 0.9709

2.0081 0.7288 0.8891

其中：Ex为原信号均值

相比于原信号的数据幅度|max(x)-min(x)|=5168,不同时间段的采样均值与信号总均值都很接近。因此，可以近似认为原信号为平稳信号。

4.2关于原信号的周期性和信噪比的分析：

周期性：

取采样序列数据长度N=512，求其自相关函数如下：

由上图可看出，自相关函数具有某种正负的周期性。由于自相关函数能保持原信号的周期性质不变，因此可以认为原信号具有类似的周期性。

信噪比：

对全部数据长度求FFT变换，得到下图：

原信号为桥梁振动信号，理论上应有n根幅值谱线分别对应不同模态的共振频率，并且每根谱线的幅值依次下降。由上图可见，幅值谱的谱线高度明显大于相邻两根谱线之间的幅值高度（即噪声产生的频谱谱线）。因此可以认为原信号的信噪比较高。

4.3用FFT对原信号进行自相关估计：

分别选择N=2048,4096,8192,16384,32768对信号进行FFT自相关估计：

各自相关估计的函数图如下：

从结果可以看出，当N变大时，自相关估计逐渐逼近原信号全长度的自相关函数（即估计方差变小，估计质量变高）。

4.4有偏估计与无偏估计的对比

有偏估计可以采用matlab自带函数Rx=xcorr（x,’biased’）来实现,也可以采用第3节（算法讨论与分析）中提到的有偏估计的算法来实现。

无偏估计采用matlab自带函数Rx=xcorr（x,’unbiased’）。

结果如下：

可见：尽管有偏估计是有偏差的，但是其方差比无偏估计好。当N变大时，有偏估计的偏差会变小，而无偏估计的方差会变得更差。因此综合来看，有偏估计要优于无偏估计。

5.原程序清单

clc

clear all

%读取数据点：

x=textread('E:\第一学期课程\测试信号分析与处理\仿真实验资料\Qjt.dat');

%========关于原信号是否可作为平稳信号的分析：=========

dL=4096;%分段长度

x1=x(1:dL);

x2=x(dL+1:2*dL);

x3=x(2*dL+1:3*dL);

x4=x(3*dL+1:4*dL);

x5=x(4*dL+1:5*dL);%取前5段数据

Ex1=mean(x1);

Ex2=mean(x2);

Ex3=mean(x3);

Ex4=mean(x4);

Ex5=mean(x5);%求各段均值

Ex=mean(x);%求原信号均值

%========关于原信号的周期性和信噪比的分析：=========

%=======周期性分析：=======

clear all

%读取数据点：

x=textread('E:\第一学期课程\测试信号分析与处理\仿真实验资料\Qjt.dat');

N=512;%采样点数

L=2*N;y=zeros(L);%将采样信号补零至2N项

for i=1:L

if i

y(i)=x(i);

else y(i)=0;

end

end

Y=fft(y,L);%对y进行FFT运算

P=(abs(Y).^2).*2/L;%求序列Y的功率谱密度函数

Ro=ifft(P,L);%对功率谱密度函数P求FFT逆变换

Rx=[Ro(N+2:L),Ro(1:N)];%对Ro移位得到Rx，Rx只有2N-1项

%显示自相关函数：

figure(1)

i=1:2*N-1;