ofdmoffset_sc_dx.m

这是一个完整的ofdm系统仿真matlab程序

% Program 4-1
% ofdm.m
%
% Simulation program to realize OFDM transmission system
%
% programmed by T.Yamamura and H.Harada
%
%
%仿真SAC算法 （定时与小数倍频偏）
function [derr,ferr]=ofdmoffset_sc_dx(ebn0,gilen)
%********************** preparation part ***************************

para=128;   % Number of parallel channel to transmit (points)
fftlen=128; % FFT length
noc=128;    % 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符号总长
%ebn0=10;     %Eb/N0
snr=10^(ebn0/10);
deltad=20;  %时延
deltaf=2.45;  %频偏
deltaderr=0;%时延均方误差
deltaferr=0;%频偏均方误差
%f=-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;
    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
        
 
    
    

%b=0.9*
%m=find(abs(p-b)<0.04) ; % 求最大值所在位置find(x< a & x> 0.8)
%m1=m;
%m1=find((m1-mean(m))<0);
%temp1=mean(m(1:size(m1)));
%temp2=mean(m(size(m1)+1:end));
%temp=floor((temp1+temp2)/2);
%b1=find(p==0.9*a) ; 
%b1=find(p< 0.901*a & p> 0.9*a);
%b2=sum(b1);
%b3=length(b1);
%m=b2./b3;
%b1
f=-angle(p1(j5))/pi;
%f
%p
% r
 %m
%这是一种 最大似然方法联合实现符号定时同步和载波同步的方法
% 具体的推导 见 《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);
f
%********************频偏的估计值***********************
l=g+f; %频偏估计值

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

deltaderr=deltaderr+(j5-16-deltad)^2;
deltaferr=deltaferr+(f-0.45)^2;
    
end

% derr=sqrt(deltaderr/(nloop*para*nd*ml));
% ferr=sqrt(deltaferr/(nloop*para*nd*ml));

derr=sqrt(deltaderr/(nloop));
ferr=sqrt(deltaferr/(nloop));
%derr=deltaderr/(nloop);
%ferr=deltaferr/(nloop);
%ferr=sqrt(deltaferr/(nloop));



% %****************** Guard interval removal *********
% 
% [ich5,qch5]= girem(ich4,qch4,fftlen2,gilen,nd);
% 
% %******************  FFT  ******************
% 
% rx=ich5+qch5.*i;
% ry=fft(rx);   	% fft : built in function
% ich6=real(ry);	% real : built in function
% qch6=imag(ry);	% imag : built in function
% 
% %***************** demoduration *******************
% 
% ich7=ich6./kmod;
% qch7=qch6./kmod;
% [demodata]=qpskdemod(ich7,qch7,para,nd,ml);   
% 
% %**************  Parallel to serial conversion  *****************
% 
% demodata1=reshape(demodata,1,para*nd*ml);
% 
% % %************************** Bit Error Rate (BER) ****************************
% 
% % instantaneous number of error and data
% 
% noe2=sum(abs(demodata1-seldata));  %  sum : built in function
% nod2=length(seldata);  %  length : built in function
% 
% % cumulative the number of error and data in noe and nod
% 
% noe=noe+noe2;
% nod=nod+nod2;
% 
% % calculating PER
% 
% if noe2~=0  
%    eop=eop+1;
% else
%    eop=eop;
% end   
%    eop;
%    nop=nop+1;
%    
% 
% fprintf('%d\t%e\t%d\n',iii,noe2/nod2,eop);  %  fprintf : built in function
%    
% end
% 
% %********************** Output result ***************************
% 
% per=eop/nop;
% ber=noe/nod;
% 
% fprintf('%f\t%e\t%e\t%d\t\n',ebn0,ber,per,nloop);
% fid = fopen('BERofdm.dat','a');
% fprintf(fid,'%f\t%e\t%e\t%d\t\n',ebn0,ber,per,nloop);
% fclose(fid);
% 
% %******************** end of file ***************************