DSP数字信号处理课程设计-高铁轨道电路移频信号数字检测研究与仿真（swjtu hanser_05）

一背景介绍

1.研究背景及其意义

铁路在国民经济的发展中重要性不言而喻，它对于促进物资交流，繁荣社会经济起着十分重要的作用。随着高速铁路  技术的发展，列车进行了多次提速，为了确保铁路的安全运行，要对列车运行主要是通过移动闭塞系统进行调控，来确保行车安全。我国的轨道交通信号技术已在吸收国外新技术的基础上不断的发展,并已在我国干线铁路及城市轨道交通中得到了应用，为了保证行车安全，确保铁路运输的安全畅通，也对列车在运行时的行车安全提出了更高的要求。

轨道电路是由发送器、两根钢轨导体、轨道绝缘和接收器通过一定方式的电气连接而成的电路。它不仅能够发送轨道区段的占用空闲信息和该段钢轨线路的完整性信息，而且还可以作为一个信息发送器，向列车传输控制信息，实现车地通信。轨道电路是现代铁路信号中自动控制和远程控制系统中的重要组成元件，是符合故障一安全原则的系统。

随着我国高速铁路发展计划的实施，铁路系统对列车自动控制的要求越来越高。列车控制系统必须提高对列车的控制精度才能保证行车安全,提高运输效率。列车的控制精度与轨道电路所提供的信息密切相关，轨道电路提供的参数越多、越精确，控制精度也就越高，运输效率也随之提高。由于轨道电路所处的环境非常恶劣，需要定期检测轨道电路的相关运行参数，可以及时消除隐患，确保列车的行车安全,所以移频信号准确、实时的检测对列车行车安全有十分重要的意义。

2.国内外的研究现状

ZPW-200OA无绝缘移频自动闭塞是在法国UM71无绝缘轨道电路技术引进并国产化的基础上，结合我们国家的国情，进行技术再开发。2003年铁路自动闭塞系统采用ZPW-2000系列轨道线路，为高速铁路建设奠定了基础。截至2018年底，ZPW-2000A系列轨道线路已投入运行近3万 km，其安全性、可靠性、环境适应性和可维护性都得到了显著提高，为10年来中国高速铁路的安全稳定、快速发展提供了有力支持。随着计算机技术、通信技术、自动化技术和网落技术的发展，铁路信号技术也将由开环控制向闭环控制发展，由孤立分散的控制向网落化、区域化控制发展，由单一的信号向通信一体化发展。区间的闭塞设备又将向新的模式发展，将由固定闭塞向移动闭塞发展。随着列车运行速度的不断提高，我国高速铁路、重轨货物运输、轨道交通异常发展，现有自动闭塞系统系统已远远不能满足列车高速行驶的需要，为此，必需加大力度要发展以计算机、通信、自动化和网落化技术为一体的数字轨道电路系统和列车运行自动控制系统（ATC），以满足铁路、轨交系统各列车安全运行的需要。

自动闭塞系统

轨道电路不同的闭塞区段对应着不同的低频信息，不同的低频信息代表着不同轨道信息。

3.研究思路与创新点

• • 二.信号产生

1.信号产生原理

ZPW-2000A型轨道电路采用的是相位连续的移频键控信号。这种移频信号的产生是采用频率调制的方法，即通过低频调制信号控制载频信号，形成振幅不变而频率随低频信号做周期性变化的信号。（通信原理2FSK信号调制，p108）

调制（modulation)就是对信号源的信息进行处理加到载波上，使其变为适合于信道传输的形式的过程，就是使载波随信号而改变的技术。

一般来说，信号源的信息（也称为信源）含有直流分量和频率较低的频率分量，称为基带信号。基带信号往往不能作为传输信号，因此必须把基带信号转变为一个相对基带频率而言频率非常高的信号以适合于信道传输。这个信号叫做已调信号，而基带信号叫做调制信号。调制是通过改变高频载波即消息的载体信号的幅度、相位或者频率，使其随着基带信号幅度的变化而变化来实现的。

.移频自动闭塞概念

