载波同步算法程序

5.2.1 载波频偏的捕获 A．前导序列估计

由于发送端和接收端载波频率的不同，每一个采样信号在时间t时包含一个未知的相位因素ej2??fct，这里?fc指的是未知的载波频偏。这个未知的相位因素在接收端必须被估计和补偿，否则子载波的正交性将会被破坏。例如，当载波频率为5GHz时，那么100ppm的晶振偏移相对应的频率偏移为50kHz。若符号的间隔时期为T?3.2us，那么?fcT?1.6。

前面对IEEE802.11a的分析，我们知道在它的帧结构中包含10个完全相同的短前导序列和两个相同的背靠背的长前导序列。其中短前导序列主要用于自动增益控制、分集选择、定时估计以及粗频率估计，而长前导序列主要用于信道估计和精确的频率估计。故结合这两个序列可以较精确的估计出载波的频偏，其中具体算法主要是利用它们良好的相关性[21]。首先设???fcT，则两个长前导字的相关值为：

J??y(l)y(l?N)?e?|y(l)|*?j2??l?0l?0N?1N?12 （39）

J*?，这里的y(l)指的是接收信号。 因此我们可以估计出??1arg???2??|J|? 然而我们知道实际?的值会比1大（如前面提到的100ppm的晶振偏移对应

的?为1.6），而长前导对其估计只能限制在?0.5内，故必须使用短前导字对其进行粗频率估计。短前导字的相关值为：

N/4?1 K??y(l)y(l?N/4)?e*l?0?j2??/4N/4?1l?0?|y(l)|2 （40）

K*?，可见短前导字的估计范围扩大到长前导字估计范围故可以得到??1arg???42??|K|?的4倍，也就是说精频偏估计的精度为粗频偏估计的4倍。结合上面提到的粗估

计和精确估计，可以得到

?J*??41 ??arg????2?|J|???2??K*????? （41） |K|???这里??指的是向下取整的意思。通过仿真结果表明此算法在理想状态下可估计的最大范围为?2.0（归一化值）。 B．载波频偏综合估计算法[12]

对于复杂的无线信道来说，上述算法?2.0的频偏估计范围偏小了一些。因此在文献[12]中提出了载波频偏估计算法，它的主要思想上在捕获阶段利用报头中长短训练符号的特殊结构，在频域和时域上分别实现整数载波频偏和分数载波频偏的捕获。

载波频差Δf 可以分为整数载波频差ΔfI 部分和分数载波频差ΔfF 部分，见下式

?f??fI??fF??nI??f'F??1Tu? (42)

?f'??f/?1Tu??nI??f'F (43)

整数载波频差ΔfI 是子载波频宽1/Tu 的整数倍，nI、?f'F和?fF分别为ΔfI 、ΔfF和Δf相对于子载波频宽1/Tu的归一化值。频偏捕获阶段的任务是在尽可能短的时间内准确地建立子载波间的正交关系，即要求能够同时较快地估计出nI和

?f'F。本文提出的捕获算法同样是基于IEEE802.11aPreamble 结构分两步，首先

利用报头Preamble 中长训练符号已知和FFT 变换的频域移位特性在频域上捕获整数载波频偏nI ，然后利用长短训练符号的结构特性在时域内捕获分数载波频偏?f'F。

（1） 整数频偏估计

设rs?m,n? 和rl?m,n? 分别表（m?1,2...Ml，n?1,2...Nl）（m?1,2...MS，n?1,2...NS）示接收端FFT变换之前preamble中第m个短（长）符号的第n个采样点，假设经过高斯白噪声信道，由于存在频偏，对于接收到的2个长训练符号中的任何一个都可以写成：

rl?m,n??ei2??f'sl?m,n??n?m,n? (44) （m?1,2...M，n?1,2...N）ll假设小数部分频偏已经得到补偿，(5-12)式可以改写成：

rl?m,n??ei2?nIsl?m,n??n?m,n? (45) （m?1,2...M，n?1,2...N）ll式中n(m,n)为高斯白噪声项。对上式两边同时做FFT变换，由于FFT变换频域移位特性可知，时域上的整数频率偏差经过FFT变换以后成为频域上的圆周移位。

(46) （m?1,2...Ml，n?1,2...Nl）Zl?m,n???Xl?m,n?nI??NI?N?m,n?

