ofdmtmep.m

非常好用

%在同一个符号中选择子载波作导频,然后利用Linear或DFT插植
%Channel :  AWGN  ,Taxi  ,Arrival,En-route,Parking
%由于各部分大同小异 ,仅在Taxi1部分给出详细注释
%所有的信道仅对ChannelAeroT1做了详细的注释 
%ATTENTION!=由于后期经过修改,程序中,NumSubc指的事载有信息的子载波数,
%而NumSubc+Numpps为子载波总数,其中Numpps为导频数(需为128的约数加1)
%修改纪录: 
%(1)2008-01-21 将SNR步进改为2dB 
%(2)2008-02-28 $$$$1提供的减轻ICI和噪声的方法,启用标志sigLp
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tic;
clear,clc;
format short;
Nn=50;                          %num of circle
Ns=1;                          %num of  symbols per exp
Numpps=16+1;                    %Num_pilot_per_symbol
NF=(256-2*Numpps)*Ns;           %Num of inf bits
NumCp=16;                        %Num of circle prefix 

DatahT1=zeros(Nn,16);           %保存各种实验结果的缓冲区
DatahT2=zeros(Nn,16);
DatahP=zeros(Nn,16);
DatahAn=zeros(Nn,16);
DatahE=zeros(Nn,16);
DatahAr=zeros(Nn,16);

Taxi1=1;                        %控制各种实验的标志,'1'开启实验.'0'禁止试验
Taxi2=0;
Parking=0;
Awgn=0;
Arrival=0;
En_route=0;

sigLp=0;                        %估计值低通滤波开通信号
sigdft=0;                       %DFT插值法开通信号,0是使用线形插值
sigcom=0;                       %using or not compensation in en-route and arrival 
file_name='rayleigh.txt';  %the name of the file used to save the result data
com=0;
comwyb=0;
num_delay=0;

for haha=1:Nn;
    OP=haha

if Taxi1==1    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Taxi1  situation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for SNR=0:2:30
   
    
N=NF;
Data=floor(rand(1,N)*2);                %data source

Fh=492000;
T=1/Fh;
NumSubc=128-Numpps;                     %Num of subcarriers with data
NumQ=2;                                 %num of bits per subc
NumP=NumSubc*NumQ;                      %num of bits per symbol 
NumError=0;                             %counter of error bits 
temptemp=(NumSubc+Numpps)/(Numpps-1);   %temptemp is the num of subc between two pilots

NumFrame=ceil(N/NumP);                  %divided into symbol and adding zeros if needed
NumData=NumFrame*NumP;
DataC(1:N)=Data(1:N);
DataC(N+1:NumData)=0;

QAMTable=[-1-j -1+j 1-j 1+j];           %Table of Qpsk 

DataCF=reshape(DataC,NumP,NumFrame);    %divided into symbol
DataSer=zeros(NumSubc+Numpps+NumCp,NumFrame);      %DataSer is the data buffer for transmitting
DataSerR=zeros(NumP,NumFrame);                     %DataSerR if the data buffer for receiving

for ii=1:NumFrame
    DataT(1:NumP)=DataCF(1:NumP,ii);
    DataTtemp1=reshape(DataT,NumSubc,NumQ);        %interlever in the same time,haha! 
    DataTtemp2=bi2de(DataTtemp1,'left-msb');       %converse bi to de
    DataTtemp=QAMTable(DataTtemp2+1);              %modulation via look-up table
  
    DataTtemp=[DataTtemp,1+j];                     %the No.128 Subc is pilot
    DataTtemptemp=zeros(NumSubc+Numpps,1);         %inter pilots 
    for jj=1:Numpps-1
        stemp=(jj-1)*temptemp+1;
        DataTtemptemp(stemp)=1+j;
        DataTtemptemp(stemp+1:stemp+temptemp-1)=DataTtemp(stemp-jj+1:stemp-jj+temptemp-1);
    end

    DataIFFT=ifft(DataTtemptemp,NumSubc+Numpps);   %OFDM via IFFT
    DataIFFTCp(1:NumCp)=DataIFFT(NumSubc+Numpps-NumCp+1:NumSubc+Numpps);   %adding circle prefix
    %DataIFFTCp(1:NumCp)=0;
    DataIFFTCp(NumCp+1:NumCp+NumSubc+Numpps)=DataIFFT(1:NumSubc+Numpps);
    DataSer( : ,ii)=DataIFFTCp;       