移频自动闭塞是以轨道移频电路为基础的自动闭塞。移频自动闭塞选用频率参数作为控制信息，采用频率调制的方法，把低频信号搬移到较高频率上，从而形成振幅不变、频率随低频信号的幅度做周期变换的调频信号。将此信号利用钢轨作为传输通道来控制通过信号机的显示，达到自动指挥列车运行的目的。移频信号的波形如图所示：

                                                                   

移频信号的时域表达式为：

st=A0cosθt=A0cos⁡(ω0t+gt)

式中: A0 为移频信号的幅度; ω0为载频角频率。

gt=kf(t)dt=∆ωt                  -T4<t<T4∆ωT2-t      T4<t<3T4

式中 Δω 为移频信号的角频偏。

由 S(t)的傅里叶级数展开式可以得到：移频信号的频谱是以载频fc=2×π×ω0

为中心频率，以f1=2×π×nω1向两边扩散。

2.信号产生过程及其仿真

在低频调制方波信号 f(t)的控制下，载频信号的频率将以频偏∆f作偏移，最终形成以载频𝑓𝑐为中心,以𝑓ℎ = 𝑓𝑐 + ∆𝑓为上边频，以𝑓𝑙 = 𝑓𝑐 − ∆𝑓为下边频的移频信号。经过方波的调制，其角频率也发生了变化，具体的偏移量∆𝜔 = 𝑘𝑓(𝑡)，k 为常数，表示移频器的灵敏度。（运行代码见附页）

fc=20hz,fd=3hz，fa=5hz

根据其中一种低频调制频率和载频频率得出的图像如下：

fc=1701.4hz；           fd=15.8hz；

移频信号的时域波形和频域波形

fft变换后频谱图

ZPW-2000A中有18种低频信号，每一种低频信号都有不同的意义，若直接传输低频信号，低频信号容易被干扰，被衰弱，因此需要对低频信号进一步处理，再传输。

注：除已经定义的14个信息使用标准，未作定义：21.3Hz，23.5Hz，25.7Hz，27.9Hz

18种低频调制频率是：
10.3、11.4、12.5、13.6、14.7、15.8、16.9、18、19.1、20.2、21.3、22.4、23.5、24.6、25.7、26.8、27.9、29

低频信号的处理需要借助载频。载频主要是高频信号，在ZPW-2000A中分为四种：1700HZ、2000HZ、2300HZ、2600HZ。每一种载频又分为1型与2型。1型的频率是原频率加1.4HZ；2型的频率是原频率减1.3HZ。因此，载频共有8种。

8种载频频率是：

1701.4、1698.7、2001.4、1998.7、2301.4、2298.7、2601.4、2598.7

注：区间线路载频锁频逻辑是指机车信号设备进入区间时，与区间载频自动匹配的方法当收到1700或2300Hz的2系载频+25.7HZ时，机车信号自动切换为只接收1700/2300HZ载频的低频信息状态。

当收到2000或2600Hz的2系载频+25.7HZ时，机车信号自动切换为只接收2000/2600Hz载频的低频信息状态。

无论是侧线股道发车还是正线股道弯出发车，均采用列车压入发车进路最末一个区段时，该区段发送载频为+2的25.7Hz转频码，列车出清该区段后，恢复发送27.9Hz的低频检测信息。

fc=2001.4hz；         fd=20.2hz；

fc= 2598.7hz；           fd=25.7hz；

fc=2298.7hz；           fd=29hz；

三.欠采样原理

1.欠采样原理

欠采样是在测试设备带宽能力不足的情况下，采取的一种手段，相当于增大了测试设备的带宽，从而达到可以采样更高频率信号的能力。

ZPW2000型移频信号最高频率可达2611Hz,由Nyquist采样定理可知，采样频率应不小于5222Hz但处理器进行FFT运算时的采样点数太多会使计算时间加长，若采样点数过少会使频率分辨率(△f=fₛ/N)下降，进而导致检测精度降低。为了既能提高检测精度，又能降低采样点数，故对ZPW2000型移频信号采用欠采样技术。