其中Zl?m,n?是rl?m,n?经过FFT变换解调子载波后在频域上的表示， Xl?m,n?是发射端preamble中长训练符号sl?m,n?在FFT调制子载波前在频域的表达式，

?Xl?m,n?nI??N表示X?m,n?长度为N的圆周移位，N?m,n?是噪声项在频域的表达

Ill式。由于Xl?m,n?对接收端来说是已经信号，对任何一个经过FFT解调子载波的长训练符号，我们都用一个带偏移量i的本地长训练符号与之相关，并逐步改变i值，由(46)性质可知，最大相关时刻对应的i值就是所求的整数载波频偏。

n?I?max?Zi?m,n?X?l?m,n?i? i??Nl/2..?1,0,1...Nl/2?1 (47)

in?1Nl文献指出此算法估计范围为(-32,31)(归一化值)，下面将进行仿真来验证，并将两者进行比较。 C．两种算法的仿真比较

在给定相同条件下，即相同的一帧数据、信道和频偏，改变信道和频偏参数对两种算进行比较。

? 通常设信噪比（SNR）为25～30，这里给定为25：

5.2.2 载波频偏的跟踪

前面对频偏捕获之后，我们认为残余的频偏足够小，可以使用跟踪模式进行补偿，而为了使得跟踪算法的可靠性，捕获后残余的频偏精度必须在0.5个子载波的空间范围内。到目前为止，可使用的频偏跟踪算法总共可以分为以下三种[36]：基于导频的算法（pilot tone-aided PTA），基于循环前缀的算法（cyclic prefix-based CPB）以及直接判决法(decision-directed DD)。PTA算法估计频偏主要通过在特定的子载波位置插入导频以及把接收到的符号和已知的导频符号相关，而对于IEEE802.11a协议来说，插入导频这是没有必要的，因为它的帧结构中已经含有导频。CPB算法一般是基于一个OFDM符号内循环前缀的周期性的特性来实施的。而DD算法大部分和PTA算法是相同的，除了一点以外，DD算法是把接收到的符号和假设符号相关，而PTA是和已知的符号相关。然而一般不推荐使用DD算法，因为在时延和错误传播的环境下，它的操作性能会大大的下降。因此接下去主要讲解PTA和CPB两种算法。

先大概的讨论两种算法，然后把这两种算法放在相同的环境下比较它们的性能。这种比较有利于决定在特定的应用下选用哪一种算法。

首先我们假设当存在频偏时接收端OFDM信号可以表示成：

zk这里?f??rkej2??fk/N??k （48）

??f是信道的相关频偏（即频偏?和子载波间隔?f的比值），而?是白

色高斯噪声。而数据和噪声被假设为非相关，且他们的随机方差值分别为

?s2和?n2。在接收端，相对应于循环前缀的采样值被去除，余下的N个采样值

被用于解调。经过FFT之后，第n个子载波上的符号可表示成：

1N?1?j2?Nzke, 0?n?N?1 （49） ? Zn?Nk?0前人已研究表明循环前缀和导频符号可以用来跟踪频偏。接下来先回顾一下这两种算法并对其作一定的修改。为了能较好的估计频偏的影响，我们这里假设帧和

kn符号定时是完全正确的，且信道被假设为慢衰落，在下面的讨论中，信噪比被定义为SNR??s/?n。 A．基于导频的频偏跟踪法

假设一个OFDM符号含有N个子载波，其中Np个子载波被导频符号调制。让P表示Np个导频符号对应的索引值，在IEEE802.11a协议中，即为[12，26，40，54]。则基于导频的算法表达式如下[25][36]：

22

N?NgN1??**?D??f??arg??Zm,nZm?D,nCm,nCm?D,n? （50）

2??n?P?????这里m和m?D分别代表第m和m?D个OFDM符号，在这里我们取D为1。而

?C?和?Cm,nm?D,n?指的是在第m和m?D个OFDM符号的第n个子载波上传输的导频

符号。然而经研究发现其实?Cm,n?和?Cm?D,n?可以不用乘，因为我们只要取对应的幅角，而?Cm,n?和?Cm?D,n?是一个常数。

由于PTA算法包括解调过程，即它经过FFT的转换过程，这势必要受到ICI（载波间的干扰）。因此我们把（48）式重新写成如下式：

