echo off;clear all;close all;clc;
fprintf( 'OFDM仿真\n') ;
tic
% --------------------------------------------- %
%                   參數(shù)定義                     %
% --------------------------------------------- %
% Initialize the parameters
NumLoop = 1000;
NumSubc = 128;
NumCP = 8;
SyncDelay = 0;
% 子載波數(shù)            128
% 位數(shù)/ 符號(hào)          2
% 符號(hào)數(shù)/ 載波        1000
% 訓(xùn)練符號(hào)數(shù)          0
% 循環(huán)前綴長(zhǎng)度        8   (1/16)*T
% 調(diào)制方式            4-QAM
% 多徑信道數(shù)          3 
% IFFT Size           128 
% 信道最大時(shí)延        2
% --------------------------------------------- %
%                   QAM MODULATION              %
% --------------------------------------------- %
% Generate the random binary stream for transmit test
BitsTx = floor(rand(1,NumLoop*NumSubc)*2);
% Modulate (Generates QAM symbols)
% input: BitsTx(1,NumLoop*NumSubc); output: SymQAM(NumLoop,NumSubc/2)
SymQAMtmp = reshape(BitsTx,2,NumLoop*NumSubc/2).';
SymQAMtmptmp = bi2de(SymQAMtmp,2,'left-msb');
%--------------------------------------------------------------------
% 函數(shù)說(shuō)明：
% bin2dec(binarystr) interprets the binary string binarystr and returns the
% equivalent decimal number.
%   bi2de是把列向量的每一個(gè)元素都由2進(jìn)制變?yōu)?0進(jìn)制
%  D = BI2DE(...,MSBFLAG) uses MSBFLAG to determine the input orientation.
%     MSBFLAG has two possible values, 'right-msb' and 'left-msb'.  Giving a
%     'right-msb' MSBFLAG does not change the function's default behavior.
%   Giving a 'left-msb' MSBFLAG flips the input orientation such that the
%     MSB is on the left.
% % %   D = BI2DE(...,P) converts a base P vector to a decimal value.
% %     Examples:
% %     >> B = [0 0 1 1; 1 0 1 0];
% %     >> T = [0 1 1; 2 1 0];
% %    >> D = bi2de(B)     >> D = bi2de(B,'left-msb')     >> D = bi2de(T,3)
% %     D =                 D =                            D =
% %         12                   3                             12
% %          5                  10                              5
%--------------------------------------------------------------------
% QAM modulation 
% 00->-1-i,01->-1+i,10->1-i,11->1+i
% 利用查表法進(jìn)行QAM星座映射
QAMTable = [-1-i -1+i 1-i 1+i];
SymQAM = QAMTable(SymQAMtmptmp+1);
% --------------------------------------------- %
%                   IFFT                        %
% --------------------------------------------- %
% input: SymQAM(NumLoop,NumSubc/2); output: SymIFFT(NumSubc,NumLoop)
SymIFFT = zeros(NumSubc,NumLoop);
SymIFFTtmp = reshape(SymQAM,NumSubc/2,NumLoop);
SymIFFTtmptmp = zeros(NumSubc,NumLoop);
SymIFFTtmptmp(1,:) = real(SymIFFTtmp(1,:)); % 實(shí)數(shù)
SymIFFTtmptmp(NumSubc/2+1,:) = imag(SymIFFTtmp(1,:)); % 實(shí)數(shù)
% 這么安排矩陣的目的是為了構(gòu)造共軛對(duì)稱(chēng)矩陣
% 共軛對(duì)稱(chēng)矩陣的特點(diǎn)是 在ifft/fft的矢量上 N點(diǎn)的矢量
% 在0，N/2點(diǎn)必須是實(shí)數(shù) 一般選為0
% 1至N/2點(diǎn) 與 （N/2）+1至N-1點(diǎn)關(guān)于N/2共軛對(duì)稱(chēng)
SymIFFTtmptmp(2:NumSubc/2,:) = SymIFFTtmp(2:NumSubc/2,:);
SymIFFTtmptmp((NumSubc/2+2):NumSubc,:) = flipdim(conj(SymIFFTtmp(2:NumSubc/2,:)),1);
%--------------------------------------------------------------------
% 函數(shù)說(shuō)明：
% B = flipdim(A,dim) returns A with dimension dim flipped.
% When the value of dim is 1, the array is flipped row-wise down. When dim is 2,
% the array is flipped columnwise left to right. flipdim(A,1) is the same as
% flipud(A), and flipdim(A,2) is the same as fliplr(A).
%--------------------------------------------------------------------
% %  >> a = [1 2 3; 4 5 6; 7 8 9; 10 11 12]
% %  a =
% %      1     2     3
% %      4     5     6
% %      7     8     9
% %     10    11    12
% %  >> b = flipdim(a,1)
% %  b =
% %     10    11    12
% %      7     8     9
% %      4     5     6
% %      1     2     3
SymIFFT = ifft(SymIFFTtmptmp,NumSubc,1);
% --------------------------------------------- %
%             Add cyclic prefix                 %
% --------------------------------------------- %
% input: SymIFFT(NumSubc,NumLoop); output: SymCP(NumSubc + NumCP,NumLoop)
NumAddPrefix = NumSubc + NumCP;
SymCP = zeros(NumAddPrefix,NumLoop);
RowPrefix = (NumSubc - NumCP + 1):NumSubc;
SymCP = [SymIFFT(RowPrefix,:);SymIFFT];
% --------------------------------------------- %
%             Go through the channel            %
% --------------------------------------------- %
% input: SymCP(NumSubc + NumCP,NumLoop); output: SymCh(1,(NumSubc + NumCP)*NumLoop)
SymCh = zeros(1,(NumSubc + NumCP)*NumLoop);
SymChtmp = SymCP(:).';   % 進(jìn)行這個(gè)轉(zhuǎn)置操作之后就成了一個(gè)矢量
                      % 相當(dāng)于把矩陣的列向量依次排列 改變?yōu)橐粋€(gè)行向量
