2x2 MIMO-OFDM with ZF EQ. Matlab Code Help

hello,
I'm working in my Graduation Project on MIMO-OFDM Implementation using Matlab
i made the following code but i m a little lost , can't figure out what to do next ,
i generated the data , modulated them , then branched the data into 2 branches to be sent , but i cant figure out how to multiply the channel matrix by the modulated data matrix .....
check the code :

% 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

clear 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 = 100; % 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



% Transmitter
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.')).';
x2t = (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/nTx) + j*randn(nRx,nSym/nTx)]; % white gaussian noise, 0dB variance
H = 1/sqrt(2)*[randn(nRx,nTx,nSym/nTx) + j*randn(nRx,nTx,nSym/nTx)]; % Rayleigh fading channel


any help here please ?

Code looks fine for the most part, the comments make things pretty easy to follow =)

One thing, putting a 0 at DC is pretty non-standard as far as I know, typically subcarriers are mapped from (in this case) -26 to 25. This shouldn't be a big deal as long as you deal with it at the transmitter.

The way you are generating noise is fine. However, the channel you are generating is a flat Rayleigh fading model; there is no multipath involved. Therefore, the channel coefficient will be the same across all frequencies in your case. If you want to model it in this way, then that is fine, but one of the main advantages of OFDM is robustness to multipath.

To model this, you would need to generate a channel response (that is less than the guard interval) for each antenna pair (4 channels total). Then, when you take the FFT, the channel you see at the output will be flat, but will vary across the subcarriers.

If you want to keep the model as-is, then you can do an element-wise multiplication of the square root of the norm of each element of H with your channel. If we define H[k] as your channel matrix for symbol k, then the channel gain between receiver i and transmitter j is given by

\[H_{ij}[k] = \sqrt{\frac{x^2+y^2}{2}}, x,y\sim\, N(0,1)\]

Then, multiply the channel with the transmitted signal and add noise:

\[\left[\begin{array}{c}y_1[n] \\ y_2[n]\end{array}\right] = \left[\begin{array}{cc}h_{11}[k] & h_{12}[k] \\ h_{21}[k] & h_{22}[k]\end{array}\right]\left[\begin{array}{c}x_1[n] \\ x_2[n]\end{array}\right]+\left[\begin{array}{c}w_1[n] \\ w_2[n]\end{array}\right],n=1,\ldots,80\]

You can vary the channel from symbol to symbol, but it should stay constant across the symbol. Hope that helps you get your received signal. Once you receive the signal you can do things like zero-forcing beamforming (which is easier in this flat model than in the multipath model, where you would have to do it on a subcarrier-by-subcarrier basis).

i made some modifications to the code and it runs perfectly except that the simulated bit error rate curve is wrong and i cant figure out what is the problem .....

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Author : Abd-El-Mageed Samy
% Email : amageed.samy@yahoo.com
% 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

clear 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

for ii = 1:length(EbN0dB)

% ====================>> T R A N S M I T T E R <<====================

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.')).';
x2t = (nFFT/sqrt(nDSC))*ifft(fftshift(x2F.')).';

% 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

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

W = PseudoInv(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 = (sqrt(nDSC)/nFFT)*fftshift(fft(y1t.')).';
y2F = (sqrt(nDSC)/nFFT)*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, ];

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

ipModHat = 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
nErr(ii) = size(find(ipBit - ipBitHat),2);
end

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


close all; figure
semilogy(EbN0dB,simBer,'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')

hi
Error in program
W = PseudoInv(H); % Pseudo Inverso of the Channel Matrix (Equalization Matrix)
can you explain this error please

- - - Updated - - -

hi
Error

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

please can you explain
thank's