测试信号分析与处理作业实验一二

王锋

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

一、实验目的

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

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

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

二、实验内容

a.读入实验数据[1]。

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

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

注[1]：实验数据文件名为“Qjt.dat ”。

实验数据来源：三峡前期工程 “覃家沱大桥” 实测桥梁振动数据。 实验数据采样频率：50Hz 。

可从数据文件中任意截取几段数据进行分析，数据长度N 自定。 注[2]：采用Matlab 编程。

三、算法讨论及分析

算法为有偏估计，利用FFT 计算相关函数 Step 1: 对原序列补N 个零，得新序列x 2N (n) Step2:作FFT[x 2N (n)]得到X 2N (k)

2,0,1,......,212k k k N N

πω==-

Step 3: 取X 2N (k)的共轭2()N X k *

，得221()N X k N

Step 4: 作220

1?()()N IFFT X k R m N ???????

Step 5: 调整0?()R m 与?()Rxx m 的错位。

四、实验结果分析

1.该信号可以近似为平稳信号么？

可以近似为平稳信号，随机过程的统计特性不随样本的采样时刻而发生变化。取N=8192，分别取间隔m=500,m=700,m=1000,所得到的均值均为0.5366，方差为47369，与时间无关。

图1-1 自相关函数图

(上图表示的R0，下图为调整后的R0)

2.该信号是否具有周期性，信噪比如何？

>> load Qjt.dat; %加载数据

N=32768; %数据长度

i=1:1:N; %提取数据

plot(i,Qjt(i));

抛去几个极值点，从图1-2可以看出，数据具有一定的周期性，杂音比较少，说明信噪比较高。

图1-2 数据图

3.数据量N 对估计结果有何影响？ 在实验中，用FFT 做相关估计

1

1?()()()N m xx

n R m x n x n m N

--==+∑

，

有偏估计均值的表达式为

?()()xx xx N m

E R m R m N

-??=??

有偏估计方差的表达式为

2

2

?()1()()()xx

xx xx xx R m i R i R i m R i m N σ

∞=-∞

??≈++-??∑ 可知，N 越大，方差越小，均值越趋于()xx R m ，?()xx

R m 为无偏渐进一致估计，估计效果好。

4.采用有偏估计和无偏估计的结果有何不同？哪一种结果更好？ 无偏估计的均值为

1

1?()()N m xx xx n E R m R m N m

--=??=??-∑

无偏估计的方差为

()2

2

?()()()xx

xx xx R m i N R i R i m N m σ

∞

=-∞

??≈++??

-∑

可以看出均值为无偏估计，但是当m增大时，方差大，估计值分散，综合质量差，有偏估计的方差就没有这样的问题，估计值比较集中，故有偏估计更适合实际应用。

五、原程序清单

>> %%计算数据的自相关函数的程序

load Qjt.dat; %加载数据

N=16384; %数据长度

%N=8192;

for i=1:1:N %提取数据

X(i)=Qjt(i);

end

num=2*N; %补零措施

for i=N+1:1:num

X(i)=0;

end

Xfft=fft(X); % x2N(n)进行FFT变换X2N(k)

conXfft=conj(Xfft);%求X2N(k)的共轭

R0=ifft((Xfft.*conXfft)/N);%求自相关函数

Rxx=0*R0;

for i1=1:1:N

Rxx(i1)=R0(i1+N);

Rxx(i1+N)=R0(i1);

end

n=1:1:num;

subplot(2,1,1);

stem(n,R0);%绘制R0的图

subplot(2,1,2);

axis on

stem(fftshift(R0));%调整后的R0

%%计算数据的均值和均方差的程序

load Qjt.dat; %加载数据

%N=16384;

N=8192; %选取数据量

for i=1:1:N

X(i)=Qjt(i); %选取数据，数据范围1~8192

end

old_mean=mean(X) %求数据的均值

old_var=var(X) %求数据的方差

%选取数据，与原数据的间隔为m=500,数据范围500~8692

for i=500:1:(N+500)

X(i)=Qjt(i);

end

new500_mean=mean(X)%求数据的均值

new500_var=var(X) %求数据的方差

%选取数据，与原数据的间隔为m=700,数据范围700~8892

for i=700:1:(N+700)

X(i)=Qjt(i);

end

new700_mean=mean(X)%求数据的均值

new700_var=var(X) %求数据的方差

%选取数据，与原数据的间隔为m=1000,数据范围1000~9192

for i=1000:1:(N+1000)

X(i)=Qjt(i);

end

new1000_mean=mean(X)%求数据的均值

new1000_var=var(X) %求数据的方差

六、实验后的体会与建议

通过该实验，熟练应用了FFT作快速相关估计的算法来解决实际中的问题，更加熟悉了MATLAB相关的操作。

