vblast.m

供初学空时编码（vblast接收）的matlab仿真程序

close all;
clear all;
Tx_n = 4;
Rx_n = 4;
index = 2;%for modu index:1/2/4/6 bits per symbol
frame_length = 48;
numloop=100;
% SNR_dB=5;
SNR_dB_MIN=0;
SNR_dB_MAX=25;
SNR_dB_stmp=5;
bertmp=zeros(1,numloop);
BER=zeros(1,length(SNR_dB_MIN:SNR_dB_stmp:SNR_dB_MAX));
if Rx_n < Tx_n
	disp_str = ['error!!!the Rx_n must no less than Tx_n'];
	disp(disp_str);
	return;
end
for SNR_dB=SNR_dB_MIN:SNR_dB_stmp:SNR_dB_MAX
   for looptmp=1:numloop
source = rand(1,Tx_n*frame_length*index)>.5;
column_num = length(source) / Tx_n;
% encode; the dim of encoder_output is Tx_n*(frame_length*index)
encoder_output = zeros(Tx_n, column_num);
for col = 1:column_num
	for row = 1:Tx_n
		encoder_output(row, col) = source((col-1)*Tx_n + row);
	end
end
% modulation ; the dim of modu_output is Tx_n  fream_length
modu_output=zeros(Tx_n,frame_length);
for t=1:Tx_n
    modu_output(t,:) = modulation(encoder_output(t,:), index);
end
% go to the channel
  % the H of the channel with the dim is  Rx_n Tx_n
  H=zeros(Rx_n,Tx_n);
for rx = 1:Rx_n
	for tx = 1:Tx_n
		randn('state',sum(100*clock));
		x = randn(1); 
		y = randn(1);
		alpha = sqrt(x^2+y^2);    %rayleigh distributed
		theta = 2*pi*rand(1);   
		H(rx,tx) = alpha*exp(j*theta);
	end
end
  % the dim of channel_output is Tx_n , frame_length
channel_output = H*modu_output;
module_signal = [];
		for t=1:Tx_n
			module_signal =[ module_signal abs(channel_output(t,:))];			
		end
		signal_level = mean(abs(module_signal).^2);
		signal_level_dB = 10*log10(signal_level);
  % add the AWGN nosie
decoder_input = awgn(channel_output,SNR_dB,signal_level_dB);
delta2 = 1/(10^(SNR_dB/10));
r = decoder_input;
y = zeros(Tx_n,frame_length);

G=pinv(H);
[gk k0]=min(sum(abs(G).^2,2));
      
 for m=1:Tx_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),:)=Q(y,index);
   r = r -H(:, k1(m))*a(k1(m),:);   
   H(:,k0)=zeros(Rx_n,1);
   G=pinv(H);
  for t=1:m
   G(k1(t),:)=inf;
  end
   [gk k0]=min(sum(abs(G).^2,2))
 end
 
% demodulation
for t = 1:Tx_n
demodu_output(t,:) = demodulation(decoder_output(t,:), index);
end
sink = zeros(1,Tx_n*frame_length*index);
%demultiplex
for i = 1:Tx_n
    sink(i:Tx_n:end)=demodu_output(i,:);
end
source=double(source);
[Num,Ber]=symerr(sink,source);
bertmp(looptmp)=Ber;
end
BER((SNR_dB-SNR_dB_MIN)/SNR_dB_stmp+1)=mean(bertmp);
end
figure;
semilogy(SNR_dB_MIN:SNR_dB_stmp:SNR_dB_MAX,BER,'*:');
grid on;
xlabel('SNR(db)');ylabel('BER');
title('V-BLAST of MMSE(no order)');
legend('Tx = 4','Rx = 4','QPSK modulation');