数字信号处理第一次大作业

实验名称

实验1 基于谐波检测的移频闭塞信号数字接收

实验目的

1.掌握基于FFT 的信号频谱分析技术，认识近似分析中出现的混叠现象、泄露现象

和栅栏现象，加深理解这些现象对频谱分析精度的影响。

2.初步掌握噪声中谐波检测算法，了解影响频率估计精度的主要因素。

3.了解铁路移频闭塞系统的信号发送和接收过程。

实验内容

1.通过理论分析推导移频闭塞信号s(t)及其等效低通形式s l(t)的傅里叶变换表达式（为

离散谱），由此计算f1取不同值时国产18信息和ZPW-2000两类闭塞信号的主谐波分量频率和有效值，并统计其功率最强的几根谱线（国产18信息为6根，ZPW-2000 为3根）的功率之和占信号总功率的百分比。将计算结果制作成表格。

2.参照实验原理部分提供的波形仿真示例程序，针对国产18信息和ZPW-2000两种制

式，产生f1为标准值时移频闭塞信号等效低通信号仿真波形（可根据需要设定生成数据的长度和对采样频率等参数进行调整），由FFT分析该信号频谱并计算其主谐波的频率和有效值，将计算结果与内容1计算结果进行比较。

3.假设移频闭塞信号的参数f c和f1取标准值（即频率偏差为0），设计算法分析接收到

的等效低通信号，通过提取移频闭塞信号的主谐波分量，估计调制低频f1的最佳取值。将编写的程序对10段信号样本进行分析，输出f1的估计结果。

4.假设移频闭塞信号的参数f c和f1的偏差满足|Δf c|<5Hz，|Δf1|<0.1Hz，根据接收到的

等效低通信号设计算法估计Δf c和Δf1，并确定调制低频f1对应的最佳标准值。将编写的程序对10 段信号样本进行分析，输出Δf c、Δf1和f1的估计结果。

实验结果与分析

内容1

记

v=Δf

f1

=TΔf

由题可知

s(t)=A0∑{c n cosθ0cos[2π(f c+nf1)t]+c n sinθ0sin[2π(f c+nf1)t]} +∞

n=?∞

其中

c n=sin

π(v?n)

2

π(v?n)

+(?1)n

sin

π(v+n)

2

π(v+n)

若要求移频信号归一化功率为1，则移频信号的各个谐波分量的归一化值为

A?n=

c

√∑|c n|2

+∞

n=?∞

=

c

√|c0|2+2∑|c n|2

+∞

n=1

内容2

ZPW-2000制式，仿真结果如下

主谐波频率

f=1

T

=

1

(0.4937?0.009338)/13

≈26.8394Hz

偏移4Hz。归一化后，在4Hz±0.5f上进行积分后开方，得有效值0.6458。

国产18信息制式，仿真结果如下

f=1

T

=

1

(0.4982?0.01508)/8

≈16.5590Hz

主谐波频率为49.677Hz

偏移4Hz。归一化后，进行数值积分后开方，得有效值0.8177。

将以上结果与理论值相比较，可以看出，由于噪声作用，二者有一定误差，但仍然体现出了一致性。

内容3

对于ZPW-2000信息制式，幅频特性类似于下图

对于国产18制式，幅频特性类似于下图

从文件读入时域波形，用fft函数得到频谱。

取峰值中最大的n个（对于ZPW-2000信息制式，取n=3；对于国产18制式，取n=6），将其对应的频率放在数组f里。

以0.1为步长，不断列举f1的值，记录下使

n

S=∑|f1?f[i]|

i=1

最小的f1的值。

然后将f1与标准值相比较，选择最接近的一个。

编写程序如下

%初始化

close all;

clear all;

fs=1024;

f1=zeros(10,1);

type=6;

if (type==3)

F1Set=10.3+1.1*(0:17);

else

F1Set=[7,8,8.5,9.0,9.5,11.0,12.5,13.5,15.0,16.5,17.5,18.5,20.0,21.5,22.5,23.5,24.5,26];

End

%数据处理部分

for step=1:10

%读取数据

load 2009010990_dat3;

x=Data(:,step);

%频谱分析

X=fftshift(fft(x));

%抽取主谐波频率

[X_sorted,index_sorted]=sort(X);

