ofdmtest.m

一个很好的OFDM系统的测试程序

% 
% ofdm.m
%
% Simulation program to obtain Subcarrier Hopping OFDM transmission system
%
% programmed by Huang Yong
%
clear
clc
%********************** preparation part ***************************
noc = 8;         % Number of carrier  载波数为8



%********************** preparation part for muti-users *************
para = 8;                      % Number of parallel channel to transmit (point) for each user  信道数量
usernum = 8;                   % Number of the users
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
fftlen = usernum*para;         % FFT length   FFT的长度
gilen = usernum*para/4;        % length of guard interval (points)  循环前缀的长度

%********************* fh code generated part ***********************
%fhcode = [
%    1:8
%    2:8,1
%    3:8,1:2
%    4:8,1:3
%    5:8,1:4
%    6:8,1:5
%];
%fhcode = fhcode';
fhcode1 = monehit(3,[1,3],[0,0,1],2,3);
fhcode = fhcode1(:,1:nd);
fhcode = fhcode+1;            % Because the bound of the fhcode1 is 0~7
                              % now convert to 1~8
%********************************************************************

%******************** fading initialization *************************

tstp = 1/sr/(fftlen+gilen);                 % Time resolution
itau = [0];                                 % Arrival time for each multipath normalized by tstp
dlvl1 = [0];                                % Mean power for each multipath normalized by direct wave
n0 = [6];                                   % Number of waves to generate fading n0(1), n0(2)
th1 = [0.0];                                % Initial Phase of delayed wave
itnd1 = [1000];                             % Set fading counter
now1 = 1;                                   % Number of directwave + Number of delayed wave
fd = 150;                                   % Maximum Doppler frequency
flat = 1;                                   % Flat or not
itnd0 = nd*(fftlen+gilen)*20;               % Number of fading counter to skip
%********************************************************************

ebn0 = 0:2:50;
ber = zeros(1,21);
lll = length(ebn0);
for jjj = 1:lll
%********************* 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


    ichd = zeros(para,nd);                  % Preparation for fft
    qchd = zeros(para,nd);                  % Preparation for fft
    demodd = zeros(para*usernum,nd*ml);     % Preparation for demodulation

    for iii = 1:nloop
%******************** transmitter ***********************************
            ichmod = zeros(1,(fftlen+gilen)*nd);    % Preparation for ifft
            qchmod = zeros(1,(fftlen+gilen)*nd);    % Preparation for ifft

    
%******************** Data generation for each user *****************
%//////////////////////////////////////////////////////////////////////////

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

%******************** Serial to parallel conversion ********************
        paradato = reshape(seldata,para*usernum,nd*ml);    % reshape: built in system function

%******************** assign subcarriers for each user *****************
        for xx = 1:usernum               % One loop represent one user's subcarrier assignment 

            for temp0 = 1:para
                paradata(temp0,:) = paradato(para*(xx-1)+temp0,:);
            end

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

%******************** QPSK modulation **********************************
            [ich,qch] = qpskmod(paradata,para,nd,ml);
            kmod = 1/sqrt(2);                            % sqrt: built in system function
            ich1 = ich.*kmod;
            qch1 = qch.*kmod;
%******************** Subcarriers hopping ******************************
            for ii = 1:nd
                ich2(fhcode(:,ii),ii) = ich1(:,ii);
                qch2(fhcode(:,ii),ii) = qch1(:,ii);
            end

%******************** fading channel ***********************************
% In this program, assume each subcarrier is faded  by mutipath fading
%for nnn=1:para
%    itemp = ich1(nnn,:);
%    qtemp = qch1(nnn,:);
%[ifade(nnn,:),qfade(nnn,:),ramp,rcos,rsin] = sefade(itemp,qtemp,itau,dlvl1,th1,n0,itnd1,now1,nd,tstp,fd,flat);
%itnd1 = itnd1+itnd0;                          % Updata fading counter
%ich1(nnn,:) = ifade(nnn,:);
%qch1(nnn,:) = qfade(nnn,:);
%end
%***********************************************************************

            for nnn = 1:para
                [ifade(nnn,:),qfade(nnn,:)] = myfading(ich2(nnn,:),qch2(nnn,:),nd);
            end
            ich1 = ifade;
            qch1 = qfade;
