www.pudn.com > space-timecodingBrankaVuceticJinhongYuan.rar > vblast_est.m, change:2004-01-07,size:6749b
function z=vblast_est(alg,modulation,corr,alpha,system)
%**************************************************************************
% This program implements 1x2 and 2x2 Vblast based on preamble type of
% channel estimation using MMSE estimation, on a frame of 130 symbols.
% z=vblast_est(alg,modulation,corr,alpha,system) where
% alg -> algorithm for vblast receiver.Choose from 'ZF' (for zero-forcing),
% 'MM'(MMSE estimation)
% modulation -> 'BPSK,'QPSK','16QAM','64QAM'
% corr -> 1 for correlation at the receiver(2x2 correlation only implemented),
% and 0 for no correlation
% alpha -> correlation coefficient value from 0 -> 1 and 0 in case corr = 0;
% system -> 12 for a 1x2 system and 22 for a 2x2 system
% EXAMPLE: vblast_est('ZF','QPSK',1,0.5,22)plots ZF curve for QPSK with 0.5
% correlation coefficient at the receiver and channel estimation for a 2x2
% system
%**************************************************************************
% identify bits for each type of modulation
switch modulation
case 'BPSK'
BITS=1;
case 'QPSK'
BITS=2;
case '16QAM'
BITS=4;
case '64QAM'
BITS=6;
end
%define frame
K=130;
% define SNR range
EbNo=[0:2:30];
%define plotting axis
SNR_axis=[];
BER_axis=[];
%set initial count
idx=1;
%define numbers of antennas
if system==12
M=2;%rx antennas
N=1;%tx antennas
elseif system==22
M=2;%rx antennas
N=2;%tx antennas
end
%number of monte-carlo runs
Num=10;
%parameters for wait bar
h = waitbar(0,'Please wait...');
wb=6.25;
%clear BER register
BER=[];
%commence SNR loop
for SNR=EbNo
errors=0;
%define standard deviation, sigma, for noise
sigma=0.5/(sqrt(N)*10^(SNR/10));
for iter=1:Num %commence iteration loop
%define a random Rayleigh channel
H=(randn(M,N)+j*randn(M,N))/sqrt(2);
%exercise correlation option, if any
if corr & system==22
R=chol([1 alpha;alpha 1]);%cholesky factor of correlation matrix, i.e. 'R' square root of correlation matrix
H=R*H; % R^0.5 *H -> correlation at receiver
end
H_save=H;%assign H to unalterable value
% modulated input data
tx_bits=randn(K,N,BITS)>0;
temp1=[];
for i=1:K
d1=tx_modulate(tx_bits(i,:),modulation);
temp1=[temp1; d1];
end
d=temp1;
%insert pilots
if system==12
pilots=[0 0 0 0 0 0 0 0 0 0 1 -1 -1 1 -1 -1 1 -1 1 1];
d=[pilots.'; temp1];
elseif system==22
pilots=[0 0 0 0 0 0 0 0 0 0 1 -1 -1 1 -1 -1 1 -1 1 1;1 -1 -1 -1 1 -1 -1 1 -1 1 0 0 0 0 0 0 0 0 0 0];
d=[pilots.'; temp1];
end
%AWGN noise
AWGN_noise = sqrt(sigma)*(randn(150, M)+j*randn(150, 2));
%receiver signal vector added to AWGN noise
r = ((H_save*d.')/sqrt(N)).' + AWGN_noise;
pilots=pilots.';
tr_symbols=r(1:20,:);
R=r(21:end,:).';
if system==12
h11=mean(tr_symbols(:,1).*conj(pilots(:,1)));
h12=mean(tr_symbols(:,2).*conj(pilots(:,1)));
est_coefs=[h11;h12];
elseif system==22
h11=mean(tr_symbols(:,1).*conj(pilots(:,1)));
h12=mean(tr_symbols(:,1).*conj(pilots(:,2)));
h21=mean(tr_symbols(:,2).*conj(pilots(:,1)));
h22=mean(tr_symbols(:,2).*conj(pilots(:,2)));
est_coefs=[[h11;h21] [h12; h22]];
end
H_save=est_coefs;
for i=1:K
r= R(:,i);
X=temp1.';
X=X(:,i);
H=H_save;
%Zero-forcing algorithm
if alg=='ZF'
%initialization
G=pinv(H);
[gk k0]=min(sum(abs(G).^2,2));
for m=1:N % This FOR loop determines the ordering sequence k1 and determines the 'a' matrix.
%This is just one run,i.e. one for each H matrix.
%The 'a' matrix is automatically sorted as [a1 a2...aM]
k1(m)=k0;
w(m,:)=G(k1(m),:);
y=w(m,:)*r;
a(k1(m),1)=Q(y,modulation);
r = r - a(k1(m)) * H_save(:, k1(m));
H(:,k0)=zeros(M,1);
G=pinv(H);
for t=1:m
G(k1(t),:)=inf;
end
[gk k0]=min(sum(abs(G).^2,2));
end
%MMSE algorithm
elseif alg=='MM'
%initialisation
G=inv(H'*H+N/(10^(0.1*SNR))*eye(N))*H';
[gk k0]=min(sum(abs(G).^2,2));
for m=1:N % This FOR loop determines the ordering sequence k1 and determines the 'a' matrix.
% This is just one run,i.e.one for each H matrix.The 'a' matrix is automatically sorted
% as [a1 a2...aM]
k1(m)=k0;
w(m,:)=G(k1(m),:);
y=w(m,:)*r;
a(k1(m),1)=Q(y,modulation);
r = r - a(k1(m)) * H_save(:, k1(m));
H(:,k0)=zeros(M,1);
G=inv(H'*H+N/(10^(0.1*SNR))*eye(N))*H';
for t=1:m
G(k1(t),:)=inf;
end
[gk k0]=min(sum(abs(G).^2,2));
end
end
%count errors
if BITS==1 %for BPSK modulation
errors(iter) = sum((sign(real(a))~=sign(real(X))));
else %for QPSK,16QAM and 64QAM modulations
errors(iter) = sum((sign(real(a))~=sign(real(X))) | sign(imag(a))~=sign(imag(X)));
end
end %end of iteration loop Loop
end
BER(idx)=sum(errors)/(Num) ; % Calculate BER after completion of 'Num' runs
SER(idx)=BER(idx)*BITS; %calculate symbol error rate
idx=idx + 1; %increment count
waitbar(wb/100);
wb=wb+6.25;%increment wait bar
end %end of SNR loop
close(h);%terminate wait bar
SNR_axis=EbNo;
BER_axis=[BER_axis BER];
SER_axis=SER;
%plot BER
semilogy(SNR_axis,BER_axis,'b-*');
xlabel('SNR [dB]');
ylabel('BER/SER');
title('BER/SER Plots');
hold;
%plot SER
semilogy(SNR_axis,SER_axis,'b-o');
axis([0 30 1e-6 1]);
grid on;
legend('BER','SER');