n=length(index_sorted);k=type;m=0;t=zeros(type,1); while (m

a=sort(index_sorted(n-k+1:n));

m=1;

t(1)=a(1);

for i=2:k

if (a(i)~=a(i-1)+1)

m=m+1;

t(m)=a(i);

else

if (X(t(m))

t(m)=a(i);

end

end

end

k=k+1;

end

index=(t-257)*2;

%确定f1

b=sort(abs(index));

r=zeros(type,1);

delta=0;

while (delta<100)

for i=floor(min(F1Set))*10:b(round(type/2))*10 tmp=0;

for j=1:type

r(j)=b(j)/(i*0.1);

tmp=tmp+abs(r(j)-round(r(j)));

end

if (tmp<=0.02*delta)

f1(step)=i*0.1;

end

end

if (f1(step)<7)

delta=delta+1;

else

delta=100;

end

end

%将f1与标准值匹配

mindelta=abs(f1(step)-F1Set(1));

f1tmp=F1Set(1);

for i=2:length(F1Set)

if (mindelta>abs(f1(step)-F1Set(i)))

mindelta=abs(f1(step)-F1Set(i));

f1tmp=F1Set(i);

end

end

f1(step)=f1tmp;

end

%输出结果

disp(f1);

分别从2009010990_dat1和2009010990_dat3读取数据，运行程序进行分析，得到结果如下表所示。

与任务3算法原理相似，但需要加入对△f c和△f1的处理，因此程序结构上有所改变。

仍然从文件中读入采样数据，用fft()函数进行频谱分析，并抽取频谱中的n个最大值，将其对应频率放入数组f。

先以0.01为步长，不断列举△f c的值，再以0.01为步长，不断列举△f1的值，并枚举f1为标准值，记录下使

n

S=∑|f1?f[i]+△f c|

i=1

最小的△f c，△f1与f1。

编写程序如下。

close all;

clear all;

fs=1024;

type=6;

if (type==3)

F1Set=10.3+1.1*(0:17);

else

F1Set=[7,8,8.5,9.0,9.5,11.0,12.5,13.5,15.0,16.5,17.5,18.5,20.0,21.5,22.5,23.5,24.5,26];

end

for step=1:10

load 2009010990_dat4;

x=Data(:,step);

X=fftshift(fft(x));

[X_sorted,index_sorted]=sort(X);

n=length(index_sorted);k=type;m=0;t=zeros(type,1);

while (m

a=sort(index_sorted(n-k+1:n));

m=1;

t(1)=a(1);

for i=2:k

if (a(i)~=a(i-1)+1)

m=m+1;

t(m)=a(i);

else

if (X(t(m))

t(m)=a(i);

end

end

end

k=k+1;

end

index=(t-257)*2;

a=sort(abs(index));

df1=(-0.1:0.01:0.1)';

F1Set0=df1*ones(1,length(F1Set))+ones(length(df1),1)*F1Set;

deltamin=100;f1index=0;dffinal=0;

for dfc=-5:0.01:5

index0=index-dfc;

thedelta=100;

for i=1:length(F1Set)

for j=1:length(df1)

delta=sum(abs((index0/F1Set0(j,i))-round(index0/F1Set0(j,i))));

if (delta<=thedelta)

thedelta=delta;

theindex=i;

thedf1=df1(j);

end

end

end

if (thedelta<=deltamin)

deltamin=thedelta;

f1index=theindex;

dfcfinal=dfc;

df1final=thedf1;

end

end

disp([dfcfinal,df1final,F1Set(f1index)]);

end

分别从2009010990_dat2和2009010990_dat4读取数据，运行程序进行分析，得到结果

其f1应为30左右，远超出了标准值的取值范围，因此程序识别错误，认为f1=15Hz。我查看了其他同学的一些数据，没有发现类似这样的数据。

事实上在任务3中，制式为国产18信息的第9条数据，同样出现了f1=30Hz，只不过任务3中采用了另一种算法，程序没有出错。

误差分析

由于频谱分析后频率最小间隔为2Hz，因此对于单个峰值，误差在1Hz之内。而谐波分析过程中，采用了多个峰值所对应的频率，造成误差应当在0.1Hz数量级。需要注意的是，这一误差可以通过增加采样频率来减少。

忽略频谱离散造成的误差，该算法中△f c，△f1，f1的误差都在0.01Hz以内，并可以通过缩减步长而减小。