实验名称实验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=Δff1=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)nsinπ(v+n)2π(v+n)若要求移频信号归一化功率为1，则移频信号的各个谐波分量的归一化值为A̅n=c√∑|c n|2+∞n=−∞=c√|c0|2+2∑|c n|2+∞n=1内容2ZPW-2000制式，仿真结果如下主谐波频率f=1T=1(0.4937−0.009338)/13≈26.8394Hz偏移4Hz。
归一化后，在4Hz±0.5f上进行积分后开方，得有效值0.6458。
国产18信息制式，仿真结果如下f=1T=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的值，记录下使nS=∑|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);elseF1Set=[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<type)a=sort(index_sorted(n-k+1:n));m=1;t(1)=a(1);for i=2:kif (a(i)~=a(i-1)+1)m=m+1;t(m)=a(i);elseif (X(t(m))<X(a(i)))t(m)=a(i);endendendk=k+1;endindex=(t-257)*2;%确定f1b=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:typer(j)=b(j)/(i*0.1);tmp=tmp+abs(r(j)-round(r(j)));endif (tmp<=0.02*delta)f1(step)=i*0.1;endendif (f1(step)<7)delta=delta+1;elsedelta=100;endend%将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);endendf1(step)=f1tmp;end%输出结果disp(f1);分别从2009010990_dat1和2009010990_dat3读取数据，运行程序进行分析，得到结果如下表所示。
与任务3算法原理相似，但需要加入对△f c和△f1的处理，因此程序结构上有所改变。
仍然从文件中读入采样数据，用fft()函数进行频谱分析，并抽取频谱中的n个最大值，将其对应频率放入数组f。
先以0.01为步长，不断列举△f c的值，再以0.01为步长，不断列举△f1的值，并枚举f1为标准值，记录下使nS=∑|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);elseF1Set=[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];endfor step=1:10load 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<type)a=sort(index_sorted(n-k+1:n));m=1;t(1)=a(1);for i=2:kif (a(i)~=a(i-1)+1)m=m+1;t(m)=a(i);elseif (X(t(m))<X(a(i)))t(m)=a(i);endendendk=k+1;endindex=(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:5index0=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);endendendif (thedelta<=deltamin)deltamin=thedelta;f1index=theindex;dfcfinal=dfc;df1final=thedf1;endenddisp([dfcfinal,df1final,F1Set(f1index)]);end分别从2009010990_dat2和2009010990_dat4读取数据，运行程序进行分析，得到结果电94 胡天骐2009010990其f1应为30左右，远超出了标准值的取值范围，因此程序识别错误，认为f1=15Hz。
我查看了其他同学的一些数据，没有发现类似这样的数据。
事实上在任务3中，制式为国产18信息的第9条数据，同样出现了f1=30Hz，只不过任务3中采用了另一种算法，程序没有出错。
误差分析由于频谱分析后频率最小间隔为2Hz，因此对于单个峰值，误差在1Hz之内。
而谐波分析过程中，采用了多个峰值所对应的频率，造成误差应当在0.1Hz数量级。
需要注意的是，这一误差可以通过增加采样频率来减少。
忽略频谱离散造成的误差，该算法中△f c，△f1，f1的误差都在0.01Hz以内，并可以通过缩减步长而减小。
11。