1?N?1j2?k(n??f)/N?zk?XHe??k,0?k?N?1 ?nn??N?n?0? （51）

这里Hn是第n个载波频率上对应的信道转换函数，相应的（49）式可改写为：

Zn?(XnHn)?sin(??f)Nsin(??f/N)?ej??f(N?1)/N?In?Wn （52）

可见，这里的Zn包括三个部分，第一部分是经过信道转换函数修改过的Xn，从上面的式子可见这部分由于受频偏的影响经历了幅度的衰减和相位的偏移。第二部分是由频偏引起的ICI，具体表达式为： In?l?0,l?n??XH?Nsin(?(l?n??llN?1sin(??f)f)/N)?ej??f(N?1)/N?e?j?(l?n)/N （53）

假设这里用到的数据具有零均值和互不相关的特性，Moose已证明了下面性质

E{In}?0以及E|In|2?0.5947|X|2|H|2sin2(??f),|?f|?0.5.可见，In和频偏成一定的比例，随着频偏的增大，PTA算法的性能逐渐下降。 B．基于循环前缀的频偏跟踪法

从公式（48）中可见，载波频率偏差?f的存在导致了接收采样值z(k)的相位偏差p(k)，

2??fkN p(k)? （54）

因此两个采样值z(k1)和z(k2)之间相位误差的差值其实是一个有关它们频率偏差和时延差的函数。如果在z(k1)和z(k2)之间最初的相位差值是已知的，而且所有其他的相位失真都不存在，则我们就可以从相位误差的差值中得到相对频偏?f。

我们知道循环前缀在时域的采样值其实就是OFDM符号最后Ng数据采样值的复制值，当OFDM符号通过一个长度为L的信道时，接收到的Ng?L?1采样值依然是原样，只是相位发生了2??f的偏差。因此频偏可以通过使用Ng?L?1个采样值中的每一个去估计，为了提高估计的精度，可以把Ng?L?1个估计值取平均。

当然在很多文献中是使用Ng个估计值的平均值去估计频偏的，然而我们注意到接收到的OFDM符号最前面的L?1个采样值通常受到前面一个OFDM符号的干扰，也就是我们说的符号间的干扰（ISI）。已经证实了使用Ng?L?1个采样值更能可靠的估计频偏。

下面给出根据循环前缀相关性的频偏跟踪算法：

Ng?L p(k)?arg(m?0?z(k?m)?z(k?m?N)) （55）

*L?8，其中k是最近输入采样值的索引，针对IEEE802.11a协议，这里取Ng?16，

Ng?LN?64。若定时正确，则

m?0?z(k?m)?z(k?m?N)*的相关峰值可以达到，则其

对应的相位可以记为p(k)max，它就等于保护间隔内的采样值和对应的OFDM数据采样值之间的平均相位偏移。由于相关的一对采样值被N个采样值隔开，从而可以得到精确的频偏估计：

p(k)max （56） 2? ?f?C．PTA和CPB两种算法的比较

在给定相同的条件下比较它们的性能。这种比较有利于决定在特定的应用下选用哪一种算法。下面设定剩余微小的频偏为0.02，用红线表示CPB算法的跟踪曲线，由用蓝线表示PTA算法的跟踪曲线图。

图19 两种算法跟踪结果的比较图

由仿真结果可以看出，CPB算法实现了快速准确地跟踪目的。然而PTA算法却没有实现预想的效果。从两者的图可以比较和分析两种算法。可以看出PTA的这种算法不适合用于IEEE802.11a协议，首先因为基于导频的算法多用于信道估计，其次该协议的规定中只是插入了四个的导频，如果增加插入的导频点数该算法也

是同样有效。相比较而言CPB算法更适用于该协议的载波同步。

附录（一）载波同步算法仿真程序

程序一：前导序列算法

function fre_offset=carry_acq(input,channel_fre_offset) input_1(1:160)=fft(input(1:160)); input_1(161:320)=fft(input(161:320)); figure(2)

title('前导序列算法比较图') subplot(311)

plot(abs(input_1))

xlabel('发送端的前导序列的频谱');

in=AWAG_channel_backup(input,25,channel_fre_offset); in_1(1:160)=fft(in(1:160));

