载波同步算法程序

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

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

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

∑∑-=-=-=+=

1

2

10

2*

|

)(|)()(N l N l j l y e N l y l y J πρ

（39）

因此我们可以估计出??

????

=||arg 21*J J πρ，这里的)(l y 指的是接收信号。 然而我们知道实际ρ的值会比1大（如前面提到的100ppm 的晶振偏移对应

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

∑∑-=--==+=

14/0

2

4

/21

4/0

*

|

)(|)4/()(N l j N l l y e

N l y l y K πρ （40）

故可以得到?

?

????=||arg 214

*K K π

ρ，可见短前导字的估计范围扩大到长前导字估计范围的4倍，也就是说精频偏估计的精度为粗频偏估计的4倍。结合上面提到的粗估

计和精确估计，可以得到

??

?

?

????????+??????=||24||arg 21

**K K J J π

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

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

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

()()u F I F I T 1f'n f f f ?+=?+?=? (42)

()F I u 'f n 1/f 'f ?+=?=? (43)

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

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

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

（1） 整数频偏估计

设()n ,m r s ）

，（S S N ...2,1n M ...1,2m ==和()n ,m r l ），（l l N ...2,1n M ...1,2m ==分别表

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

()()()n ,m n n ,m s e n ,m r l 'f i2l +=?π ）

，（l

l

N ...2,1n M ...1,2m == (44) 假设小数部分频偏已经得到补偿，(5-12)式可以改写成：

()()()n ,m n n ,m s e n ,m r l n i2l I +=π ）

，（l

l

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

()()()()n ,m N n n ,m X n ,m Z I N I l l +-=

），（l l N ...2,1n M ...1,2m == (46) 其中()n ,m Z l 是()n ,m r l 经过FFT 变换解调子载波后在频域上的表示， ()n ,m X l 是发射端preamble 中长训练符号()n ,m s l 在FFT 调制子载波前在频域的表达式，

()()I

N I l n n ,m X -表示()n ,m X l 长度为l

N 的圆周移位，()n ,m N 是噪声项在频域的表达

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

()()i n ,m X n ,m Z max n ?l N 1n i i

I l

+=*=∑ 12/N ...1,0,1..2/N i l l ---= (47)

文献指出此算法估计范围为(-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 信号可以表示成：

k N

k j k k

f e

r z ωπδ+=/2 （48）

这里f

f ?=

ε

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

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

22n s σσ和。在接收端，相对应于循环前缀的采样值被去除，余下的N 个采样值

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

,110

2∑-=-=N k N

kn

j k n e z N Z π 10-≤≤N n （49） 前人已研究表明循环前缀和导频符号可以用来跟踪频偏。接下来先回顾一下这两种算法并对其作一定的修改。为了能较好的估计频偏的影响，我们这里假设帧和

符号定时是完全正确的，且信道被假设为慢衰落，在下面的讨论中，信噪比被定义为2

2

/n s SNR σσ=。 A ．基于导频的频偏跟踪法

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

()()

?

??????=??+∑∈++P n n D m n m n D m n m f g

C C Z Z

D N

N N ,,*

,*,arg 21πδ （50）

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

{}{}n

D m n

m C

C ,,+和指的是在第

D m m +和个OFDM 符号的第n 个子载波上传输的导频

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

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

10,110/)(2-≤≤+???

???=∑-=+N k e H X N z k N n N n k j n n k f ωδπ （51）

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

n n N

N j f f n n n W I e

N N H X Z f ++??

=-/)1()

/sin()sin()(πδπδπδ （52）

可见，这里的n Z 包括三个部分，第一部分是经过信道转换函数修改过的n X ，从上面的式子可见这部分由于受频偏的影响经历了幅度的衰减和相位的偏移。第二部分是由频偏引起的ICI ，具体表达式为： ()N n l j N

N j f

f N n

l l l

l

n e e

N n l N H X I f /)(/)1(1

,0)

/)(sin()

sin(----≠=??+-=

∑ππδδ

ππδ （53）

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

0}{=n I E 以及.5.0||),(sin ||||5947.0||2222≤≤f f n H X I E δπδ可见，n I 和频偏成一定的比例，随着频偏的增大，PTA 算法的性能逐渐下降。 B ．基于循环前缀的频偏跟踪法

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

N

k k p f πδ2)(= （54）

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

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

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

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

))()(arg(

)(0

*

∑-=--?-=L

N m g N m k z m k z k p （55）

其中k 是最近输入采样值的索引，针对IEEE802.11a 协议，这里取16=g N ，

8=L ，64=N 。若定时正确，则

∑-=--?-L

N m g N m k z m k z 0

*

)

()(的相关峰值可以达到，则其

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

π

δ2)(max

k p f = （56） 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-1 60-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)*8 0+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+sq rt(-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