基于MATLAB的MIMO通信系统仿真.doc

目 录
（一）基于MATLAB的MIMO通信系统仿真…………………………
一、基本原理………………………………………………………
二、仿真……………………………………………………………
三、仿真成果………………………………………………………
四、仿真成果分析…………………………………………………
（二）自选习题部分…………………………………………………
（三）总结与体会……………………………………………………
（四）参照文献……………………………………………………
实训报告
（一）基于MATLAB的MIMO通信系统仿真
一、 基本原理
二、 仿真
三、 仿真成果
四、 仿真成果分析

OFDM技术通过将频率选择性多径衰落信道在频域内转换为平坦信道，减小了多径衰落的影响。OFDM技术如果要提高传播速率，则要增长带宽、发送功率、子载波数目，这对于频谱资源紧张的无线通信时不现实的。
MIMO可以在空间中产生独立并行信道同步传播多路数据流，即传播速率很高。这些增长的信道容量可以用来提高信息传播速率，也可以通过增长信息冗余来提高通信系统的传播可靠性。但是MIMO却不可以克服频率选择性深衰落。
因此OFDM和MIMO这一对互补的技术自然走到了一起，目前是3G，将来也是4G，以及新一代WLAN技术的核心。总之，是核心物理层技术之一。

1、MIMO系统理论：
核心思想：时间上空时信号解决同空间上分集结合。
时间上空时通过在发送端采用空时码实现： 空时分组、空时格码，分层空时码。
空间上分集通过增长空间上天线分布实现。此举可以把本来对顾客来说是有害的无线电波多径传播转变为对顾客有利。

2、MIMO系统模型：

可以看到，MIMO模型中有一种空时编码器，有多根天线，其系统模型和上述MIMO系统理论一致。为什么说nt>nr，由于一般来说，移动终端所支持的天线数目总是比基站端要少。
接受矢量为：，即接受信号为信道衰落系数X发射信号+接受端噪声

3、MIMO系统容量分析：
（附MIMO系统容量分析程序）
香农公式的信道容量（即信息传送速率）为：

4、在MIMO中计算信道容量分两种状况：
未知CSI和已知CSI（CSI即为信道状态信息），其公式推导较为复杂，推导成果为信道容量是信噪比与接受、发射天线的函数。
在推导已知CSI中，常用的有waterfilling，即出名的注水原理。但是，根据有关文献资料，一般状况下CSI可以当做已知，由于发送，接受端会根据具体信道状况估算CSI的有关参数。
在这里对注水原理做一种简朴简介：之因此成为注水原理是由于抱负的注水原理是在噪声大的时候少分派功率，噪声小时多分派功率，最后噪声+功率=定值，这如果用图形来表达，则类似于给水池注水的时候，水池低的地方就多注水，也就是噪声小分派的功率就多，故称这种达到容量的功率分派方式叫做注水原理。通过给各个天线分派不同的发射功率，增长系统容量。核心思想就是上面所论述的，信道条件好，则分派更多功率；信道条件差，则分派较少的功率。
在MIMO的信道容量当中要注意几种问题：（下面说已知CSI都是加入了估计CSI的算法，并且采用了注水原理。）
1. 已知CSI的状况下的信道容量要比发送端未知CSI的状况下的信道容量高，这是由于当发送端已知CSI的时候，发送端可以优化发送信号的协方差矩阵。也就是可以通过注水原理使得信道容量达到最大。因此在实际系统当中，发射端必须有效运用CSI（如上述阐明的一般采用估算的措施），这样可以优化发送信号。
2. 如果信道信噪比足够大的时候，这时已知和未知CSI相差不大。由于已知CSI状况下此时所有的池子都是满的，此时发送天线上的功率得到了最大分派。和未知CSI状况同样。
3. 如果收发天线数相等，这时已知CSI比未知CSI的优势并不明显；但是，当发射天线数不小于接受天线数时，已知CSI的信道容量要明显高于未知CSI。

5、MIMO和OFDM结合使用，即MIMO-OFDM系统：

运用MIMO和OFDM结合，可以大大提高无线通信系统的信道容量和传播速率，有效抵御信道衰落和克制干扰，被觉得是构建宽带无线通信系统最核心的物理层传播方案。
从图中可以看出，MIMO-OFDM系统中，每根发射天线的通路上均有一种OFDM调制器，每根接受天线上均有一种OFDM解调器。

空时编码技术之STBC简介：
在上文的论述中讲到了MIMO通过时间上空时信号解决和空间上分集结合实现。空时信号解决，即空时编码技术目前研究较多的是分层空时编码（BLAST），空时网格编码（STTC）以及空时分组编码（STBC）。其中，空时分组编码（STBC）与1998年发明，应用最为广泛。