Ch = [1 1/2 1/4];
SymChtmptmp = filter(Ch,1,SymChtmp);
%--------------------------------------------------------------------
% 函數(shù)說(shuō)明：
% Firlter data with an infinite impulse response (IIR) or finite impulse response
% (FIR) filter
% y = filter(b,a,X) filters the data in vector X with the filter described by
% numerator coefficient vector b and denominator coefficient vector a. If a(1) is
% not equal to 1, filter normalizes the filter coefficients by a(1). If a(1) equals
% 0, filter returns an error.
%--------------------------------------------------------------------
% If X is a matrix, filter operates on the columns of X. If X is a multidimensional
% array, filter operates on the first nonsingleton dimension.
%--------------------------------------------------------------------
% Add the AWGN
BerSnrTable = zeros(20,3);
for snr=0:19;  % = SNR + 10*log10(log2(2));
    BerSnrTable(snr+1,1) = snr;
SymCh = awgn(SymChtmptmp,snr,'measured');
%--------------------------------------------------------------------
% 函數(shù)說(shuō)明：
%  AWGN Add white Gaussian noise to a signal.
%     Y = AWGN(X,SNR) adds white Gaussian noise to X.  The SNR is in dB.
%     The power of X is assumed to be 0 dBW.  If X is complex, then 
%     AWGN adds complex noise.
%  ------------------------------------------------------------------
%     Y = AWGN(X,SNR,SIGPOWER) when SIGPOWER is numeric, it represents 
%     the signal power in dBW. When SIGPOWER is 'measured', AWGN measures
%     the signal power before adding noise.
% ---------------------------------------------------------------------
%     Y = AWGN(X,SNR,SIGPOWER,STATE) resets the state of RANDN to STATE.
%  
%     Y = AWGN(..., POWERTYPE) specifies the units of SNR and SIGPOWER.
%     POWERTYPE can be 'db' or 'linear'.  If POWERTYPE is 'db', then SNR
%     is measured in dB and SIGPOWER is measured in dBW.  If POWERTYPE is
%     'linear', then SNR is measured as a ratio and SIGPOWER is measured
%     in Watts.
%  
%     Example: To specify the power of X to be 0 dBW and add noise to produce
%              an SNR of 10dB, use:
%              X = sqrt(2)*sin(0:pi/8:6*pi);
%              Y = AWGN(X,10,0);
%  
%     Example: To specify the power of X to be 0 dBW, set RANDN to the 1234th
%              state and add noise to produce an SNR of 10dB, use:
%              X = sqrt(2)*sin(0:pi/8:6*pi);
%              Y = AWGN(X,10,0,1234);
%  
%     Example: To specify the power of X to be 3 Watts and add noise to
%              produce a linear SNR of 4, use:
%              X = sqrt(2)*sin(0:pi/8:6*pi);
%              Y = AWGN(X,4,3,'linear');
%  
%     Example: To cause AWGN to measure the power of X, set RANDN to the 
%              1234th state and add noise to produce a linear SNR of 4, use:
%              X = sqrt(2)*sin(0:pi/8:6*pi);
%              Y = AWGN(X,4,'measured',1234,'linear');
% --------------------------------------------- %
%            Remove Guard Intervals             %
% --------------------------------------------- %
% input: SymCh(1,(NumSubc + NumCP)*NumLoop); output: SymDeCP(NumSubc,NumLoop)
SymDeCP = zeros(NumSubc,NumLoop);
SymDeCPtmp = reshape(SymCh,NumSubc + NumCP,NumLoop);
SymDeCP = SymDeCPtmp((NumCP+1+SyncDelay):NumAddPrefix+SyncDelay,:);
% --------------------------------------------- %
%                     FFT                       %
% --------------------------------------------- %
% input: SymDeCP(NumSubc,NumLoop); output: SymFFT(NumSubc,NumLoop)
SymFFT = fft(SymDeCP,NumSubc,1);
% --------------------------------------------- %
%        Make Decision(Include DeQAM)           %
% --------------------------------------------- %
% SymFFT(NumSubc,NumLoop); output: SymDec(NumSubc,NumLoop)
SymDec = zeros(NumSubc,NumLoop);
SymEqtmp(1,:) = SymFFT(1,:)+i*SymFFT(NumSubc/2+1,:);
SymEqtmp(2:NumSubc/2,:) = SymFFT(2:NumSubc/2,:);
for m = 1:NumLoop
    for n = 1:NumSubc/2
        Real = real(SymEqtmp(n,m));
        Imag = imag(SymEqtmp(n,m));
        
        if( abs((Real -1)) < abs((Real +1)))
            SymDec(2*n-1,m) = 1;
        else
            SymDec(2*n-1,m) = 0;
        end
        if( abs((Imag -1)) < abs((Imag +1  )) )  
            SymDec(2*n,m) = 1;
        else
            SymDec(2*n,m) = 0;
        end
    end
