rician_channel_right.m

无线通信系统中莱斯信道仿真

% aa是 efso_afa, aa2=efso_afa^2;
clc
clear
close all
%%%%   rician fading channel a^2=1; k=5 10 15 
%%%%   bw_lpf1=6hz, bw_lpf2=150 hz fs=fm*4
%%%                            
aa2=0;
%angle=60;  use angle to generate aa , uz ,az
%b1=[ 0.133925,-1.6196*10^-3,1.0385*10^-5];
%b2=[0.204926,-3.9083*10^-3,2.3221*10^-5];
%if angle>40
%    aa2=b2(1)+b2(2)*angle+b2(3)*angle^2;
%else
%    aa2=b1(1)+b1(2)*angle+b1(3)*angle^2;
%end

ka=[5,10,15];  %%%% rician factor in dB
ka=ka./10;
for i=1:3      %% caculate k莱斯因子
  ka(i)=10^ka(i);
end


aa2=1;       %% a=1
A=sqrt(ka.*2*aa2);  %%
r=0:0.05:9;     %% r is  幅度
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  theory value
for k=1:3
    for i=1:length(r)
        %fr(k,i)=r(i)/aa2*exp(-(r(i)^2+A(k)^2)/(2*aa2))*besseli(0,A(k)*r(i)/aa2);
        fr(k,i)=r(i)*2/aa2*exp(-r(i)^2/aa2-ka(k))*besseli(0,2*r(i)*sqrt(ka(k)/aa2));
    end
end
figure(3)
plot(r,fr(1,:),'m-');
hold on;
plot(r,fr(2,:),'r-.');
hold on;
plot(r,fr(3,:),'--');
hold on
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  lpf2
fm=150;
f2=0:1:149;
lenf2=length(f2);
for i=1:lenf2
    Hw(i)=(1-(f2(i)/fm)^2)^(-0.5);%%%  lpf2
end
w=zeros(1,lenf2*4); %%  fs=fm*4
w(1:lenf2)=Hw;
for i=1:lenf2*2
    w(lenf2*2+i)=w(lenf2*2-i+1);
end
hn=ifft(w,lenf2*4);
hn=real(hn);                      %%% lpf2
hw=fft(hn);
figure(5)
plot(Hw,'ro');
hold on;
plot(abs(hw),'*');
energy=sum(hn.^2);    %%%%  能量归一化
hn=hn./sqrt(energy);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for k=1:3
    len=200000;
    grv2=zeros(1,len);
    grv3=zeros(1,len);
    grv2=randn(1,len);
    grv2=sqrt(aa2/2)*grv2;
    grv3=randn(1,len);
    grv3=sqrt(aa2/2)*grv3;
    %%%%%%%%%%%%%%%%%%%%%% use fft and then ifft to caculate conv
    grv2=[grv2,zeros(1,length(hn)-1)];
    grv3=[grv3,zeros(1,length(hn)-1)];%%  use to conv
    reg2=zeros(1,length(grv2));
    reg3=zeros(1,length(grv3));
    reg2(1:length(hn))=hn;            %%% use to conv
    reg3(1:length(hn))=hn;
    lpf_grv2=zeros(1,length(grv2));
    lpf_grv3=zeros(1,length(grv3));
    fre_grv2=fft(grv2);
    fre_grv3=fft(grv3);
    fre_hn=fft(reg2);
    lpf_grv2=fre_grv2.*fre_hn;
    lpf_grv2=real(ifft(lpf_grv2));
    lpf_grv3=fre_grv3.*fre_hn;
    lpf_grv3=real(ifft(lpf_grv3));

    %%%for i=1:length(grv2)
    %%%    reg2(2:length(hn))=reg2(1:length(hn)-1);
    %%%    reg2(1)=grv2(i);
    %%%    reg3(2:length(hn))=reg3(1:length(hn)-1);
    %%%    reg3(1)=grv3(i);
    %%%    lpf_grv2(i)=0;
    %%%    lpf_grv3(i)=0;
    %%%    for j=1:length(hn)
    %%%    lpf_grv2(i)=hn(j)*reg2(j)+lpf_grv2(i);
    %%%    lpf_grv3(i)=hn(j)*reg3(j)+lpf_grv3(i);
    %%%   end
    %%%end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    delay2=round(length(hn)/2);%% delay in the result
    los=ka(k)*aa2;             %% los signal
    e_grv1=sqrt(los);

    output=zeros(1,len);
    outputj=zeros(1,len);
    output=lpf_grv2(delay2:len+delay2-1)+e_grv1;
    outputj=lpf_grv3(delay2:len+delay2-1)*j;
    output=output+outputj;
    range=zeros(1,len);
    range=abs(output);

    range=quant(range,0.05);%%%%%   simulate r
    range1=round(range/0.05); %%% there are some value  not int,which lead to error when compute pdf;
    num=zeros(1,length(r));
    for i=1:len
        m=range1(i);
        if m>160
            m=160;
        end
        num(m+1)=num(m+1)+1;
    end
    %%%  pdf
    rate(k,:)=zeros(1,length(r));
    rate(k,:)=num./len*20;
end
figure(3)
plot(r,rate(1,:),'m*');
hold on;
plot(r,rate(2,:),'ro');
hold on;
plot(r,rate(3,:),'+');
hold on;
figure
dop = abs(output);
dop = 10*log(dop);
grv = abs(grv2+grv3*j);
grv = 10*log(grv);
subplot(2,1,1)
plot(dop);
subplot(2,1,2)
plot(grv);