欠采样技术就是指以低于奈奎斯特采样频率 K倍的采样频率进行无失真的采样过程。从频域上分析,信号的采样过程其实就是原信号频谱沿频率轴的搬移过程,要使信号不失真,则要求采样信号频谱在整个频域内不重叠。所以欠采样频率 fs 应满足如下关系：

｛fL - k ( fs /2) > 0且fH - k ( fs /2) ≤( fs /2)

即｛kfs < 2fL且( k + 1) fs ≥2fH

即 2fH/（k + 1）≤fs<2fL/k，其中k = 1, 2, 3, …, K; K =【fL /B 】向下取整。

这样就可保证原信号频谱中的上下截止频率都不与上边带重叠,即采样信号的频谱无混叠现象出现。

欠采样前后中心频率的关系式 ：

f0 =[k + 1/2 + ( - 1)^( k -1)×1/2] × fs/2 + ( - 1)^k ×f′0

其中 , fs 为采样频率; f′0 为欠采样后的中心频率; k 为选择 fs 时所对应的 k值。

根据欠采样无混叠条件计算：

2.运行效果

 

N=40960，满足奈奎斯特定理fs≥2fc

N=2048,，满足欠采样定理

对移频信号进行 FFT 快速傅里叶变换、欠采样技术，计算所需时间。对两个时间进行对比，证明采用欠采样技术可以有效减少移频信号的处理时间。

（运行代码见附页）

四.Zoomfft原理

1.引入Zoom_FFT算法的原因

FFT虽然迅速方便地实现了数字频谱分析的功能，但是对于希望在某些感兴趣的频率范围内以较高的分辨率进行频域分析，传统的 FFT 就做不到了。因为普通 FFT在频域的分辨率是有一定限度的，其分析频率范围总是从0Hz 开始到 1/2倍的采样频率为最高分析频率。设采样频率为fs，FFT点数为NFFT，频率分辨率△f的关系式为△f =fs/NFFT。可以看出，一旦采样频率fs和NFFT确定后，频率分辨率是一个固定值。

由于频率分辨率的限制，那么相邻两个离散点之间的频率无法获得，会导致部分有用的频率被漏掉，这就是所谓的栅栏效应。

解决这种问题的一种方法是可以增加数据的点数或者提高傅里叶变换的点数，但是这样会增加系统的硬件成本。不过可以通过某些方法在不改变硬件的基础上提高频率的分辨率，下面简单介绍一种经典的方法Zoom-FFT算法。

2.原理介绍

在实际的应用中，只需要对感兴趣的某段频率的信号进行分析，下面介绍一种基于复调制的频谱细化方法。

基于复调制的频谱细化方法

频移细化FFT的方法是基于离散傅立叶变幻的频移原理。设fk为所需要细化的频带的中心频率，对信号乘以exp(-j2pifkt)进行数字化平移后，原fk处的谱线移至频率轴的0处，用低通滤波器滤除所需细化的频段外的频率成分，并以细化倍数为间隔进行重采样，最后对重采样后的数据做FFT变换，并对变换结果重新排序。

3.细化FFT的具体步骤

1）欠采样

由奈奎斯特采样定理得，采样频率是信号最高频率的两倍以上，信号才能无失真恢复。这样，采样频率高于信号最高频率的两倍，一般称为过采样；采样频率低于信号最高频率的2倍，一般称为欠采样。

2）复调制移频

将采样信号进行频移（复调制），即乘以单位旋转因子exp(-j2pifkt)，这样就把频率原点由0处移到所需要细化的频率fk处，频率分量fk停留在频率为零的位置上，形成了一个以fk为频率零点的新的信号xk(t)。

3）低通滤波

遵循采样定理的原则，为防止采样信号的频率混淆，首先需要通过模拟低通抗混叠滤波器滤波或这顶足够高的采样频率fs。然后需要采集足够长度的信号数据，数据的长度为细化倍数D与NFFT的乘积，即为DNFFT。

4）重新采样

信号经过移频、低通滤波后，分析信号点数变少，再以较低的采样频率进行重新采样，在通过补零保证相同的采样点数时，样本的总长度加大，频谱的分辨率也就得到了提高。