星座映射的解释如下：
数字调制用"星座图"来描述，星座图中定义了一种调制技术的两个基本参数：1）信号分布；2）与调制数字比特之间的映射关系。星座图中规定了星座点与传播比特间的相应关系，这种关系称为"映射"，一种调制技术的特性可由信号分布和映射完全定义，即可由星座图来完全定义。将输入的串行二进制信息序列经串－并变换，变成m=log2M个并行数据流，M为星座图的星座点数目，每一路的数据率是R/m,R是串行输入码的数据率。每m个比特相应星座图上的的一种星座点，例如BPSK调制，每1个比特相应一种星座点；QPSK调制，每2个比特相应一种星座点；16QAM调制，每4个比特相应一种星座点。一般映射采用格雷码为映射规则。
在此例中，信源发送的二进制信息比特一方面进行星座映射。假设采用4进制的调制星座，有。把从信源来的二进制信息比特每2个比特分为一组（此例中是x1和x2），对持续的两组比特进行星座映射，得到两个调制符号x1,x2。把这两个符号送入编码器，并按照如下方式编码：

在第一种发送时刻，符号在天线1上发送出去，符号在天线2上发送出去。第二个时刻，符号在天线1上发送出去，符号在天线2上发送出去。可以看出，两幅发送天线上发送信号批次存在着一定的关系，因此这种空时码是基于发送分集的。
两幅发送天线上发送的信号满足正交特性。

考虑两个发送天线，一种接受天线的状况：
假设接受端可以完全精确地估计出信道的衰落系数和，在接受端采用最大似然估计，从星座中找出一对符号，该符号也就是最后接受端觉得发送端发送的符号。在程序中，其判决式为：

其中，是根据信道衰落系数和接受信号进行合并得到的信号。

考虑多接受天线的状况：
多天线系统中，发送端的编码与传播方案和单接受天线系统同样。只是在接受端的解决变得复杂，需要对不同接受天线上接受到的信号进行合并解决。多接受天线下的判决度量可以通过把各副接受天线上的接受信号得到的判决度量线性合并得到。判决式如下（该判决式在程序中有体现）：

最后的系统构造图：

程序：
信道容量：
close all;
clear all;
clc;
N_loop=100;
EbN0=0:5:25;
N0=1;
Eb=10.^(EbN0/10)*N0;
normalized=sqrt(1/2);
C_norm=zeros(1,length(EbN0));
D=zeros(1,N_loop);
n=[1,2,4,8];
W=3*10^4;
for i=1:4
for ee=1:length(EbN0)
for s=1:N_loop
H=randn(n(i),n(i))+j*randn(n(i),n(i));
Q=H'*H;
m=n(i);
I=eye(m,m);
SNR=Eb(ee)/N0;
C=W*log2(det(I+SNR/n(i)*Q));
D(s)=C;
end
C_avg=sum(D)/N_loop;
C_normal=C_avg/W;
C_norm(ee)=C_normal;
disp(['When SNR is ',num2str(EbN0(ee)),'dB: C_avg',int2str(i),' = ',num2str(C_avg),'; C_normal ',int2str(i),' = ',num2str(C_normal)]);
end
P=plot(EbN0,C_norm,'r-d');
set(P,'linewidth',[2]);
axis([0 30 0 100]);
AX=gca;
set(AX,'fontsize',14);
title('\fontsize{12}\bfMIMO ÐŵÀÈÝÁ¿');
X=xlabel('Eb/N_0 [dB]£¨ÐÅÔë±È£©');
set(X,'fontsize',12);
Y=ylabel('Capacity bit/s/Hz£¨ÐŵÀÈÝÁ¿£©');
set(Y,'fontsize',12);
hold on;
grid on;
end
legend('n1=1','n2=2','n3=4','n4=8');

误码率与信噪比关系：
clear all
close all
clc

%+++++++++++++++++++++++++++±äÁ¿++++++++++++++++++++++++++++++
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 -3;2+j 1+j 0;3+j 0 1+j;0 -3+j 2+j];
co_time=size(O,1);
Nt=size(O,2);
Nr=2;
disp('--------------start-------------------');
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
for con_co=1:Nt
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));
conjugate_carriers=IFFT_bin_length-carriers+2;
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);
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(:)';
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
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=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
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);
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);
y_sym=graph_inf_sym(:,2,rev);
subplot(Nr,1,rev);
semilogy(x_sym,y_sym,'b-*');
axis([2 16 0.0001 1]);
xlabel('ÐÅÔë±È/dB');
ylabel('ÎóÂëÂÊ');
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
disp('--------------end-------------------');
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% ************************beginning of file*****************************
%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
%************************end of file**********************************

（二）自选习题部分
实例二：
clear;
reg=ones(1,9);
coeff=[1,0,0,1,0,1,1,0,0,1];
N=2^length(reg)-1;
for k=1:N
a_n=mod(sum(reg.*coeff(1:length(coeff)-1)),2);
reg=[reg(2:length(reg)),a_n];
out(k)=reg(1);
end
out=2*out-1;
for j=0:N-1
rho(j+1)=sum(out.*[out(1+j:N),out(1:j)])/N;
end
j=-N+1:N-1;
rho=[fliplr(rho(2:N)),rho];
plot(j,rho);axis([-10 10 -0.1 1.2]);
运营成果：

实例三：
clear;
reg=ones(1,6); 寄存器初始状态，全1，寄存器级数为6
coeff=[1,0,0,0,0,1,1]; 抽头数cr,---,c1,c0取决于特性多项式
N=2^length(reg)-1; 周期
for k=1:N 计算一种周期的m序列输出
a_n=mod(sum(reg.*coeff(1:length(coeff)-1)),2); 反馈
reg=[reg(2:length(reg)),a_n]; 寄存器移位，反馈
out