end

%%%Passing the channel%%%%%%%%%%%%%%%%%%%%%%%%%
p=(NumSubc+Numpps+NumCp)*NumFrame;                  %toltal num of sampling point
DataSerShape=reshape(DataSer,1,p);                  %becoming serial  data
data=DataSerShape;
%chan=rayleighchan(T,[250],[0.6e-15 0.7e-15],[-3 -3]);
%chan.NormalizePathGains=1;
%chan.StoreHistory=1;
%temp=filter(chan,[zeros(1,10),data]);
%DataChannel=temp(11:length(data)+10);
%DopCom=1;
%mean(data.*data),
%mean(DataChannel.*DataChannel)
%df=sdlfkjs
[DataChannel,DopCom]=ChannelAeroT1(data);            %passing time-variable channel

temp=10.^(SNR/10);                %converse dB to linear
SNRr=10*log10(temp*2);
DataChannel=awgn([DataChannel],SNRr,'measured');     %adding Gauss Noise

%data=DataChannel;
%[DataChannel,DopCom]=ChannelAeroT1_com(data);  
%data=data-DataChannel;
%[DataChannel,DopCom]=ChannelAeroT1_com(data);  
%data=data-DataChannel;
%[DataChannel,DopCom]=ChannelAeroT1_com(data);  
%data=data-DataChannel;
%DataChannel=data;

%[out,DopCom]=Syn(DataChannel,NumSubc+Numpps,NumCp,SNR); %estimate the doppler frenquency
DataChannel=(DopCom).*DataChannel;  

