4x4 alamouti matlab code MIMO system

i need a matlab code for 4x4 alamouti scheme...so if any one has any idea regarding to it than please share it. thanks. waiting for ur help.

This a code 2*2 Alamounti used a 16QAM modulation:
clc
clear all;
close all;
M=16;
N = 10^6; % number of bits or symbols
snr=linspace(0,20,21);
h11 = 1/sqrt(2)*(randn(1,N/2) + 1i*randn(1,N/2));
h12 = 1/sqrt(2)*(randn(1,N/2) + 1i*randn(1,N/2));
h21 = 1/sqrt(2)*(randn(1,N/2) + 1i*randn(1,N/2));
h22 = 1/sqrt(2)*(randn(1,N/2) + 1i*randn(1,N/2));
h=[h11;h12;h21;h22];
h = reshape(h,2,N);
HEq = zeros(2,2*N);
HEq,(1:4:end)) = h,(1:2:end));
HEq,(2:4:end)) = h,(2:2:end));
HEq,(3:4:end)) = flipud(h,(1:2:end)));
HEq,(4:4:end)) = flipud(h,(2:2:end)));
HEq(2,=conj(HEq(2,);
HEq(2,(3:4:end))=-1*HEq(2,(3:4:end));
HEq(2,(4:4:end))=-1*HEq(2,(4:4:end));
HEq = reshape(HEq,4,N);
hEq = zeros(2,N);
hEq,(1:2:end)) = reshape(h,(1:2:end)),2,N/2); % [h1 0 ... ; h2 0...]
hEq,(2:2:end)) = kron(ones(1,N/2),[1;-1]).*flipud(reshape(h,(1:2:end)),2,N/2)); % [h1 h2 ... ; h2 -h1 ...]
hEq(1, = conj(hEq(1,); % [h1* h2* ... ; h2 -h1 .... ]
hEqPower3 = sum(hEq.*conj(hEq),1);% mod(h1)^2+mod(h2)^2
hEq1 = zeros(2,N);
hEq1,(1:2:end)) = reshape(h,(2:2:end)),2,N/2); % [h1 0 ... ; h2 0...]
hEq1,(2:2:end)) = kron(ones(1,N/2),[1;-1]).*flipud(reshape(h,(2:2:end)),2,N/2)); % [h1 h2 ... ; h2 -h1 ...]
hEq1(1, = conj(hEq1(1,); % [h3* h4* ... ; h4 -h3 .... ]
hEqPower1 = sum(hEq1.*conj(hEq1),1);% mod(h3)^2+mod(h4)^2
hEqPower3=hEqPower3+hEqPower1;
hEqPower=zeros(1,N/2);
hEqPower=hEqPower3,(1:2:end));
% Transmitter
x = randint(N,1,M);
x5=reshape(x,1,N);
y3=modulate(modem.qammod(M),x);
y3=reshape(y3,1,N);
y5 = zeros(2,N);
y5,1:2:end) = reshape(y3,2,N/2); % [x1 x2; ...]
y5,2:2:end) =(kron(ones(1,N/2),[-1;1]).*flipud(reshape(conj(y3),2,N/2)));% [-x2* x1*; ....]
h1=h,(1:2:end));
h2=h,(2:2:end));
H1 = kron(reshape(h1,2,N/2),ones(1,2));
H2 = kron(reshape(h2,2,N/2),ones(1,2));
ser=zeros(1,length(snr));
tg1=sum(H1.*y5,1);
tg2=sum(H2.*y5,1);
hEqPower6=[hEqPower;hEqPower];
HEq1=conj(HEq);
for i = 1:length(snr)
N11=1/sqrt(2)*(randn(1,N/2)+1i*randn(1,N/2));
N12=1/sqrt(2)*(randn(1,N/2)+1i*randn(1,N/2));
N1=[N11;N12];
N1=reshape(N1,1,N);
y1=tg1+10^(-(snr(i)-10*log10(20))/20)*N1;
N11=1/sqrt(2)*(randn(1,N/2)+1i*randn(1,N/2));
N12=1/sqrt(2)*(randn(1,N/2)+1i*randn(1,N/2));
N1=[N11;N12];
N1=reshape(N1,1,N);
y2=tg2+10^(-(snr(i)-10*log10(20))/20)*N1;
y1=reshape(y1,2,N/2);
y2=reshape(y2,2,N/2);
y=[y1;y2];
y((2:2:end),=conj(y((2:2:end),);
y9=zeros(4,2*N);
y9=[y;y];
y9=reshape(y9,4,N);
some=sum(y9.*HEq1,1);
some=reshape(some,2,N/2)./sqrt(hEqPower6);
dd3=[];
for g=1:2
some1=some(g,;
j0=abs(some1-sqrt(hEqPower).*modulate(modem.qammod(M),0));
j1=abs(some1-sqrt(hEqPower).*modulate(modem.qammod(M),1));
j2=abs(some1-sqrt(hEqPower).*modulate(modem.qammod(M),2));
j3=abs(some1-sqrt(hEqPower).*modulate(modem.qammod(M),3));
j4=abs(some1-sqrt(hEqPower).*modulate(modem.qammod(M),4));
j5=abs(some1-sqrt(hEqPower).*modulate(modem.qammod(M),5));
j6=abs(some1-sqrt(hEqPower).*modulate(modem.qammod(M),6));
j7=abs(some1-sqrt(hEqPower).*modulate(modem.qammod(M),7));
j8=abs(some1-sqrt(hEqPower).*modulate(modem.qammod(M),8));
j9=abs(some1-sqrt(hEqPower).*modulate(modem.qammod(M),9));
j10=abs(some1-sqrt(hEqPower).*modulate(modem.qammod(M),10));
j11=abs(some1-sqrt(hEqPower).*modulate(modem.qammod(M),11));
j12=abs(some1-sqrt(hEqPower).*modulate(modem.qammod(M),12));
j13=abs(some1-sqrt(hEqPower).*modulate(modem.qammod(M),13));
j14=abs(some1-sqrt(hEqPower).*modulate(modem.qammod(M),14));
j15=abs(some1-sqrt(hEqPower).*modulate(modem.qammod(M),15));
rVec = [j0;j1;j2;j3;j4;j5;j6;j7;j8;j9;j10;j11;j12;j13;j14;j15];
[jj, dd] = min(rVec,[],1);
dd3=[dd3;dd];
dd=[];
end
dd3=dd3-1;
dd3=reshape(dd3,1,N);
[num ty]=symerr(x5,dd3);
ser(i)=ty;
dd3=[];
y=[];
end
semilogy(snr,ser,'g-+','linewidth',1);
grid on;hold on;
title('Alamouti nTx=2, nRx=2, for 16-QAM','Color','k','FontSize',13);
legend('sim (nTx=2, nRx=2, Alamouti(16-QAM))','location','southwest');
xlabel('SNR(dB) ---->','Color','k','FontSize',11);Ylabel('Symbol Error Rate ---->','Color','k','FontSize',11);