% ofdm.m
% 實現OFDM傳輸的仿真程序
function [ber] = ofdm1(ebn0,nloop)
%********************** 準備部分 ***************************

para=256;   % 并行子信道的個數
fftlen=256; % FFT的長度
noc=256;    % 載波的個數（這里我認為其數目和子信道數目相同）
nd=16;       % 循環一次OFDM的符號個數
ml=2;       % 采用QPSK調制
sr=250000;  % 符號速率
br=sr.*ml;  % 每個載波的比特速率
gilen=40;   % 保護間隔的長度


%************************** 主循環部分 **************************

% nloop=5000;  % 仿真循環次數
noe = 0;    % 錯誤數據的個數
nod = 0;    % 傳輸數據的個數
eop=0;      % 錯誤包的個數
nop=0;      % 傳輸包的個數

 for iii=1:nloop

%****************** 發射部分 ******************************************************

%****************** 產生數據 **********************

%********destination USER

numcarrier=8;
numuser=para/numcarrier;
userdata=rand(1,numcarrier*nd*ml)>0.5;           %  rand : built in function
userpara=reshape(userdata,numcarrier,nd*ml);     %  reshape : built in function

%********other users
paradata=zeros(para,nd*ml);% 初始化進行IFFT矩陣
for i=1:numcarrier%把IFFT矩陣的每隔numuser列設為目標用戶的信息
    paradata((1+(i-1)*numuser),:)=userpara(i,:);
end
for ii=1:(numuser-1)%在IFFT矩陣中加入其他用戶信息
    %隨機產生其他某個用戶的信息并加入IFFT矩陣矩陣
    disturbdata=rand(1,numcarrier*nd*ml)>0.5;
    disturbpara=reshape(disturbdata,numcarrier,nd*ml); 
    for jj=1:numcarrier%把IFFT矩陣的每隔numuser列設為干擾用戶的信息
    paradata(((ii+1)+(jj-1)*numuser),:)=disturbpara(jj,:);
    end

end

%****************** QPSK調制 ***********************
paradata=paradata.*2-1;
inputdata=zeros(para,nd);
count=0;
j=sqrt(-1);
for jj=1:nd
    inputdata(:,jj)=paradata(:,(count+1))+j.*paradata(:,(count+2));
    count=count+2;
%     count
end
%****************** IFFT(離散傅里葉反變換) *********
y=ifft(inputdata);
%****************** 插入保護間隔 *******************
y=reshape(y,1,para*nd);
ty=zeros(1,(length(y)+nd*gilen));
for ii=1:nd
    a=para+gilen;
    ty(((ii-1)*a+gilen+1):ii*a)=y(((ii-1)*para+1):ii*para);
    ty(((ii-1)*a+1):((ii-1)*a+gilen))=y((ii*para-gilen+1):ii*para);
end

spow=sum(abs(ty).^2)/(nd*para);        %  sum : built in function
attn=0.5*spow*(10^(-ebn0/10));
attn=sqrt(attn);


%****************** 信道 ***************************************************** 
[y1]=comb1(ty,attn);

%****************** 接收部分  ******************************************************

%****************** 去除保護間隔 ********************
y2=zeros(1,length(y));
for ii=1:nd
    a=para+gilen;
    y2(((ii-1)*para+1):ii*para)=y1(((ii-1)*a+gilen+1):ii*a);
%     y2((ii*para-gilen+1):ii*para)=y1(((ii-1)*a+1):((ii-1)*a+gilen));
end
ry=reshape(y2,para,nd);

%****************** FFT(離散傅里葉變換)  ************
outy=fft(ry);

%***************** QPSK解調 **************************
demodata=zeros(para,nd*ml);
demodata(:,1:2:(nd*ml-1))=0.5*(outy+conj(outy))>0;
demodata(:,2:2:nd*ml)=0.5*j*(-outy+conj(outy))>0;



%**************  并串變換  ***************************


reuserpara=zeros(numcarrier,nd*ml);
for ii=1:numcarrier%把IFFT矩陣的每隔numuser列設為目標用戶的信息
    reuserpara(ii,:)=demodata((1+(ii-1)*numuser),:);
end
reuserdata=reshape(reuserpara,1,numcarrier*nd*ml);




%*************** 比特誤碼率 (BER) *********************

% instantaneous number of error and data
eeee=0;
noe2=0;
nod2=0;

noe2=sum(abs(reuserdata-userdata));         %  sum : built in function
nod2=length(userdata);  %  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

%****************** 結果輸出 ***************************
% 
% 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 ******************************************************