in_1(161:320)=fft(in(161:320)); subplot(312) plot(abs(in_1))

xlabel('经过信道前导序列的频谱'); start_time=321; corr1=0; for i=1:16

corr1=corr1+in(start_time-160-i)*conj(in(start_time-160-i-16)); end

phase1=angle(corr1); corr2=0; for i=1:64

corr2=corr2+in(start_time-i)*conj(in(start_time-i-64)); end

phase2=angle(corr2);

phase=round(phase1*4/pi/2)*2*pi+phase2; fre_offset=phase/2.0/pi;

v=fre_offset-channel_fre_offset;

temp1=AWAG_channel_backup(input(1:320),1000,v); temp2(1:160)=fft(temp1(1:160)); temp2(161:320)=fft(temp1(161:320)); subplot(313)

plot(abs(temp2))

xlabel('经过频偏补偿后的前导序列频谱');

程序二：频偏综合估计算法

function fre_offset=carry_acq_improve(in,input,freoffset) in_1(1:160)=fft(in(1:160));

in_1(161:320)=fft(in(161:320)); subplot(311) plot(abs(in_1))

xlabel('发送端的前导序列的频谱');

z_stream=AWAG_channel_backup(in,25,freoffset); z_stream_1(1:160)=fft(z_stream(1:160)); z_stream_1(161:320)=fft(z_stream(161:320)); subplot(312)

plot(abs(z_stream_1))

xlabel('经过信道后前导序列的频谱'); Z_temp(1:64)=z_stream(193:256);

Z_stream=b_fft(Z_temp); % 将接收到的信号的长前导字序列变换到频域，以便在频域上捕获大频偏 X_stream=input;

Z_stream_double(1:64)=Z_stream(1:64); Z_stream_double(65:128)=Z_stream(1:64); Z_stream_double(129:192)=Z_stream(1:64); spot=0; max=0; i=0;

for i=-32:31

corr1(i+33)=sum(conj(X_stream(1:64)).*Z_stream_double(i+65:128+i)); % 通过计算本地前导序列和接收到信号的前导序列的相关性来找出最大的n值 spot_temp(i+33)=abs(corr1(i+33)); if spot_temp(i+33)>max max=spot_temp(i+33); n=i; end end

fre_offset_int=n; % 最大的n值就是频偏的整数部分 N=64; L=16;

length_temp=size(in); length=length_temp(2); for i=1:length

z_stream_c(i)= z_stream(i)*exp(-j*2*pi*fre_offset_int*i/N); %将捕获到的整数部分先补偿回去再进行小数部分的估计 end

start_time=321; corr2=0;

for i=1:16

corr2=corr2+z_stream_c(start_time-160-i)*conj(z_stream_c(start_time-160-i-16)); end

phase1=angle(corr2); corr3=0; for i=1:64

corr3=corr3+z_stream_c(start_time-i)*conj(z_stream_c(start_time-i-64)); end

phase2=angle(corr3);

phase=round(phase1*4/pi/2)*2*pi+phase2; fre_offset=phase/2.0/pi+fre_offset_int; v=freoffset-fre_offset;

temp1=AWAG_channel_backup(in(1:320),1000,v); temp2(1:160)=fft(temp1(1:160)); temp2(161:320)=fft(temp1(161:320)); subplot(313)

plot(abs(temp2))

xlabel('经过频偏补偿后的前导序列频谱');

程序（三）：基于循环前缀算法

function final1(frameout, offset) N=64; L=16;

length_temp=size(frameout); length=length_temp(2);

frameout_data=frameout(401:length); for i=1:length-400

in(i)=frameout_data(i)*exp(j*2*pi*offset*i/N); end

num_sym=length/80-5; ri=in;

T=4; % symble period (us)

T_fft=T*N/(L+N); % T_fft=3.2 (us) Kd=1; % awgn Ko=1;

K=Kd*Ko; % K=Kd*Ko , assuming any value , but it is determined by loop BL=0.033/T; kesai=0.707;

C1=(2*kesai/K)*(2*BL*T/(kesai+1/4/kesai));

C2=(2*BL*T/(kesai+1/4/kesai))^2/K; fai=0;

deltafai(1)=0; for n=1:num_sym if n==1

fo(1)=0; % (C1+C2)*deltafai(1) else