DataRtempSer=reshape(DataChannel,NumSubc+Numpps+NumCp,NumFrame);   %divided into symbol
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ii=1:NumFrame
    DataRtempCp=DataRtempSer( : ,ii);
    DataRtemp=DataRtempCp(NumCp+1:NumCp+NumSubc+Numpps);
    
    if comwyb==1
        
        DataRtemp_temp=DataRtemp;
        DataRtemp=delay_fun(DataRtemp_temp,num_delay,NumCp);
    end
    
    DataRtemp1=fft(DataRtemp,NumSubc+Numpps); 
    
    t=1:temptemp:NumSubc+Numpps;                      %linear inter base pilot 
    t=[t,NumSubc+Numpps];
    ffttemp(1:Numpps)=(DataRtemp1(t)./(1+j));
   
   
    %$$$$1提供的减轻ICI和噪声的方法
    if sigLp==1;    
    temp=fft(ffttemp);
    temp(2:16)=0;
    ffttemp=ifft(temp);
    end
    
    if com==1
    if sigdft==0;
    [ph]=linear_inter1(ffttemp,NumSubc+Numpps);
    else 
        temp=ifft(ffttemp);
        ph=fft(temp,NumSubc+Numpps);
    end
    
    Comtemp=ph;            
    DataRtemp1=DataRtemp1./(Comtemp.');
    end
        
    for jj=1:Numpps-1                                 %remove the pilot
        stemp=(jj-1)*temptemp+1;        
        DataRtemp3temp(stemp-jj+1:stemp-jj+temptemp-1)= DataRtemp1(stemp+1:stemp+temptemp-1);
    end
    DataRtemp1=DataRtemp3temp(1:NumSubc);
    
    DataRtemp3( : ,1)=real(DataRtemp1)>0;               %making Qpsk decision
    DataRtemp3( : ,2)=imag(DataRtemp1)>0;  
    DataSerR( : ,ii)=reshape(DataRtemp3,1,NumQ*(NumSubc));
end
DataR=reshape(DataSerR,1,NumData);
re=DataC-DataR;
NumError=sum(abs(re));
ErrorN(SNR/2+1)=NumError;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end
DatahT1(haha, : )=DatahT1(haha, : )+ErrorN;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end

if Taxi2==1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Taxi2  situation%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for SNR=0:2:30
    
N=NF;
Data=floor(rand(1,N)*2);               %data source

Fh=984000;
T=1/Fh;
NumSubc=128-Numpps;                    %Num of subcarriers with data
NumQ=2;                                 %num of bits per subc
NumP=NumSubc*NumQ;
NumError=0;
temptemp=(NumSubc+Numpps)/(Numpps-1);

NumFrame=ceil(N/NumP);                  %divided into symbol and adding zeros if needed
NumData=NumFrame*NumP;
DataC(1:N)=Data(1:N);
DataC(N+1:NumData)=0;
QAMTable=[-1-j -1+j 1-j 1+j];

DataCF=reshape(DataC,NumP,NumFrame);
DataSer=zeros(NumSubc+Numpps+NumCp,NumFrame);
DataSerR=zeros(NumP,NumFrame);

for ii=1:NumFrame
    DataT(1:NumP)=DataCF(1:NumP,ii);
    DataTtemp1=reshape(DataT,NumSubc,NumQ);        %interlever in the same time,haha! 
    DataTtemp2=bi2de(DataTtemp1,'left-msb');
    DataTtemp=QAMTable(DataTtemp2+1);
   
    DataTtemp=[DataTtemp,1+j];
    DataTtemptemp=zeros(NumSubc+Numpps,1);       %inter pilots
    for jj=1:Numpps-1
        stemp=(jj-1)*temptemp+1;
        DataTtemptemp(stemp)=1+j;
        DataTtemptemp(stemp+1:stemp+temptemp-1)=DataTtemp(stemp-jj+1:stemp-jj+temptemp-1);
    end
    
    DataIFFT=ifft(DataTtemptemp,NumSubc+Numpps);  
    DataIFFTCp(1:NumCp)=DataIFFT(NumSubc+Numpps-NumCp+1:NumSubc+Numpps);
    DataIFFTCp(NumCp+1:NumCp+NumSubc+Numpps)=DataIFFT(1:NumSubc+Numpps);
    DataSer( : ,ii)=DataIFFTCp;       
end

p=(NumSubc+Numpps+NumCp)*NumFrame;
DataSerShape=reshape(DataSer,1,p);
data=DataSerShape;
[DataChannel,DopCom]=ChannelAeroT2(data); 
temp=10.^(SNR/10);
SNRr=10*log10(temp*2);
DataChannel=awgn(DataChannel,SNRr,'measured');
DataChannel=(DopCom).*DataChannel;    
DataRtempSer=reshape(DataChannel,NumSubc+Numpps+NumCp,NumFrame);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ii=1:NumFrame
    DataRtempCp=DataRtempSer( : ,ii);
    DataRtemp=DataRtempCp(NumCp+1:NumCp+NumSubc+Numpps);
    DataRtemp1=fft(DataRtemp,NumSubc+Numpps); 
    
    t=1:temptemp:NumSubc+Numpps;
    t=[t,NumSubc+Numpps];
    ffttemp(1:Numpps)=(DataRtemp1(t)./(1+j));
    
    t=1:temptemp:NumSubc+Numpps;
    t=[t,NumSubc+Numpps];
    ffttemp(1:Numpps)=(DataRtemp1(t)./(1+j));
    
    %$$$$1提供的减轻ICI和噪声的方法
    if sigLp==1;    
    temp=fft(ffttemp);
    temp(2:16)=0;
    ffttemp=ifft(temp);
    end
    
    if sigdft==0;
    [ph]=linear_inter1(ffttemp,NumSubc+Numpps);
    else 
        temp=ifft(ffttemp);
        ph=fft(temp,NumSubc+Numpps);
    end
    
    Comtemp=ph;            
    DataRtemp1=DataRtemp1./(Comtemp.');
        
    for jj=1:Numpps-1
        stemp=(jj-1)*temptemp+1;        
        DataRtemp3temp(stemp-jj+1:stemp-jj+temptemp-1)= DataRtemp1(stemp+1:stemp+temptemp-1);
    end
    DataRtemp1=DataRtemp3temp(1:NumSubc);
    
    DataRtemp3( : ,1)=real(DataRtemp1)>0;
    DataRtemp3( : ,2)=imag(DataRtemp1)>0;  
    DataSerR( : ,ii)=reshape(DataRtemp3,1,NumQ*(NumSubc));
end
DataR=reshape(DataSerR,1,NumData);
re=DataC-DataR;
NumError=sum(abs(re));
ErrorN(SNR/2+1)=NumError;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end
DatahT2(haha, : )=DatahT2(haha, : )+ErrorN;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end

if Parking==1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Parking
for SNR=0:2:30
    
N=NF;
Data=floor(rand(1,N)*2);               %data source

Fh=984000;
T=1/Fh;
NumSubc=128-Numpps;                    %Num of subcarriers with data
NumQ=2;                                 %num of bits per subc
NumP=NumSubc*NumQ;
NumError=0;
temptemp=(NumSubc+Numpps)/(Numpps-1);

NumFrame=ceil(N/NumP);                  %divided into symbol and adding zeros if needed
NumData=NumFrame*NumP;
DataC(1:N)=Data(1:N);
DataC(N+1:NumData)=0;
QAMTable=[-1-j -1+j 1-j 1+j];

DataCF=reshape(DataC,NumP,NumFrame);
DataSer=zeros(NumSubc+Numpps+NumCp,NumFrame);
DataSerR=zeros(NumP,NumFrame);

for ii=1:NumFrame
    DataT(1:NumP)=DataCF(1:NumP,ii);
    DataTtemp1=reshape(DataT,NumSubc,NumQ);        %interlever in the same time,haha! 
    DataTtemp2=bi2de(DataTtemp1,'left-msb');
    DataTtemp=QAMTable(DataTtemp2+1);
    
    DataTtemp=[DataTtemp,1+j];
    DataTtemptemp=zeros(NumSubc+Numpps,1);       %inter pilots
    for jj=1:Numpps-1
        stemp=(jj-1)*temptemp+1;
        DataTtemptemp(stemp)=1+j;
        DataTtemptemp(stemp+1:stemp+temptemp-1)=DataTtemp(stemp-jj+1:stemp-jj+temptemp-1);
    end
      
    DataIFFT=ifft(DataTtemptemp,NumSubc+Numpps);  
    
    DataIFFTCp(1:NumCp)=DataIFFT(NumSubc+Numpps-NumCp+1:NumSubc+Numpps);
    DataIFFTCp(NumCp+1:NumCp+NumSubc+Numpps)=DataIFFT(1:NumSubc+Numpps);
    DataSer( : ,ii)=DataIFFTCp;       
end

p=(NumSubc+Numpps+NumCp)*NumFrame;
DataSerShape=reshape(DataSer,1,p);
data=DataSerShape;
[DataChannel,DopCom]=ChannelAeroP(data); 
temp=10.^(SNR/10);
SNRr=10*log10(temp*2);
DataChannel=awgn(DataChannel,SNRr,'measured');
%[out,DopCom]=Syn(DataChannel,NumSubc+Numpps,NumCp,SNR);
DataChannel=(DopCom).*DataChannel;   
DataRtempSer=reshape(DataChannel,NumSubc+Numpps+NumCp,NumFrame);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ii=1:NumFrame
    DataRtempCp=DataRtempSer( : ,ii);
    DataRtemp=DataRtempCp(NumCp+1:NumCp+NumSubc+Numpps);
    DataRtemp1=fft(DataRtemp,NumSubc+Numpps); 
    
    t=1:temptemp:NumSubc+Numpps;
    t=[t,NumSubc+Numpps];
    ffttemp(1:Numpps)=(DataRtemp1(t)./(1+j));
    
    t=1:temptemp:NumSubc+Numpps;
    t=[t,NumSubc+Numpps];
    ffttemp(1:Numpps)=(DataRtemp1(t)./(1+j));
    
    %$$$$1提供的减轻ICI和噪声的方法
    if sigLp==1;    
    temp=fft(ffttemp);
    temp(4:14)=0;
    ffttemp=ifft(temp);
    end
    
    if sigdft==0;
    [ph]=linear_inter1(ffttemp,NumSubc+Numpps);
    else