tu1-31.m

我自己在学习空时编码一书时编写的第一章的仿真程序

%=================关于本程序的说明和学习心得====================================================================
%本程序中有对波长λ的计算,其实根据eq1-103关于相关系数的计算中,被积函数exp(-j*2*pi*d_mk/λ*sin(α))中,d_mk/λ可
%以直接把λ约掉,但是为什么这里还是要计算λ呢?
%注意到的是，本程序中还有对散射体相关矩阵的计算,散射体的个数和接收天线的位置个数是不一样的，导致了接受天线相关矩
%阵和发射天线相关矩阵与散射体相关矩阵的维数不一样,我们是用同一个函数来计算这些相关矩阵的,为了编写计算相关矩阵的函
%数(sitafun())的方便，所以把他们各自的长度表示出来,这也就是为什么程序中要计算波长λ,继而计算接收天线和发射天线位
%置分布(程序中第22,23,25行)和散射体的位置分布的原因,位置分布为随机产生的均匀分布变量,应该要通过求统计平均的方式计算,
%程序见tu131(1).m和tu131(2).m(注意:函数文件sitafun()中用x代表位置分布)
%==========================================================================
clear
Dr=50;%50m
Dt=Dr;
Rr=Dr;
Rt=Dr;
R=1000e3;%1000km
SNRdB=20;
SNR=10.^(SNRdB./10);
alfa_r=2*atan(Dr./Rr);       %接收端角度扩展
alfa_t=2*atan(Dt./Rt);       %发射端角度扩展
alfa_s=2*atan(Dt./R);        %接收散射体的角度扩展
f=900e6;%900MHz
lamda=3e8/f;%wave length       %波长λ的计算
Transantannaposition=0:3*lamda:9*lamda;   %传输天线位置分布（0   3λ   6λ   9λ）
Receiantannaposition=0:3*lamda:9*lamda;   %接收天线位置
S=32;%Scatter bodies number               %散射体数量
Scatterposition=rand(1,S).*2*Dt;          %散射体位置分布   满足均匀分布
Sita_r=sitafun(Receiantannaposition,lamda,alfa_r);    %接收天线相关矩阵
Sita_t=sitafun(Transantannaposition,lamda,alfa_t);    %传输天线的相关矩阵
Sita_s=sitafun(Scatterposition,lamda,alfa_s);         %散射体相关矩阵
Kr=(llh(Sita_r));                                    %对接收天线的相关矩阵进行复Cholesky分解，自编函数llh(),
                                                      %如果是matlab7.2版本中的chol函数就可以直接完成这个功能
Kt=(llh(Sita_t));
Ks=(llh(Sita_s));
simutimes=10000
for sim=1:simutimes
    Gr=(randn(4,S)+j*randn(4,S))./sqrt(2);           %独立的复高斯随机变量
    Gt=(randn(S,4)+j*randn(S,4))./sqrt(2);
    H=1./sqrt(S).*(Kr*Gr*Ks*Gt*Kt);%eq1-141
    C(sim)=log2(det(eye(4,4)+SNR./4.*H'*H));%1-30
end
dx=0.01;
[c,xout]=hist(real(C),0:dx:30);
pdf=c./simutimes./dx; %pdf fun
plot(xout,1-cumsum(pdf).*dx);
axis([0 30 0 1]);grid on;hold on;
%================================================(d)
SNRdB=20;
SNR=10.^(SNRdB./10);
simutimes=10000
for sim=1:simutimes
H=(randn(4,4)+j*randn(4,4))./sqrt(2);
C(sim)=log2(det(eye(4,4)+SNR./4.*H'*H));%1-30
end
dx=0.01;
[c,xout]=hist(real(C),0:dx:30);
pdf=c./simutimes./dx; %pdf fun
plot(xout,1-cumsum(pdf).*dx);
hold on; grid on;
xlabel('容量(bit/s/Hz)');
ylabel('大于横坐标的概率 [容量>横坐标的概率]');
text(7,0.52,'a');
text(12,0.52,'b');
text(21,0.52,'c');
text(22.5,0.52,'d');
%===============================================