mimo code
Hi, try this program:
clear all; close all;
N=20; % number of users
M=20; % spreading gain
L=4; % channel order
Z=4; % number of antenna
Q=M-L;
%Symbols=500;% no of data vectors
cx=hadamard(2*N);
cn=cx(2:N+1,1:N);clear cx;
SNR=25; % signal to noise ratio in dB
%%% autocorrelation matrix of the spreading code
Rcc=cn*cn';
po=51;%number of points for the plot
nj=50:5:300;
% po=4;%number of points for the plot
% nj=50:50:200;
%nj=[50 100 150 200 250 300];%size of data vectors
% n=200;%number of generated symbols
ni=po;%number of iterations
%var=logspace(-3,-1,ni);%noise variance
var=0.07;
sig=sqrt(var);
%%% parameters for the adaptive algorithms
beta=0.05;gm=0.4;
% Difference=0;
%====================
%Channel
%====================
chn=[
-0.0490+0.3590i
0.4820-0.5690i
-0.5560+0.5870i
1.0000
-0.1710+0.0610i
0.4430-0.0364i
1.0000
0.9210-0.1940i
0.1890-0.2080i
-0.0870-0.0540i
-0.2110-0.3220i
-0.1990+0.9180i
1.0000
-0.2840-0.5240i
0.1360-0.1900i
0.4170+0.0300i
1.0000
0.8730 + 0.1450i
0.2850 + 0.3090i
-0.0490 + 0.1610i];
% channel r=1
c1=[-0.049+0.359i 0.482-0.569i -0.556+0.587i 1 -0.171+0.061i zeros(1,M-L-1)];
r1=[-0.049+0.359i zeros(1,M-1)];
H1=toeplitz(c1,r1);
% channel r=2
c2=[0.443-0.0364i 1 0.921-0.194i 0.189-0.208i -0.087-0.054i zeros(1,M-L-1)];
r2=[0.443-0.0364i zeros(1,M-1)];
H2=toeplitz(c2,r2);
% channel r=3
c3=[-0.211-0.322i -0.199+0.918i 1 -0.284-0.524i 0.136-0.190i zeros(1,M-L-1)];
r3=[-0.211-0.322i zeros(1,M-1)];
H3=toeplitz(c3,r3);
%channel r=4
c4=[0.417+0.030i 1 0.873+0.145i 0.285+0.309i -0.049+0.161i zeros(1,M-L-1)];
r4=[0.417+0.030i zeros(1,M-1)];
H4=toeplitz(c4,r4);
R=[zeros(M-L,L) eye(M-L)];%to generate truncate matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% repeat the process for a variety of SNR
K=50;%the number of averaging= block size
snr=zeros(1,ni);
mci=50;%no of MC runs
er1=zeros(ni,mci); er2=er1; msx=er1;
mse=zeros(ni,1);nm=mse;
it=1;
for p=1:ni,
for mc=1:mci,%100 Monte Carlo runs
fprintf('%d th MC run\n',it);it=it+1;
n=nj(p);
he=zeros(4*(L+1),n);
szf=zeros(1,n);%detected transmitted symbols
smm=zeros(1,n);%detected transmitted symbols
u=zeros(M,n);
s=randsrc(N,n);%generate tx BPSK symbols
u=cn*s;%size M x n
Ruu=zeros(M);
%generate correlation matrix
for m=1:K,
Ruu=Ruu+[u
,m)*u
,m)'];
end;
Ruu=Ruu/K;
%%%%%%%%%%%%%%%%%%%%%%%%%% 1st antenna %%%%%%%%%%%%%%%%%%%%%%%%
n1=sig*(1/sqrt(2))*[randn(Q,n)+sqrt(-1)*randn(Q,n)];%generate random noise
x1=(R*H1*cn*s)+n1;%received vector for the 1st antenna, size of Q x n
%%%%%%%%%%%%%%%%%%%%%%%%%% 2nd antenna %%%%%%%%%%%%%%%%%%%%%%%%
n2=sig*(1/sqrt(2))*[randn(Q,n)+sqrt(-1)*randn(Q,n)];%generate random noise
x2=(R*H2*cn*s)+n2;%received vector for the 2nd antenna, size of Q x n
%%%%%%%%%%%%%%%%%%%%%%%%%% 3rd antenna %%%%%%%%%%%%%%%%%%%%%%%%
n3=sig*(1/sqrt(2))*[randn(Q,n)+sqrt(-1)*randn(Q,n)];%generate random noise
x3=(R*H3*cn*s)+n3;%received vector for the 3rd antenna, size of Q x n
%%%%%%%%%%%%%%%%%%%%%%%%%% 4th antenna %%%%%%%%%%%%%%%%%%%%%%%%
n4=sig*(1/sqrt(2))*[randn(Q,n)+sqrt(-1)*randn(Q,n)];%generate random noise
x4=(R*H4*cn*s)+n4;%received vector for the 4th antenna, size of Q x n
%%%%%%%%%%%%%%
x=[x1.' x2.' x3.' x4.'].';%stack all received vectors from the 4 antennas, size 4Q x n
Nos=sum(sum((abs(n1).^2),1))/(Q*n);%average noise power
sn=sum(sum((abs(x1-n1).^2),1))/(n*N);%energy per-bit, BPSK= 1 bit per-symbol
snr(p)=10*log10(sn/Nos);
%%%%% Blind Subspace Channel Estimation %%%%%%%%%%
Rxx=zeros(4*Q,4*Q);
%generate autocorrelation matrix
for m=1:K,
Rxx=Rxx+[x
,m)*x
,m)'];
end;
Rxx=Rxx/K;
[V,D]=eig(Rxx);%eigen value decomposition, size of V and D is 4Q x 4Q
Gn=V
,1
4*Q)-M);%noise subspace
Gj1=zeros(L+1,M*((4*Q)-M));Gj2=Gj1;Gj3=Gj1;Gj4=Gj1;
Gj=zeros(Z*(L+1),M*((4*Q)-M));
%%%%%%%%
c=0;
for k=1
4*Q)-M,
b=0;
%%% noise subspace from the first antenna
for i=1:L+1,
a=Q;
for j=1:Q, Gj1(i,j+b+c)=Gn(a,k); a=a-1; end;
b=b+1;
end;
b=0;
%%% noise subspace from the 2nd antenna
for i=1:L+1,
a=2*Q;
for j=1:Q, Gj2(i,j+b+c)=Gn(a,k); a=a-1; end;
b=b+1;
end;
b=0;
%%% noise subspace from the third antenna
for i=1:L+1,
a=3*Q;
for j=1:Q, Gj3(i,j+b+c)=Gn(a,k); a=a-1; end;
b=b+1;
end;
b=0;
%%% noise subspace from the fourth antenna
for i=1:L+1,
a=4*Q;
for j=1:Q, Gj4(i,j+b+c)=Gn(a,k); a=a-1; end;
b=b+1;
end;
c=c+M;
end;
Gj=[Gj1.' Gj2.' Gj3.' Gj4.'].';%size Z*(L+1) x M*((4*Q)-M)
Qx=zeros(Z*(L+1));clear a;a=1:M;
for i=1
4*Q)-M, Qx=Qx+(Gj
,a)*cn*cn'*Gj
,a)');a=a+M;end;
%%%%%% Obtaining channel estimate, hest
[hest,A]=eigs(Qx,1,'SM');%return eigenvector corresponds to the minimum eigenvalue of Q, size 4(L+1) x 1
hest=hest*(c1(1)/hest(1));%phase ambiguity removal
%%%%%%%%%% Blind Adaptive Subspace Tracking Algorithm %%%%%%%%%%%
%compute the channel estimate for the whole tx symbols
Hq1=zeros(Q,M);Hq2=Hq1;Hq3=Hq1;Hq4=Hq1;HQ=zeros(4*Q,M);
for m=1:n,
%NOOJA Algorithm
%noise subspace, algorithm at i interval
fn=Gn'*x
,m);
gn=Gn*fn;
pn=x
,m)-gn;
bop=((beta)/[norm(x
,m))-norm(fn)+gm]);
phi=1/sqrt(1+(bop^2*norm(pn)*norm(fn)));
tau=(phi-1)/norm(fn);
pi=((-tau*gn)/bop)+(phi*pn);
chi=pi/norm(pi,1);
ka=Gn'*chi;
Gnn=Gn-(2*chi*ka');%noise subspace for the next symbol
%%%%%%%% construct matrix Q based on the estimate Gn
c=M;
for k=1
4*Q)-M,
b=0;
%%% noise subspace from the first antenna
for i=1:L+1,
a=Q;
for j=1:Q, Gj1(i,j+b+c)=Gnn(a,k); a=a-1; end;
b=b+1;
end;
b=0;
%%% noise subspace from the 2nd antenna
for i=1:L+1,
a=2*Q;
for j=1:Q, Gj2(i,j+b+c)=Gnn(a,k); a=a-1; end;
b=b+1;
end;
b=0;
%%% noise subspace from the third antenna
for i=1:L+1,
a=3*Q;
for j=1:Q, Gj3(i,j+b+c)=Gnn(a,k); a=a-1; end;
b=b+1;
end;
b=0;
%%% noise subspace from the fourth antenna
for i=1:L+1,
a=4*Q;
for j=1:Q, Gj4(i,j+b+c)=Gnn(a,k); a=a-1; end;
b=b+1;
end;
c=c+M;
end;
Gj=[Gj1.' Gj2.' Gj3.' Gj4.'].';%size Z*(L+1) x M*((4*Q)-M)
Qx=zeros(Z*(L+1));clear a;a=1:M;
for i=1
4*Q)-M, Qx=Qx+(Gj
,a)*cn*cn'*Gj
,a)');a=a+M;end;
%%%% compute the channel estimate for the next symbol
mu=0.0001;
ben=hest'*Qx*hest;
hen=hest-(mu*[(norm(hest)*Qx)-(ben*eye(4*(L+1)))]*hest);%the updated channel estimate
%%%% store all channel estimates into he
if m==1, he
,m)=hest; else he
,m)=hen; end;
%%%%%% Construct the Hq matrix for each antenna %%%%%%%
b=0;
%%% the first antenna
for i=1:Q,
a=L+1;
for j=1:L+1, Hq1(i,j+b)=he(a,m); a=a-1; end;
b=b+1;
end;
b=0;
%%% the 2nd antenna
for i=1:Q,
a=2*(L+1);
for j=1:L+1, Hq2(i,j+b)=he(a,m); a=a-1; end;
b=b+1;
end;
b=0;
%%% the third antenna
for i=1:Q,
a=3*(L+1);
for j=1:L+1, Hq3(i,j+b)=he(a,m); a=a-1; end;
b=b+1;
end;
b=0;
%%% the fourth antenna
for i=1:Q,
a=4*(L+1);
for j=1:L+1, Hq4(i,j+b)=he(a,m); a=a-1; end;
b=b+1;
end;
HQ=[Hq1.' Hq2.' Hq3.' Hq4.'].';%the channel estimates, size 4Q x M
%%%%%%%%%% Equalisation using ZF and MMSE equalisers
Gzf=inv(HQ'*HQ)*HQ';%ZF, size M x 4Q
% clear tmp;tmp=Gzf;Gzf=tmp.';%change the matrix orientation
% Rxx=zeros(4*Q,4*Q);
% %generate autocorrelation matrix
% for m=1:K,
% Rxx=Rxx+[x
,m)*x
,m)'];
% end;
% Rxx=Rxx/K;
% Gmmse=Rux*inv(Rxx);%mmse
Gmmse=Ruu*HQ'*pinv(Rxx);%mmse, size M x 4Q
%%%Detection for the 1st user
%%%using ZF detector/equaliser
sx
,m)=sign(real(cn'*Gzf*x
,m)));%size M x n
%%%using MMSE detector/equaliser
sy
,m)=sign(real(cn'*Gmmse*x
,m)));%size M x n
end;%end of loop for m
%%%% obtain the estimated transmitted data for user 1
szf=sx(1,
;
smm=sy(1,
;
%%% to calculate the BER
%ZF
for k=1:n,
if szf(k)==s(1,k), b=0;
else b=1;
end;
er1(p,mc)=er1(p,mc)+b;
end;
er1(p,mc)=er1(p,mc)/n;
%MMSE
for k=1:n,
if smm(k)==s(1,k), b=0;
else b=1;
end;
er2(p,mc)=er2(p,mc)+b;
end;
er2(p,mc)=er2(p,mc)/n;
%%% compute NRMSE
clear tmp; tmp=0;
for i=1:n, tmp=tmp+norm(he
,i)-chn); end;
msx(p,mc)=tmp/n;%MSE
end;%end of 'mc' loop
mse(p)=mean(msx(p,
);
nm(p)=(1/norm(chn))*sqrt(mse(p));%NRMSE
end;%end of 'p' loop
% figure;
% semilogy(snr,er1,'-*',snr,er2,'-o');
% xlabel('SNR (dB)');
% ylabel('BER');
% legend('ZF','MMSE');
% title('Performance of the Detectors for the DS-CDMA');
%
% figure;
% semilogy(snr,mse,'-*');
% xlabel('SNR (dB)');
% ylabel('MSE');
% legend('DP');
% title('Performance of the Channel Estimator for the DS-CDMA');
save dat5 nj nm;
figure;
plot(nj,nm,'-o');
xlabel('Number of iterations');
ylabel('NRMSE');
legend('DP');
title('Performance of the Channel Estimator for the DS-CDMA');