建议：课上可以多教授一些关于MATLAB的相关函数的用法，以及常见的误区。

实验二：采用Welch法作谱估计

一、实验目的

a. 掌握We l c h算法的概念、应用及特点；

b. 了解谱估计在信号分析中的作用；

c. 能够利用We l c h法对信号作谱估计，对信号的特点加以分析。

二、实验内容

a. 读入实验数据。

b. 编写一利用We l c h法作谱估计的算法程序。

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

三、算法讨论及分析

四、实验结果及分析

图2-1 实验结果图

1 不同的分段点数对Welch谱估计结果有何影响？

分段数K增大，每组数据量M减小，则会造成谱估计的偏差增大，曲线平滑；反之, 曲线起伏激烈, 谱线的尖峰保留。从图2-2中的K值大于图2-3的K 值，图2-2中的曲线比图2-3的曲线平滑，但是图2-3的尖峰效果比图2-2好。

图2-2 K =N/1100时的频谱图

图2-3 K=N/3000时的频谱图

2 不同的数据重叠长度Welch 对谱估计结果有何影响？

如果数据交叠使用，分段数K 比平均周期图法增大，则估计方差比平均周期图法小，数据交叠长度越大，分段数K 越大，则估计方差会越小，图形会变得平滑，如图2-4所示，上图数据重叠度大，K 大，则曲线相对平滑。从式（1）也可以看出，式（1）为谱估计的方差。

()()2

4sin ?var 1sin x xx M S K M ωσωω????????=+ ??????

???（1）

图2-4 使用不同的重叠度得到的谱估计图

（上：重叠度为97%，下：重叠度为3%）

3 不同的窗函数对Welch 谱估计结果有何影响？

使用窗函数，会使Welch 谱估计的图形变得平滑，但是，不同的窗函数不同，则其主瓣宽度不同。主瓣宽度越窄，对信号的主瓣频率辨识效果越好，分辨率较高。从图2-5可以看出，使用矩形窗的分辨率高于汉宁窗。

图2-5 使用不同的窗函数得到的谱估计

（上：矩形窗，下：汉宁窗）

4周期图法与Welch 法的谱估计结果有何不同？

周期图法本身的误差就很大，不是无偏估计，而且偏差比较大；welch 方法应该会准确一些，它频谱分辨率可能低一些，但估计的方差会小些，所以功率谱看起来平滑一些，另外一个原因是welch 方法你是加了窗的。两个估计都是跟准确值有偏差的，只是后者统计方差更小一些，相对更加可靠一点。

五、原程序清单

load Qjt.dat %加载数据

N=length(Qjt); %获得数据的长度

L=1800; %确定每段数据的长度

K=N/L; %在总数据中共取K段，没有重合

%K=2*N/L-1; %在总数据中共取K段，有重合

K=fix(K);

for num=1:1:K %使用循环，提取每段数据

for i=1:1:L

A(i,num)=Qjt((num-1)*L+i); %每段数据占据A阵的一列，没有重合

% A(i,num)=Qjt((num-1)*L/2+i); %每段数据占据A阵的一列，有重合

end

for i2=L+1:1:2048 % 使用循环，补零

A(i2,num)=0;

end

end

for num=1:1:K %使用循环，对每段数据使用矩形窗截取

tempA=A(:,num); %

tempw=boxcar(2048); % 获取窗函数

for i3=1:1:2048 %对每段数据使用矩形窗截取

A(i3,num)=tempA(i3)*tempw(i3);

end

end

for num=1:1:K %对截取之后的数据进行FFT变换，并获得傅里叶变换后的数据模

tempA=A(:,num);

tempFFT=fft(tempA); %对每一段数据进行FFT变换

for i3=1:1:2048 %获取数据FFT变换后的模

Sxi(i3,num)=tempFFT(i3)*conj(tempFFT(i3))/2048;

end

end %循环结束

U=dot(boxcar(2048),boxcar(2048))/2048; %获得修正系数U

tempSUM=0*ones(2048,1);

for num=1:1:K

tempSUM=Sxi(:,num)+tempSUM; %获得数据FFT变换后的功率谱和end

Sxx=tempSUM/(K*U); %调整系数

subplot(3,1,1);

Sxxi=ones(1024,1);

for inum=1:1:1024

Sxxi(inum)=Sxx(inum);

end

plot(Sxx); %绘制0~2048点的图

[Px,F]=pwelch(Qjt,boxcar(2048),0,2048); %使用pwelch函数计算

subplot(3,1,2);

plot(Sxxi);

subplot(3,1,3);

plot(Px); %绘制pwelch所得的图

六、实验后的体会和建议

通过此次实验，掌握了We l c h算法的概念、应用及特点，更加熟悉了MATLAB 相关的操作。

建议：课上可以多教授一些关于MATLAB的相关函数的用法，以及常见的误区。