end
% ---------------------------------------------------------------------
% Test by lavabin
% Another way to DeQAM
%  QAMTable = [-1-i -1+i 1-i 1+i];
%  00->-1-i,01->-1+i,10->1-i,11->1+i
TestSymDec = zeros(NumSubc,NumLoop);
TestSymEqtmp(1,:) = SymFFT(1,:)+i*SymFFT(NumSubc/2+1,:);
TestSymEqtmp(2:NumSubc/2,:) = SymFFT(2:NumSubc/2,:);
TestSymEqtmp1 = reshape(TestSymEqtmp,1,NumSubc*NumLoop/2);
min_d = zeros(size(TestSymEqtmp1));
min_ddd = zeros(1,NumSubc*NumLoop);
d = zeros(4,1);
min_index = 0;
for ii = 1:1:(NumSubc*NumLoop/2)
    for jj = 1:4
        d(jj) = abs(TestSymEqtmp(ii) - QAMTable(jj));
    end
     [min_d(ii),min_index] = min(d);
% % [Y,I] = MIN(X) returns the indices of the minimum values in vector I.    
      switch min_index
  case 1
     min_ddd(2*ii-1) = 0 ;
         min_ddd(2*ii) =  0 ;
  case 2
     min_ddd(2*ii-1) = 0 ;
         min_ddd(2*ii) =  1 ;
  case 3
     min_ddd(2*ii-1) = 1 ;
         min_ddd(2*ii) =  0 ;
  case 4
     min_ddd(2*ii-1) = 1 ;
         min_ddd(2*ii) =  1 ;
          otherwise
      fprintf('Impossible error!!! \n\n');
      end
end
%--------------------------------------------------------------------
% 函數(shù)說(shuō)明：
% % C = min(A) returns the smallest elements along different dimensions of an
% % array.
% % If A is a vector, min(A) returns the smallest element in A.
% % If A is a matrix, min(A) treats the columns of A as vectors, returning a row
% % vector containing the minimum element from each column.
% % [C,I] = min(...) finds the indices of the minimum values of A, and returns
% % them in output vector I. If there are several identical minimum values, the
% % index of the first one found is returned.
% Bit Error
BitsRx = zeros(1,NumSubc*NumLoop);
BitsRx = SymDec(:).';
[Num,Ber] = symerr(BitsTx,BitsRx)
BerSnrTable(snr+1,2) = Num ;
BerSnrTable(snr+1,3) = Ber ;
end
%--------------------------------------------------------------------
% Test by lavabin
if min_ddd == BitsRx 
    fprintf('DeQAM two ways the same results \n\n');
else
     fprintf('DeQAM two ways the different results');