设原采样频率为 ,采样点数为N,则频率分辨率为fs/N,现重采样频率为fs/D，当采样点数仍是N时，其分辨率为fs/(D*N)，分辨率提高了D倍。这样就在原采样频率不变的情况下得到了更高的频率分辨率。

5）复数FFT

重采样后信号实部和虚部是分开的，需要对信号进行N点复FFT，从而得出N条谱线，此时分辨率为，可见分辨率提高了D倍。

6）频率调整

经过算法运行后的谱线不为实际频率的谱线，需要将其反向搬移，转换成实际频率，进而得出细化后的频率。

fc=1701.4;  fd=11.4;  fa=11;  fs=2048时;

算得：((346.6-335.2)+(358-346.6))/2=11.4

fc=2301.4;  fd=11.4;  fa=11;  fs=2048时;

算得：((253.4-242)+(264.8-253.4))/2=11.4

fc=2298.7;  fd=13.6;  fa=11;  fs=2048时;

算得：((250.7-237.1)+(264.3-250.7))/2=13.6

（运行代码见附页）

通过以上三组数据可看出，经过zoom_fft细化后的频谱的频偏，即（（最高峰-左边的次高峰）+（右边的次高峰-最高峰））/2的结果，与所设置的低频调制信号的频率相等。由此可得，此代码正确，进行zoom_fft后确实能提高频率的分辨率，使有用的频率不被漏掉。

五、总结

在本次课程设计中，我们利用Matlab仿真工具对ZPW2000信息轨道电路移频信号检测方法进行了仿真分析。首先，我们用Matlab产生了轨道移频信号，这种移频信号的产生是采用频率调制的方法，即通过低频调制信号控制载频信号，形成振幅不变而频率随低频信号做周期性变化的信号。再利用欠采样技术，对其信号选取统一的采样频率进行采样，从而可以大大减少轨道移频信号检测的复杂性；通过分析信号的频谱特点，给出了基于FFT频谱校正的算法，即Zoom-FFT算法，来检测移频信号参数。分析表明，将Zoom-FFT频谱细化技术应用于通过轨道电路移频信号的检测是完全可行的，各个参数在误差指标的条件下都能准确、高精度地检测出来，可以完全满足较高的频率和分辨率的要求。

根据此次的研究，可以将该算法与硬件相结合，设计出一款移频信号参数测试仪表，这对于轨道移频信号检测方法的研究以及检测系统的开发也具有一定的指导意义。

参 考 文 献:

[1] 袁薇，孙景峰，樊文侠. 欠采样与ZFFT在移频信号检测中的应用.

[2] 王安，焦美鹏，罗世刚. 基于FFT频谱校正的移频信号检测算法的研究.

[3] 何林，彭莉峻. 移频信号解调算法研究与DSP实现.

[4] 樊文侠，孙景峰. 基于Matlab轨道移频信号检测的仿真分析.

六.代码部分

1.信号产生部分

