park.m

另外一个OFDM频偏估计的代码

clear all;
close all;
Numloop=500;
N=1024;
Tu=102.4e-6; %useful OFDM symbol period
T=Tu/N; %baseband elementary period
G=1/32; 
M=G*N;
delta=G*Tu; %guard band duration
Ts=delta+Tu; %total OFDM symbol period
FS=(1+G)*N; %IFFT/FFT length
q=10; %carrier period to elementary period ratio
fc=q*1/T; %carrier frequency
Rs=4*fc; %simulation period5
t=0:1/Rs:(Tu+Ts-1/Rs);
tt=0:T:Tu;
%Data generator
for SNR=0:5:30
    SNR
%rand('state',0);
BitsTx = floor(rand(1,N)*2);
SymQAMtmp = reshape(BitsTx,2,N/2).';
SymQAMtmptmp = bi2de(SymQAMtmp,2,'left-msb');
QAMTable = [1 1 -1 -1];
Ck1=QAMTable(SymQAMtmptmp+1);

r11=0;%定时偏移的均值
r22=0;%定时偏移的MSE
ffo22=0;%小数倍频偏的均值
ffo33=0;%小数倍频偏的MSE
ifo11=0;%发生整数倍频偏的概率
ifo22=0;%整数倍频偏的均值
ifo33=0;%整数倍频偏的MSE
fo11=0;%频偏的均值
fo22=0;%频偏的MSE
for k=1:Numloop
    k