%//////////////////////////////////////////////////////////////////////////

%******************* Data modulation for each user **********************
    ichc = zeros(para*usernum,nd);          % Preparation for ifft of each user
    qchc = zeros(para*usernum,nd);          % Preparation for ifft of each user
            for  temp1 = 1:para
                 ichc(para*(xx-1)+temp1,:) = ich1(temp1,:);
                 qchc(para*(xx-1)+temp1,:) = qch1(temp1,:);
            end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%******************** IFFT *********************************************

        x = ichc + qchc.*i;
        y = ifft(x);                                 % ifft: built in system function
        ich2 = real(y);                              % real: built in system function
        qch2 = imag(y);                              % imag: built in system function

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

        [ich3,qch3] = giins(ich2,qch2,fftlen,gilen,nd);  % Insert gurad interval
        ichmod = ichmod + ich3;
        qchmod = qchmod + qch3;
        
        end
  fftlen2 = fftlen +gilen;


%******************** IFFT *********************************************
%
%        x = ichc + qchc.*i;
%        y = ifft(x);                                 % ifft: built in system function
%        ich2 = real(y);                              % real: built in system function
%        qch2 = imag(y);                              % imag: built in system function
%
%******************** Gurad interval insertion *************************
%
%        [ich3,qch3] = giins(ich2,qch2,fftlen,gilen,nd);  % Insert gurad interval
%        fftlen2 = fftlen +gilen;                         % fftlen2 is the length of the symbol inserted gurad interval


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

        spow = sum(ichmod.^2+qchmod.^2)/nd./para;            % sum: built in system function

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

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

        [ich4,qch4] = comb(ichmod,qchmod,attn);              % By using "attn" and comb.m,the transmitted data was contaminated by AWGN

%******************* Receiver *******************************************
%******************* Guard interval removal *****************************

        [ich5,qch5] = girem(ich4,qch4,fftlen2,gilen,nd);

%******************** FFT ************************************************

        rx = ich5+qch5.*i;
        ry = fft(rx);                                    % fft: built in system function
        ich6 = real(ry);                                 % real: built in system function
        qch6 = imag(ry);                                 % imag: built in system function

%******************* demodulation ****************************************

        ich7 = ich6./kmod;                               % Unnormalize the data
        qch7 = qch6./kmod;
%****************** mutiuser convert to each user ************************
        for xx = 1:usernum                              % One loop represent one user's subcarrier recover and demodulation
            for temp2 = 1:para
                ichd(temp2,:) = ich7(para*(xx-1)+temp2,:);
                qchd(temp2,:) = qch7(para*(xx-1)+temp2,:);
            end
%******************* subcarriers defh ************************************
            for ii = 1:nd
                ich8(:,ii) = ichd(fhcode(:,ii),ii);
                qch8(:,ii) = qchd(fhcode(:,ii),ii);
            end
            [demodata] = qpskdemod(ich8,qch8,para,nd,ml);
%**************************************************************************
            for temp3 = 1:para
                demodd(para*(xx-1)+temp3,:) = demodata(temp3,:); 
            end
        end
%******************* Parallel to serial conversion ************************

        demodata1 = reshape(demodd,1,para*usernum*nd*ml);


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

%******************* Bit Error Rate (BER) *********************************
% instantaneous number of errors and data bits

        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;

    end

%*********************** Output result ************************************

    per = eop/nop;
    ber(jjj) = noe/nod;
    fprintf('ebn0=%d\t ber=%e\t\n',ebn0(jjj),ber(jjj));   % fprintf: built in system function
    fprintf('------------------------------------------------------\n');
end

%**************************************************************************
SNRindB2=0:0.1:50;
for i=1:length(SNRindB2)
    SNR=exp(SNRindB2(i)*log(10)/10);
    theo_err_prb(i)=0.5*(1-1/sqrt(1+1/SNR));
end
semilogy(ebn0,ber,'-*r');
hold
semilogy(SNRindB2,theo_err_prb);
grid on

%************************* end of file ************************************