各位高手你们好，在下有理了。我最近在写有关ofdm的同步的问题，主要是SC算法，

我已经把M(d)的图像仿真出来了，也找到了M(d)取最大时的那个点，还有小数频偏也算算出来了，但是我不会算'定时估计均方差'和‘频率估计均方差’，请给位高手指点一下 在下感激不尽!!!

我写的代码 ，请高手们看一下那儿出错了：

前半部分仿真可能是对的 ，但是后面求频偏估计均方误差与定时估计均方误差的图像时，就不对了

希望告诉们能给我解答一下

clear

clc

gilen=32;

para=512;   % Number of parallel channel to transmit (points)

fftlen=512; % FFT length

noc=512;    % Number of carrier

nd=6;       % Number of information OFDM symbol for one loop

ml=2;       % Modulation level : QPSK

sr=250000;  % Symbol rate

br=sr.*ml;  % Bit rate per carrier

%gilen=32;   %Length of guard interval (points)

fullen=gilen+para;  %ofdm符号总长

deltad=250;  %时延

deltaf=0.25;  %频偏

deltaderr=0;%时延均方误差

deltaferr=0;%频偏均方误差

derr=zeros(1,11);%时延均方误差

ferr=zeros(1,11);%频偏均方误差

%f=-1;

error=zeros(1,11);

o=0

%************************** main loop part **************************

for ebn0=0:2:20

o=o+1;

%************************** main loop part **************************

nloop=100;  % Number of simulation loops

noe = 0;    % Number of error data

nod = 0;    % Number of transmitted data

eop=0;      % Number of error packet

nop=0;      % Number of transmitted packet

for iii=1:nloop

%************************** transmitter*********************************

%**************************训练序列 ****************************

pn = rand(1,para)>0.5;

pn = reshape(pn,para/2,2);

[ipn0,qpn0] = qpskmod(pn,para/2,1,ml);

kmod=1/sqrt(2); %  sqrt : built in function

ipn0=ipn0.*kmod;

qpn0=qpn0.*kmod;

sym = ipn0 + i*qpn0;

symbuf = zeros(para,1);

symbuf(1:2:para,1) = sym;

%train1 = symbuf*7;                                %为什么乘以7

train1 = symbuf*15;

clear sym;

clear symbuf;

pn1 = rand(1,para)>0.5;

pn1 = reshape(pn1,para/2,2);

[ipn1,qpn1] = qpskmod(pn1,para/2,1,ml);

kmod=1/sqrt(2); %  sqrt : built in function

ipn1=ipn1.*kmod;

qpn1=qpn1.*kmod;

sym = ipn1 + i*qpn1;

symbuf = zeros(para,1);

symbuf(1:2:para,1) = sym;

clear sym;

pn2 = rand(1,para)>0.5;

pn2 = reshape(pn2,para/2,2);

[ipn2,qpn2] = qpskmod(pn2,para/2,1,ml);

kmod=1/sqrt(2); %  sqrt : built in function

ipn2=ipn2.*kmod;

qpn2=qpn2.*kmod;

sym = ipn2 + i*qpn2;

symbuf(2:2:para,1) = sym;

train2 = symbuf*5;                                 %为什么乘以5

%train2 = symbuf;

clear sym;

clear symbuf;

v = zeros(1,para);

kpois=1:2:para;

v(kpois) = sqrt(2)*train2(kpois)./train1(kpois); % 含有 训练序列对应子载波的相位差信息

%************************** Data generation

%****************************

seldata=rand(1,para*nd*ml)>0.5;  %  rand : built in function

%****************** Serial to parallel conversion

%***********************

paradata=reshape(seldata,para,nd*ml); %  reshape : built in function

%************************** QPSK modulation

[ich,qch]=qpskmod(paradata,para,nd,ml);

kmod=1/sqrt(2); %  sqrt : built in function

ich1=ich.*kmod;

qch1=qch.*kmod;

%******************* IFFT ************************

x=ich1+qch1.*i;

x(:,1)=train1;

x(:,2)=train2;