end
%--------------------------------------------------------------------
figure(1);
subplot(2,1,1);
semilogy(BerSnrTable(:,1),BerSnrTable(:,2),'o-');
subplot(2,1,2);
semilogy(BerSnrTable(:,1),BerSnrTable(:,3),'o-');
%--------------------------------------------------------------------
time_of_sim = toc
echo on;
% --------------------------------------------- %
%                   The END                     %
% --------------------------------------------- %
%    本程序中得到的收端OFDM信號(hào)的頻譜波形，是與其發(fā)端信號(hào)的排步有關(guān)的。在發(fā)端的
% 載波安排上，128個(gè)載波有前后各32個(gè)載波是null載波(如果這前后各32各載波是帶外頻段，
% 那么理論上它們都應(yīng)該是零！)，中間的64個(gè)載波是數(shù)據(jù)載波。這樣的排步很明顯就是一個(gè)
% 兩邊低，中間高的頻譜形式。所以，收端也應(yīng)該是這個(gè)輪廓。
clc;clear all;close all;echo off;tic;
% -------------------------------------------------------------------
%                      Parameter Definition
% --------------------------------------------------------------
Fd    = 1;             % symbol rate (1Hz)
Fs    = 1*Fd;          % number of sample per symbol
M     = 4;             % kind(range) of symbol (0,1,2,3)
Ndata = 1024;          % all transmitted data symbol 
Sdata = 64;            % 64 data symbol per frame to ifft
Slen  = 128;           % 128 length symbol for IFFT 
Nsym  = Ndata/Sdata;   % number of frames -> Nsym frame
GIlen = 144;           % symbol with GI insertion  GIlen = Slen + GI
GI    = 16;            % guard interval length
% ----------------------------------------------------------------
%                        Vector Initialization
% ----------------------------------------------------------------
X  = zeros(Ndata,1);
Y1 = zeros(Ndata,1);
Y2 = zeros(Ndata,1);
Y3 = zeros(Slen,1);
z0 = zeros(Slen,1);
z1 = zeros(Ndata/Sdata*Slen,1);
g  = zeros(GIlen,1);
z2 = zeros(GIlen*Nsym,1);
z3 = zeros(GIlen*Nsym,1);
% random integer generation by M kinds 
X = randint(Ndata, 1, M);
% digital symbol mapped as analog symbol
% Y1 is a Ndata-by-2 matrix, is changed into Y2 by "amodce"
Y1 = modmap(X, Fd, Fs, 'qask', M);
% covert to complex number
Y2 = amodce(Y1,1,'qam');
% figure(1);
% scatterplot(Y2,length(Y2),0,'bo');grid on;
scatterplot(Y2,Fd,0,'bo');grid on;
title('4-QAM Constellation');
Tx_spectrum = zeros(size(Y3));
for j=1:Nsym;
    for i=1:Sdata;
        Y3(i+Slen/2-Sdata/2,1)=Y2(i+(j-1)*Sdata,1);
        Tx_spectrum = Tx_spectrum + abs(Y3);
    end
    z0=ifft(Y3);
    
    for i=1:Slen;    % generate time-domain vector, z1, without GI
        z1(((j-1)*Slen)+i)=z0(i,1);
    end
    %
    for i=1:Slen;
        g(i+16)=z0(i,1);
    end
    
    for i=1:GI;
        g(i)=z0(i+Slen-GI,1);    
    end
    for i=1:GIlen;  % generate time-domain vector, z2, with GI
        z2(((j-1)*GIlen)+i)=g(i,1);
    end
    
end
% graph on time domain
figure(2);
f = linspace(-Sdata,Sdata,length(z1));
plot(f,abs(z1));
   
Y4 = fft(z1);  % FFT operation, at receiver
% if Y4 is under 0.01 Y4=0.001
for j=1:Ndata/Sdata*Slen;
if abs(Y4(j)) < 0.01
    Y4(j)=0.01;
end
end
Y4 = 10*log10(abs(Y4));  %  Y4 is used for spectrum display.
% graph on frequency domain
figure(3);
f = linspace(-Sdata,Sdata,length(Y4));
plot(f,Y4);grid on;
axis([-Slen/2 Slen/2 -20 20]);
title('Received OFDM signal spectrum');
figure(4);
f = linspace(-Sdata,Sdata,length(Y3));
plot(f,abs(Y3),'b-','LineWidth',6);grid on;
axis([-Slen/2 Slen/2 -2 2]);
title('Transmitted OFDM signal spectrum');
simulation_time = toc