%%STBC-MIMO-OFDM.m
clear all;
i=sqrt(-1);
IFFT_bin_length=512; %傅立叶变换抽样点数目
carrier_count=100; %子载波数目
symbols_per_carrier=66; %符号数/载波
cp_length=10; %循环前缀长度
addprefix_length=IFFT_bin_length+cp_length;
M_psk=4;
bits_per_symbol=log2(M_psk); %位数/符号
%O=[1 2;-2+j 1+j]; %Alamouti Scheme [x1 x2;-x2* x1*] 二天线发送矩阵
O=[1 -2 -3;2+j 1+j 0;3+j 0 1+j;0 -3+j 2+j]; %for [x1 -x2 -x3;x2* x1* 0;x3* 0 x1*;0 -x3* x2*] 三天线发送矩阵
co_time=size(O,1);
Nt=size(O,2); %发射天线数目
Nr=2; %接收天线数目
num_X=1;
for cc_ro=1:co_time
for cc_co=1:Nt
num_X=max(num_X,abs(real(O(cc_ro,cc_co))));
end
end
co_x=zeros(num_X,1);
for con_ro=1:co_time %计算delta,epsilon,eta and conj matrices
for con_co=1:Nt %用于确定矩阵“O”中元素的位置，符号以及共轭情况
if abs(real(O(con_ro,con_co)))~=0
delta(con_ro,abs(real(O(con_ro,con_co))))=sign(real(O(con_ro,con_co)));
epsilon(con_ro,abs(real(O(con_ro,con_co))))=con_co;
co_x(abs(real(O(con_ro,con_co))),1)=co_x(abs(real(O(con_ro,con_co))),1)+1;
eta(abs(real(O(con_ro,con_co))),co_x(abs(real(O(con_ro,con_co))),1))=con_ro;
coj_mt(con_ro,abs(real(O(con_ro,con_co))))=imag(O(con_ro,con_co));
end
end
end
eta=eta.';
eta=sort(eta);
eta=eta.';
carriers = (1: carrier_count) + (floor(IFFT_bin_length/4) - floor(carrier_count/2)); % 坐标： (1 to 100) + 14=(15:114)
conjugate_carriers=IFFT_bin_length-carriers+2; % 坐标 ：256 - (15:114) + 1= 257 - (15:114) = (242:143)
tx_training_symbols=training_symbol(Nt,carrier_count);
baseband_out_length = carrier_count * symbols_per_carrier;
snr_min=3; %最小信噪比
snr_max=15; %最大信噪比
graph_inf_bit=zeros(snr_max-snr_min+1,2,Nr);
%绘图信息存储矩阵
graph_inf_sym=zeros(snr_max-snr_min+1,2,Nr);
for SNR=snr_min:snr_max %开始仿真
clc
disp('Wait until SNR=');disp(snr_max);
SNR
n_err_sym=zeros(1,Nr);
n_err_bit=zeros(1,Nr);
Perr_sym=zeros(1,Nr);
Perr_bit=zeros(1,Nr);
re_met_sym_buf=zeros(carrier_count,symbols_per_carrier,Nr);
re_met_bit=zeros(baseband_out_length,bits_per_symbol,Nr);
baseband_out=round(rand(baseband_out_length,bits_per_symbol)); %生成随机数用于仿真
de_data=bi2de(baseband_out); %二进制向十进制转换
data_buf=pskmod(de_data,M_psk,0); %PSK调制
carrier_matrix=reshape(data_buf,carrier_count,symbols_per_carrier);
for tt=1:Nt:symbols_per_carrier %取数为空时编码做准备，此处每次取每个子载波上连续的两个数
data=[];
for ii=1:Nt
tx_buf_buf=carrier_matrix(:,tt+ii-1);
data=[data;tx_buf_buf];
end
XX=zeros(co_time*carrier_count,Nt);
for con_r=1:co_time %进行空时编码
for con_c=1:Nt
if abs(real(O(con_r,con_c)))~=0
if imag(O(con_r,con_c))==0
XX((con_r-1)*carrier_count+1:con_r*carrier_count,con_c)=data((abs(real(O(con_r,con_c)))-1)*carrier_count+1:abs(real(O(con_r,con_c)))...
*carrier_count,1)*sign(real(O(con_r,con_c)));
else
XX((con_r-1)*carrier_count+1:con_r*carrier_count,con_c)=conj(data((abs(real(O(con_r,con_c)))-1)*carrier_count+1:abs(real(O(con_r,con_c)))...
*carrier_count,1))*sign(real(O(con_r,con_c)));
end
end
end
end %空时编码结束
XX=[tx_training_symbols;XX]; %添加训练序列
rx_buf=zeros(1,addprefix_length*(co_time+1),Nr);
for rev=1:Nr
for ii=1:Nt
tx_buf=reshape(XX(:,ii),carrier_count,co_time+1);
IFFT_tx_buf=zeros(IFFT_bin_length,co_time+1);
IFFT_tx_buf(carriers,:)=tx_buf(1:carrier_count,:);
IFFT_tx_buf(conjugate_carriers,:)=conj(tx_buf(1:carrier_count,:));
time_matrix=ifft(IFFT_tx_buf);
time_matrix=[time_matrix((IFFT_bin_length-cp_length+1):IFFT_bin_length,:);time_matrix];
tx=time_matrix(:)';
%------------------------------------------------------------------------
%d=randint(1,4,[1,7]);
%4多经信道模拟
%a=randint(1,4,[2,7])/10;
tx_tmp=tx;
d=[4,5,6,2;4,5,6,2;4,5,6,2;4,5,6,2];
a=[0.2,0.3,0.4,0.5;0.2,0.3,0.4,0.5;0.2,0.3,0.4,0.5;0.2,0.3,0.4,0.5];
for jj=1:size(d,2)
copy=zeros(size(tx)) ;
for kk = 1 + d(ii,jj): length(tx)
copy(kk) = a(ii,jj)*tx(kk - d(ii,jj)) ;
end
tx_tmp=tx_tmp+copy;
end
%------------------------------------------------------------------------
txch=awgn(tx_tmp,SNR,'measured'); %添加高斯白噪声
rx_buf(1,:,rev)=rx_buf(1,:,rev)+txch;
end
rx_spectrum=reshape(rx_buf(1,:,rev),addprefix_length,co_time+1);
rx_spectrum=rx_spectrum(cp_length+1:addprefix_length,:);
FFT_tx_buf=zeros(IFFT_bin_length,co_time+1);
FFT_tx_buf=fft(rx_spectrum);
spectrum_matrix=FFT_tx_buf(carriers,:);
Y_buf=(spectrum_matrix(:,2:co_time+1));
Y_buf=conj(Y_buf');
spectrum_matrix1=spectrum_matrix(:,1);
Wk=exp((-2*pi/carrier_count)*i);
L=10;
p=zeros(L*Nt,1);
for jj=1:Nt
for l=0:L-1
for kk=0:carrier_count-1
p(l+(jj-1)*L+1,1)=p(l+(jj-1)*L+1,1)+spectrum_matrix1(kk+1,1)*conj(tx_training_symbols(kk+1,jj))*Wk^(-(kk*l));
end
end
end
%q=zeros(L*Nt,L*Nt);
%for ii=1:Nt
% for jj=1:Nt
% for l1=0:L-1
% for l2=0:L-1
% for kk=0:carrier_count-1
% q(l2+(ii-1)*L+1,l1+(jj-1)*L+1)= q(l2+(ii-1)*L+1,l1+(jj-1)*L+1)+tx_training_symbols(kk+1,ii)*conj(tx_training_symbols(kk+1,jj))*Wk^(-(kk*(-l1+l2)));
% end
% end
% end
% end
%end
%h=inv(q)*p;
h=p/carrier_count;
H_buf=zeros(carrier_count,Nt);
for ii=1:Nt
for kk=0:carrier_count-1
for l=0:L-1
H_buf(kk+1,ii)=H_buf(kk+1,ii)+h(l+(ii-1)*L+1,1)*Wk^(kk*l);
end
end
end
H_buf=conj(H_buf');
RRR=[];
for kk=1:carrier_count
Y=Y_buf(:,kk);
H=H_buf(:,kk);
for co_ii=1:num_X
for co_tt=1:size(eta,2)
if eta(co_ii,co_tt)~=0
if coj_mt(eta(co_ii,co_tt),co_ii)==0
r_til(eta(co_ii,co_tt),:,co_ii)=Y(eta(co_ii,co_tt),:);
a_til(eta(co_ii,co_tt),:,co_ii)=conj(H(epsilon(eta(co_ii,co_tt),co_ii),:));
else
r_til(eta(co_ii,co_tt),:,co_ii)=conj(Y(eta(co_ii,co_tt),:));
a_til(eta(co_ii,co_tt),:,co_ii)=H(epsilon(eta(co_ii,co_tt),co_ii),:);
end
end
end
end
RR=zeros(num_X,1);
for iii=1:num_X %Generate decision statistics for the transmitted signal "xi"
for ttt=1:size(eta,2)
if eta(iii,ttt)~=0
RR(iii,1)=RR(iii,1)+r_til(eta(iii,ttt),1,iii)*a_til(eta(iii,ttt),1,iii)*delta(eta(iii,ttt),iii);
end
end
end
RRR=[RRR;conj(RR')];
end
r_sym=pskdemod(RRR,M_psk,0);
re_met_sym_buf(:,tt:tt+Nt-1,rev)=r_sym;
end
end
re_met_sym=zeros(baseband_out_length,1,Nr);
for rev=1:Nr
re_met_sym_buf_buf=re_met_sym_buf(:,:,rev);
re_met_sym(:,1,rev)= re_met_sym_buf_buf(:);
re_met_bit(:,:,rev)=de2bi(re_met_sym(:,1,rev));
for con_dec_ro=1:baseband_out_length
if re_met_sym(con_dec_ro,1,rev)~=de_data(con_dec_ro,1)
n_err_sym(1,rev)=n_err_sym(1,rev)+1;
for con_dec_co=1:bits_per_symbol
if re_met_bit(con_dec_ro,con_dec_co,rev)~=baseband_out(con_dec_ro,con_dec_co)
n_err_bit(1,rev)=n_err_bit(1,rev)+1;
end
end
end
end
graph_inf_sym(SNR-snr_min+1,1,rev)=SNR;
graph_inf_bit(SNR-snr_min+1,1,rev)=SNR;
Perr_sym(1,rev)=n_err_sym(1,rev)/(baseband_out_length); %Count number of error bits and symbols
graph_inf_sym(SNR-snr_min+1,2,rev)=Perr_sym(1,rev);
Perr_bit(1,rev)=n_err_bit(1,rev)/(baseband_out_length*bits_per_symbol);
graph_inf_bit(SNR-snr_min+1,2,rev)=Perr_bit(1,rev);
end
end
for rev=1:rev
x_sym=graph_inf_sym(:,1,rev); %Generate plot
y_sym=graph_inf_sym(:,2,rev);
subplot(Nr,1,rev);
semilogy(x_sym,y_sym,'k-v');
axis([2 16 0.0001 1]);
xlabel('SNR, [dB]');
ylabel('Symbol Error Probability');
grid on
%hold on
end
%hold off
%for rev=1:rev
%x_bit=graph_inf_bit(:,1,rev);
%y_bit=graph_inf_bit(:,2,rev);
%subplot(2,1,2);
%semilogy(x_bit,y_bit,'k-v');
%axis([2 16 0.0001 1]);
%xlabel('SNR, [dB]');
%ylabel('Bit Error Probability');
%grid on
%hold on
%end
%hold off
%%training_symbol.M
function tx_training_symbols=training_symbol(Nt,carrier_count)
j=sqrt(-1);
Wk=exp((-2*pi/carrier_count)*i);
training_symbols= [ 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 ...
-j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 ...
j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 ]';
tx_training_symbols=[];
for ii=1:carrier_count
training_symbols_buf=[];
for jj=1:Nt
training_symbols_buf=[training_symbols_buf,Wk^(-floor(carrier_count/Nt)*(jj-1)*ii)
*training_symbols(ii,1)];
end
tx_training_symbols=[tx_training_symbols;training_symbols_buf];
end