fo(n)=fo(n-1)+(C1+C2)*deltafai(n)-C1*deltafai(n-1); end

for m=1:L+N if n==1

rio(m)=ri(m)*exp(-j*fai); %out(n,m)=rio(m); else

rio(m)=ri((L+N)*(n-1)+m)*exp(-j*fai); %out(n,m)=rio(m);

fai=Ko*2*pi*T_fft/N*fo(n-1)+fai;

end end

est(n)=sum(conj(rio(4:L)).*rio(4+N:L+N)); if n==1

deltafai(n+1)=imag(est(n)); else

deltafai(n+1)=deltafai(n)+imag(est(n)); end end

hold on

n=1:num_sym-1-1;

plot(Ko*fo(n)*T_fft,'r') xlabel('OFDM符号数'); ylabel('频率误差');

程序（四）：基于导频算法 function pta(in,offset)

a1=AWAG_channel_backup(in,25,offset); length1=length(a1); a2=a1(193:length1);

length2=fix((length(a2)-128)/80); a3(1:64)=b_fft(a2(1:64)); a3(65:128)=b_fft(a2(65:128)); for i=1:length2

a3(128+(i-1)*64+1:128+(i-1)*64+64)=b_fft(a2((i-1)*80+1+128+16:(i-1)*80+80+128)); end

Bl=0.00008; kesai=0.707; k=1;

C1=(4*kesai*Bl/(kesai+0.25/kesai)); C2=((2*Bl/(kesai+0.25/kesai))^2); e1(1:length2+1)=0; e2(1:length2+1)=0; e3=0;

sum_angle(1:length2+1)=0;

for i=3:length2; e3=0;

for n=1:64

e3=mod(i*e2(i-1)+e3,2*pi);

a3((i)*64+n)=(a3((i)*64+n)*(exp(-j*((e3))))); end

R1(i)=(a3((i)*64+12))*(conj(a3((i-1)*64+12))); R2(i)=(a3((i)*64+26))*(conj(a3((i-1)*64+26))); R3(i)=(a3((i)*64+40))*(conj(a3((i-1)*64+40))); R4(i)=(a3((i)*64+54))*(conj(a3((i-1)*64+54)));

sum_angle(i)=(((imag(R1(i))+imag(R2(i))+imag(R3(i))+imag(R4(i)))))+ sum_angle(i-1);

e3(i)=(C2*sum_angle(i));%这个值很小可以忽略不计

e1(i)=(C1+C2)*sum_angle(i)-C1*sum_angle(i-1)+e1(i-1); e2(i)=(e1(i))+e2(i-1); end

plot(e2/(2*pi));

xlabel('OFDM符号数'); ylabel('频率误差');

程序（五）：加噪

function [gsrv1, gsrv2]=gngauss(m, sgma) if nargin==0, m=0; sgma=1;

elseif nargin==1, sgma=m;

m=0; end; u=rand;

z=sgma*(sqrt(2*log(1/(1-u)))); u=rand;

gsrv1=m+z*cos(2*pi*u); gsrv2=m+z*sin(2*pi*u);

程序（六）：信道 function

[channel_out]=AWAG_channel_backup(in,snr_in_dB,channel_fre_offset) echo off

%snr_in_dB = 60; a=size(in); len=a(2);

%--------------Channel------------ % Only BPSK Eb = 0.0125;

SNR = 10^(snr_in_dB/10); sgma=sqrt(Eb/(2*SNR));

for i=1:len

[ni,nq]=gngauss(sgma);

channel_out(i)=in(i)*exp(sqrt(-1)*2*pi*channel_fre_offset*i/64)+ni+sqrt(-1)*nq; end

程序（七）：64点FFT变换 function out=b_fft(in) for n=1:1:64 out(n)=0;

for k=-32:1:31

out(n)=out(n)+in(k+33)*exp(-1*sqrt(-1)*2*pi*k*(n-33)/64); end end

程序（八）：64点IFFT变换 function out=b_ifft(in) for n=1:1:64 out(n)=0;

for k=-32:1:31

out(n)=out(n)+in(k+33)*exp(sqrt(-1)*2*pi*k*(n-33)/64); end

out(n)=out(n)/64; end

本文来源：https://www.bwwdw.com/article/gla8.html