MiMo-OFDM ( I got this code from tis website but the result is wrong) why?

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%
% Author : Abd-El-Mageed Samy
% Email : [email]amageed.samy@yahoo.com[/email]
% Version : 1.2
% Date : 14th March 2013
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%

% Script for computing the BER of 2x2 MIMO OFDM System
% in Rayleigh Fading Channel with Zero Force Receiver.

% Note : Simulation Parameters are OFDM Specs for the IEEE 802.11a Standard
clc
clear all;
close all;
nFFT = 64; % fft size
nDSC = 52; % number of data subcarriers
nBitPerSym = 52; % number of bits per OFDM symbol (same as the number of subcarriers for BPSK)
nSym = 10^4; % number of symbols
totalBits = 2*nBitPerSym*nSym; % total number of Data Bits
bitsPerBranch = nBitPerSym*nSym; % total number of data Bits in each of the two branches
nTx = 2 ; nRx = 2 ; % 2x2 MIMO System
EbN0dB = [0:10]; % bit to noise ratio
EsN0dB = EbN0dB + 10*log10(nDSC/nFFT) + 10*log10(64/80); % converting Bit Energy into Symbol energy

no_of_error=[];
ratio=[];

for ii = 1:length(EbN0dB)

% ====================>> T R A N S M I T T E R <<====================
cc=rand(1,totalBits);
ipBit = rand(1,totalBits)> 0.5 ; % random 1's and 0's , multiplied by 2 for the two branches of the transmitting Antennas
ipMod = 2*ipBit-1; % BPSK modulation 0 --> -1, 1 --> +1
ipMod1 = ipMod(1,1:bitsPerBranch); % 1st Branch Data


ipMod2 = ipMod(bitsPerBranch+1:totalBits); % 2nd Branch Data
ipMod1 = reshape(ipMod1,nBitPerSym,nSym).'; % 1st Branch grouping into multiple symbols and transposing as Pre-IFFT Stage
ipMod2 = reshape(ipMod2,nBitPerSym,nSym).'; % 2nd Branch grouping into multiple symbols and transposing as Pre-IFFT Stage

% Assigning modulated symbols to subcarriers from [-26 to -1, +1 to +26]
x1F = [zeros(nSym,6) ipMod1(:,[1:nBitPerSym/2]) zeros(nSym,1) ipMod1(:,[nBitPerSym/2+1:nBitPerSym]) zeros(nSym,5)] ;
x2F = [zeros(nSym,6) ipMod2(:,[1:nBitPerSym/2]) zeros(nSym,1) ipMod2(:,[nBitPerSym/2+1:nBitPerSym]) zeros(nSym,5)] ;
% Taking FFT, the term (nFFT/sqrt(nDSC)) is for normalizing the power of transmit symbol to 1 
%x1t = (nFFT/sqrt(nDSC))*ifft(fftshift(x1F.')).';
a=max(max(abs(x1F)));
b=max(max(abs(x2F)));
x1t = ifft(fftshift(x1F.')).';
x1t=x1t./a;
%x2t = (nFFT/sqrt(nDSC))*ifft(fftshift(x2F.')).';
x2t = ifft(fftshift(x2F.')).';
x2t=x2t./b;
c=(nFFT/sqrt(nDSC))*ifft(fftshift(x1F.')).';
% Appending cylic prefix
x1t = [x1t(:,[49:64]) x1t];
x2t = [x2t(:,[49:64]) x2t];

% Generating White Gaussian Noise and the Rayleigh fading Channel parameters
N = 1/sqrt(2)*[randn(nRx,nSym*80) + j*randn(nRx,nSym*80)]; % white gaussian noise, 0dB variance
H = 1/sqrt(2)*[randn(nRx,nTx) + j*randn(nRx,nTx)]; % Rayleigh fading channel

% Adding Noise and Channel
xt = [reshape(x1t,1,nSym*80) ; reshape(x2t,1,nSym*80)]; % Formating & Concatenating the data from both antennas
%yt = sqrt(80/64)*H*xt + 10^(-EsN0dB(ii)/20)*N ; % y = Hx + n // the term sqrt(80/64) is to account for the cyclic prefix
yt = awgn(xt,ii);

% =======================>> R E C E I V E R <<=======================

W = pinv(H); % Pseudo Inverso of the Channel Matrix (Equalization Matrix)

% Pseudo Inverse Matrix is a 2x2 Matrix.
% yt is a 2x8000 Matrix where y1(1x8000) is the first row & y2(2x8000) is the 2nd row.

% Splitting into two Branches and reshaping for FFT
y1t = yt(1,:); 
y2t = yt(2,:);
y1t = reshape(y1t,nFFT+nFFT/4,nSym).'; 
y2t = reshape(y2t,nFFT+nFFT/4,nSym).'; 


% Removing Cyclic Prefix
y1t = y1t(:,[1:nFFT]); % 100 x 64
y2t = y2t(:,[1:nFFT]); % 100 x 64
% Applying FFT and adding the term (sqrt(nDSC)/nFFT) for normalization
y1F = a*fftshift(fft(y1t.')).'; 
y2F = b*fftshift(fft(y2t.')).'; 

% Removing Zeros
y1F = [ y1F(:,7:32) y1F(:,34:59) ];
y2F = [ y2F(:,7:32) y2F(:,34:59) ];

% ---------------------------- EQUALIZATION ----------------------------

% Reshaping to prepare for Equalization
y1F = reshape(y1F.',1,nBitPerSym*nSym);
y2F = reshape(y2F.',1,nBitPerSym*nSym);


Wf = fft(W); % Frequency Response of the Channel's Inverse
yF = [ y1F ; y2F ] ; % Reshaping into a 2x5200 Matrix for Equalization

yHat = Wf*yF ; % Equalization ==> Y = WX

% ---------------------------- DeModulation ----------------------------
ipModHat = real(yHat);
ipModHat = [ ipModHat(1,:) ipModHat(2,:) ];
f=ipModHat;

% +ve value --> 1, -ve value --> -1

ipModHat =  2*floor(ipModHat/2) + 1;
e = 2*floor(ipModHat/2) + 1;
ipModHat(find(ipModHat>1)) = +1;
ipModHat(find(ipModHat<-1)) = -1;

% converting modulated values into bits
ipBitHat = (ipModHat+1)/2;

% counting the errors

[no_of_error(ii),ratio(ii)]=biterr(ipBit , ipBitHat);
end


theoryBer = (1/2)*erfc(sqrt(10.^(EbN0dB/10)));


close all; figure
semilogy(EbN0dB,ratio,'mo-','LineWidth',2);
hold on
semilogy(EbN0dB,theoryBer,'bs-','LineWidth',2);



axis([0 10 10^-5 1])
grid on
legend('simulation' , 'theory');
xlabel('Eb/No, dB')
ylabel('Bit Error Rate')