clear,close all;
clc
tic%计算时间：开始
%8种载频频率
fc1=1701.4;
fc2=1698.7;
fc3=2001.4;
fc4=1998.7;
fc5=2301.4;
fc6=2298.7;
fc7=2601.4;
fc8=2598.7;
%18种低频调制频率
f01=10.3;
f02=11.4;
f03=12.5;
f04=13.6;
f05=14.7;
f06=15.8;
f07=16.9;
f08=18;
f09=19.1;
f010=20.2;
f011=21.3;
f012=22.4;
f013=23.5;
f014=24.6;
f015=25.7;
f016=26.8;
f017=27.9;
f018=29;
%ZPW-2000轨道移频信号设计
dt=0.0001;
t=0:dt:1;
fc=fc1;%载波频率========================可人工修改行
fd=f06;%调制频率========================可人工修改行
fa=11;%频偏
%产生一个以1/f01为周期的调制方波g
g=square(2*pi*fd*t,50);%周期为1/f0,占空比为50%的方波调制信号
figure;
subplot(3,1,1);
plot(t,g);
axis([0 0.1 -1.2 1.2]);
title(strcat('低频调制信号，频率为fd=',int2str(fd),'HZ'));
xlabel('时间');
%为了消除频率上的跳变点，对方波进行积分，使方波变成周期三角波信号，得到gs
gs=cumsum(g)*dt;%取积分
rfsk=cos(2*pi*fc*t+2*pi*fa*gs);%构建频域时移表达式
subplot(3,1,2);
plot(t,rfsk);
axis([0,0.1,-1.2,1.2]);
title('ZPW-2000轨道移频信号');
xlabel('时间');
ylabel('幅度');
%频域
df=0.01;
f=fc-3*fd:df:fc+3*fd;%频域显示范围
w=2*pi*f;
S=rfsk*exp(-1j*t'*w)*dt;
S1=abs(S);
subplot(3,1,3);
plot(f,S1);
title('ZPW-2000轨道移频信号的频域波形');
xlabel('f(Hz)');
ylabel('幅度');
fs=10000;
N=40960;%N点傅里叶变换
Y1K=fft(rfsk,N); %对移频信号进行FFT快速傅里叶变换
mag=abs(Y1K); %对Y1K进行取绝对值
df1 = fs/N; %频率分辨率
Y2K=mag(1:N/2)*2/N; %频谱图是对称图形，取N/2，乘2/N能得真正的幅度值 
f=(1:N/2)*df1; %每个序列号对应的频率
figure;
title('移频信号频谱');
xlabel('频率(Hz)');
stem(f,Y2K,'.');
toc%结束
disp(['直接运行FFT时间:',num2str(toc)]);

2.欠采样部分（时间比较）

dt=0.0001;
t=0:dt:1;
fc=1701.4;%载波频率
fd=11.4;%调制频率
fa=11;%频偏
%产生一个以1/f01为周期的调制方波g
g=square(2*pi*fd*t,50);%周期为1/f0,占空比为50%的方波调制信号
%为了消除频率上的跳变点，对方波进行积分，使方波变成周期三角波信号，得到gs
gs=cumsum(g)*dt;%取积分
rfsk=cos(2*pi*fc*t+2*pi*fa*gs);%构建频域时移表达式
fs=10000;
N=2048;%N点傅里叶变换
Y1K=fft(rfsk,N); %对移频信号进行FFT快速傅里叶变换
mag=abs(Y1K); %对Y1K进行取绝对值
df1 = fs/N; %频率分辨率
Y2K=mag(1:N/2)*2/N; %频谱图是对称图形，取N/2，乘2/N能得真正的幅度值 
f=(1:N/2)*df1; %每个序列号对应的频率
figure;
title('移频信号频谱');
xlabel('频率(Hz)');
stem(f,Y2K,'.');
toc%结束
disp(['欠采样所需时间:',num2str(toc)]);

3.Zoomfft部分

方案一

fe=325; %中心频率
tic
k=1:M; 
w=0.5+0.5*cos(pi*k/M); %Hanning函数，其中`k`是窗口内的采样点序号，`w`是Hanning窗的加权系数。
fl=max(fe-fs/(4*D),-fs/2.2); %频率下限
fh=min(fe+fs/(4*D),fs/2.2); %频率上限
yf=D*fl; %移频量
df=fs/D/N; %分辨率
f=fl:df:fl+(N/2-1)*df;%频率向量
xz=zeros(1,N/2);%幅度谱
wl=2*pi*fl/fs; %归一化角频率
wh=2*pi*fh/fs;
hr(1)=(wl-wh)/pi;
hr(2:M+1)=(sin(wl*k)-sin(wh*k))./(pi*k).*w;
hi(1)=0;
hi(2:M+1)=(cos(wl*k)-cos(wh*k))./(pi*k).*w;
%{这一段代码计算Zoom-FFT的窗口函数，
其中`wl`和`wh`是Zoom-FFT的下限频率和上限频率对应的角频率，
`hr`和`hi`是Zoom-FFT的实部和虚部窗口函数。}%
p=0:N-1;
w=0.5-0.5*cos(2*pi*p/N);%定义FFT的窗口函数，其中`w`是FFT窗口函数的加权系数。
xrz=zeros(1,N/2);
xiz=zeros(1,N/2);
L=10; %循环次数
for i=1:L% 进行Zoom-FFT的循环计算
for k=1:N% 对每个时域序列进行计算
kk=(k-1)*D+M;% 计算当前时域序列所在的位置
xrz(k)=x(kk+1)*hr(1)+sum(hr(2:M+1).*(x(kk+2:kk+M+1)+x(kk:-1:kk-M+1)));% 低通滤波
xiz(k)=x(kk+1)*hi(1)+sum(hi(2:M+1).*(x(kk+2:kk+M+1)-x(kk:-1:kk-M+1)));% 高通滤波
end
xzt=(xrz+1j*xiz).*exp(-1j*2*pi*(0:N-1)*D*fl/fs);
%复调制移频,使得(xrz+1j*xiz)的fl搬到xzt的零频点,乘以D将其细化D倍
xzt=xzt.*w;%为了保证重新采样后的信号在频谱分析时不发生频谱混叠，需进行抗混叠滤波，滤出需要分析的频段信号
xzt=xzt-sum(xzt)/N;
xzt=fft(xzt);%复数FFT
xz=xz+(abs(xzt(1:N/2))/N*2).^2;%幅度谱计算
end
xz=(xz./L).^0.5;%求平均
%{这一段代码是Zoom-FFT的主要计算过程，
其中包括了对信号进行移频、应用Zoom-FFT窗口函数、FFT计算和幅度谱的计算。
循环次数为Zoom-FFT的计算次数`L`。
最后，将幅度谱取平均并开方即可得到Zoom-FFT的结果。}%
toc
subplot(2,1,2);
plot(f,xz);%重新采样
axis([280 390 0 4])
grid on
xlabel('f(Hz)');ylabel('Amplitude')
title('Zoom-FFT细化后的频谱')
toc
disp(['移频信号进行 FFT 快速傅里叶变换、欠采样技术、Zoom-FFT 技术所需要的时间',num2str(toc)]);

方案二

fe=325; %中心频率
fs=2048; %采样频率
N=2048; %傅里叶变换点数
D=10; %细化倍数
M=200; %滤波器阶数
fc=1701.4;
fd=15.8;
fa=11;
dt=1/fs;
t=(0:N*D+2*M)/fs;
g=square(2*pi*fd*t,50);
gs=cumsum(g)*dt;
x=3*cos(2*pi*fc*t+2*pi*fa*gs);
%%复调制频率细化分析,调用格式：[y,freq]=exzfft_ma(x,fe,fs,N,D);
[y,freq]=exzfft_ma(x,fe,fs,N,D);
%{说明：其中输入变量x是被测信号，x的长度要大于等于nfft*D，fe是细化区间的中心频率，fs是采样频率，nfft是细化FFT的长度，D是细化倍数。输出变量y是细化FFT后的输出，是一个复数序列；freq是细化FFT后的频率刻度。}%
subplot(2,1,2);
plot(freq,abs(y),'k');
axis([280 390 0 4])%设定横、纵坐标显示范围
grid on
xlabel('f(Hz)');ylabel('Amplitude')
title('Zoom-FFT细化后的频谱')
%函数程序如下
function [y,freq,c]=exzfft_ma(x,fe,fs,nfft,D)
nt=length(x);           % 计算读入数据长度
fi=fe-fs/D/2;           % 计算细化截止频率下限
fa=fi+fs/D;             % 计算细化截止频率上限
na=round(0.5 * nt/D+1); % 确定低通滤波器截止频率对应的谱线条数
% 频移
n=0: nt-1;              % 序列索引号
b=n*pi* (fi+fa)/fs;     % 设置单位旋转因子
y=x.*exp(-1i*b);        % 进行频移 
b= fft(y, nt);          % FFT
% 低通滤波和下采样
a(1: na) =b(1: na);     % 取正频率部分的低频成分
a(nt-na+2 : nt) =b(nt-na+2 : nt); % 取负频率部分的低频成分
b=ifft(a, nt);          % IFFT 
c=b(1 : D: nt);         % 下采样
% 求细化频谱
y=fft(c, nfft) * 2/nfft;% 再一次FFT
y=fftshift(y);          % 重新排列
freq=fi+(0:nfft-1)*fs/D/nfft; % 频率设置
end