y=ifft(x); %  ifft : built in function

ich2=real(y);   %  real : built in function

qch2=imag(y);   %  imag : built in function

%********* Gurad interval insertion **********

[ich3,qch3]= giins(ich2,qch2,fftlen,gilen,nd);

fftlen2=fftlen+gilen;

%********* Attenuation Calculation *********

spow=sum(ich3.^2+qch3.^2)/nd./para;  %  sum : built in function

attn=0.5*spow*sr/br*10.^(-ebn0/10);

attn=sqrt(attn);

%***************************  Receiver  *****************************

%*************************加频率偏移和时延***************************

[ich3,qch3]=delay(ich3,qch3,length(ich3),deltad);

y=ich3+qch3.*i;

for k=1:length(ich3);

y(k)=y(k)*exp(i*2*pi*deltaf*k/para);

end

ich3=real(y);

qch3=imag(y);

%***************** AWGN addition *********

[ich4,qch4]=comb(ich3,qch3,attn);

y=ich4+qch4.*i;

%*****************小数频率偏移和时延估计*********************************

for d=1:para

for n=1:para/2

z1(n)=conj(y(d+n+para/2))*y(d+n);

end

p1(d)=sum(z1);

for n=1:para/2

z2(n)=abs(y(d+n+para/2)).^2;

end

p2(d)=sum(z2);

end

p=abs(p1).^2./(p2.^2);

plot(p);

%c=linspace(-150,150,256);

%stem(c,p);

%grid on;

[a a1]=max(p); % 求最大值

th=a*0.9;

xulie1=p(1:a1-1);

xulie2=p(a1+1:end);

[j1 j2]=min(abs(xulie1-th));

[j3 j4]=min(abs(xulie2-th));

if(mod(j2+j4+a1,2)==0)

j5=(j2+j4+a1)/2;

else

j51=(j2+j4+a1-1)/2;

j52=(j2+j4+a1+1)/2;

[j53 j54]=max([p(j51) p(j52)]);

j5=j54+j51-1;

end

%********************求小数陪频偏***********************

f=-angle(p1(j5))/pi;

%这是一种 最大似然方法联合实现符号定时同步和载波同步的方法

% 具体的推导 见 《OFDM移动通信技术原理与应用》 4.6节 P104

%m =m-1;%  公式推导中的d实际上是delay 了 d-1 个样值   wm added

%********************定时估计方差***********************

%********************小数倍频偏补偿***********************

%r1=y;

%ich5=real(r1);

%qch5=imag(r2);

for k=1:length(ich4);

y(k)=y(k)*exp(i*2*pi*(-f)*k/para);

end

%********************整数倍频偏估计***********************

r1 = fft(y(gilen+deltad+1:para+gilen+deltad));

r2 = fft(y(2*gilen+deltad+para+1:2*para+2*gilen+deltad));

for n=2:2:para/2-1

for k=1:2:para/2-n

z11(k)=conj(conj(r1(k+n))*r2(k+n))*v(k);

end

r11(n)=sum(z11);

for k=para/2-n+2:2:para/2

z22(k)=conj(conj(r1(k+n-para/2))*r2(k+n-para/2))*v(k);

end

r22(n)=sum(z22);

x=r11+r22;

end

for k=1:2:para/2

z33(k)=2*abs(r2(k)).^2;

end

q=sum(z33);

F=abs(x).^2./(q.^2);

b=max(F);   % 求最大值

g=find(F==b) ;   % 求最大值所在位置

%plot(abs(x));

%plot(F);

%********************频偏的估计值***********************

l=g+f; %频偏估计值

%********************同步估计均方误差计算***********************

deltaderr=deltaderr+(j5-16-deltad)^2;

deltaferr=deltaferr+(f-deltaf)^2;

end

derr(o)=sqrt(deltaderr/(nloop))

ferr(o)=sqrt(deltaferr/(nloop))

end

figure(1);

grid on

semilogy(0:2:20,ferr,'x--b');

hold on

semilogy(0:2:20,ferr,'x--b');

hold on