rand('state',sum(100*clock));
a(1)=-1+2*round(rand(1,1));%+i*(-1+2*round(rand(1,1)));
b=-1+2*round(rand(1,N))+i*(-1+2*round(rand(1,N)));
r=floor(rand*(N-1));
info=zeros(1,N);
info(1)=a(1);
for l=2:N/2
a(l)=a(l-1)*Ck1(l-1);
info(2*l-1) =sqrt(2)*a(l);
end
Ck1(N/4)=a(1)/a(N/4);
carriers1=N.*ifft(info,N);
carriers2=N.*ifft(b,N);
carrierscp=zeros(1,(FS+N));
carrierscp(1:r)=carriers2(1:r);
carrierscp((r+1):(M+r))=carriers1((N/2-M+1):N/2);
carrierscp((r+M+1):(FS+r))=carriers1;
carrierscp((FS+r+1):(FS+N))=carriers2((1+r):N);
%Upconverter
L = length(carrierscp);
chips = [ carrierscp;zeros((4*q)-1,L)];
p=1/Rs:1/Rs:T;
g=ones(length(p),1);
dummy=conv(g,chips(:));
u=[dummy; zeros(45,1)];
[b,aa] = butter(13,1/20);
uoft = filter(b,aa,u);
delay=64; %Reconstruction filter delay
s_tilde=(uoft(delay+(1:length(t))).').*exp(1i*2*pi*fc*t);
s=real(s_tilde);
%OFDM RECEPTION
Ch=zeros(1,705);
Ch(1)=1;Ch(5)=0.81;Ch(9)=0.68;Ch(17)=0.55;Ch(29)=0.45;Ch(41)=0.37;Ch(65)=0.3;
Ch(77)=0.25;Ch(97)=0.2;Ch(129)=0.17;Ch(173)=0.34;Ch(225)=0.18;Ch(285)=0.1;
Ch(353)=0.06;Ch(429)=0.04;Ch(513)=0.02;Ch(605)=0.01;Ch(705)=0.002;
s1 = filter(Ch,1,s);
s_tidle = awgn(s1,SNR,'measured');
%Downconversion
f1=rand*(N-2);%2-1;
f2=floor(f1/2)*2;;%rand*(N/2-1);
f11=f1/Tu;
r_tilde1=exp(-1i*2*pi*fc*t).*s_tidle; %(F)
r_tilde=r_tilde1.*exp(1i*2*pi*f11*t);
%Carrier suppression
[B,AA] = butter(3,1/2);
r_info=2*filter(B,AA,r_tilde); %Baseband signal continuous-time (G)
%Sampling
r_datacp=real(r_info(1:(4*q):length(t)))... %Baseband signal, discretetime
+1i*imag(r_info(1:(4*q):length(t))); % (H)
S=zeros(1,length(r_datacp)-N+1);
for n=1:(length(r_datacp)-N+1)
    r_datacp1=r_datacp(n:(n+N-1));
    rtmp2(1:N/2)=r_datacp1((2):(N/2+1));
    for l=1:N/2
        rtmptmp2(l)=rtmp2(N/2-l+1);
    end
    rtmp3=r_datacp1((1+N/2):N);
    S(n)=rtmptmp2*conj(rtmp3)';
    P(n)=rtmp3*rtmp3';
    S1(n)=abs(S(n))/(P(n));
    S1(n)=S1(n)*S1(n);
end
[numtmp11 numtmp12]=max(S1);
%S1(numtmp12)=0;
%[numtmp21 numtmp22]=max(S1);
%numtmp12=min(numtmp12,numtmp22);
r1=numtmp12-(M+r+1);
r11=r11+r1;
r1=r1*r1;
r22=r1+r22;
%ffo
%numtmp12=(M+r+1);
    r_datacp1=r_datacp(numtmp12:(numtmp12+N-1));
    rtmp=r_datacp1(1:N/2);
    rtmptmp=r_datacp1((N/2+1):N);
    Ss=rtmptmp*rtmp';
ffo=angle(Ss)/(pi);   
if ffo<=0
    f2=f2+2;
end
f22=f1-f2;
ffo2=ffo-f22;
ffo22=ffo2+ffo22;
ffo2=ffo2*conj(ffo2);
Err1(k)=ffo2;
ffo33=ffo2+ffo33;

%ifo
r_tilde=r_tilde.*exp(-1i*2*pi*ffo/Tu*t);
[B,AA] = butter(3,1/2);
r_info=2*filter(B,AA,r_tilde);
r_datacp=real(r_info(1:(4*q):length(t)))... %Baseband signal, discretetime
+1i*imag(r_info(1:(4*q):length(t)));
r_data(1:N)=r_datacp(numtmp12:(N+numtmp12-1));
info_2N=(1/N).*fft(r_data,N);%fft
y1=zeros((N/2),1);
y2=zeros((N/2),1);
%y=zeros((N/2),1);
for l=1:(N/2-1)
    y1(l)=info_2N(2*l-1);
    y2(l)=info_2N(2*l+1);
    y(l)=y2(l)/y1(l);
end
y1(N/2)=y2(N/2-1);
y2(N/2)=y1(1);
%y(N/2)=info_2N(1)/info_2N(N-1);
for c=0:(N/2-1);
    for l=(c+1):(N/2)%(c+1+N/2-2)
        C(l-c)=y1(l)*conj(y2(l));%(l-floor((l)/(N/2-1))*(N/2-1));
    end     
    for l=(N/2+1):(c+1+N/2-1)
        C(l-c)=y1(l-(N/2))*conj(y2(l-(N/2)));
    end
    C=conj(C);
    J(c+1)=abs(Ck1*(C'));
end
[numtmp31 numtmp32]=max(J);
ifo=2*(numtmp32-1);
if ifo~=f2
    ifo11=ifo11+1;
end
ifo2=ifo-f2;
Err(k)=ifo2;
ifo22=ifo22+ifo2;
ifo2=ifo2*conj(ifo2);
ifo33=ifo33+ifo2;
fo=ifo+ffo-f1;
fo11=fo+fo11;
fo=fo*conj(fo);
fo22=fo+fo22;

end
r11=r11/Numloop;%定时偏移的均值
r22=r22/Numloop;%定时偏移的MSE
ffo22=ffo22/Numloop;%小数倍频偏的均值
ffo33=ffo33/Numloop;%小数倍频偏的MSE
ifo11=ifo11/Numloop;%发生整数倍频偏的概率
ifo22=ifo22/Numloop;%整数倍频偏的均值
ifo33=ifo33/Numloop;%整数倍频偏的MSE
fo11=fo11/Numloop;%频偏的均值
fo22=fo22/Numloop%频偏的MSE
ErrSym(:,SNR/5+1)=[r11;r22;ffo22;ffo33;ifo11;ifo22;ifo33;fo11;fo22];
end

n=1:(length(r_datacp)-N+1);
figure(5);
%subplot(211);
plot(n,S1);
%subplot(212);
%plot(n,MM);