distrib.m

来自「本函数计算一个大区域上的信号干扰分布情况。」· M 代码 · 共 79 行

%（论文）这个程序的功能是在已知了几个原有用户位置的情况下，计算考察频点在区域上的分布（自由空间传播模型）
%lo表示原有用户的位置信息
%p_f行表示原有用户，列表示考察频段
function distrib(lo,p_f)   
LoN=size(lo,1);   
D=6e3;               %区域半径
fs=1e9;
f0=4e8;             %初始频率
B=2e5;                %考察小频段为200K
% deltf=fs/(2.^16);
span=50;
Q=zeros(2*D/span+1,2*D/span+1);
% k=1;
% 
% %     A=(sqrt((x-lo(1,1)).^2+(y-lo(1,2)).^2))>2e1;
% %     for i=2:LoN
% %         B=(sqrt((x-lo(i,1)).^2+(y-lo(i,2)).^2))>2e1;
% %         A=A&B;
% %     end
    for x=-D:span:D
        for y=-D:span:D
            for n=1:LoN
                 if((x~=lo(n,1))&(y~=lo(n,2)))
%                      am=(1.1*log10((f0+k*B)/1e6)-0.7)*3-(1.56*log10((f0+k*B)/1e6)-0.8);                    
%                      Lm((x+D)/100+1,(y+D)/100+1)=69.55+26.16*log10((f0+k*B)/1e6)-13.82*2-am...
%                          +(44.9-6.55*2)*log10(((sqrt((x-lo(n,1)).^2+(y-lo(n,2)).^2)))/1000);
%                      Lm((x+D)/100+1,(y+D)/100+1)=10.^(Lm((x+D)/100+1,(y+D)/100+1)/10);
%                      P((x+D)/100+1,(y+D)/100+1,n)=p_f(n,k)./Lm((x+D)/100+1,(y+D)/100+1);
                         P((x+D)/span+1,(y+D)/span+1,n)=p_f(n)*(3e8/4/pi).^2/((lo(n,1)-x)^2+(lo(n,2)-y)^2)/(f0^2);
                 else
%                      am=(1.1*log10((f0+k*B)/1e6)-0.7)*3-(1.56*log10((f0+k*B)/1e6)-0.8);                    
%                      Lm((x+D)/span+1,(y+D)/100+1)=69.55+26.16*log10((f0+k*B)/1e6)-13.82*2-am...
%                          +(44.9-6.55*2)*log10(((sqrt((x+10-lo(n,1)).^2+(y+10-lo(n,2)).^2))+100)/1000);  %+90是为了不使log函数的参数为0
%                      Lm((x+D)/100+1,(y+D)/100+1)=10.^(Lm((x+D)/100+1,(y+D)/100+1)/10);
%                      P((x+D)/100+1,(y+D)/100+1,n)=p_f(n,k)./Lm((x+D)/100+1,(y+D)/100+1);
                         P((x+D)/span+1,(y+D)/span+1,n)=p_f(n)*(3e8/4/pi).^2/((lo(n,1)-x)^2+(lo(n,2)-y)^2+80)/(f0^2);
%                          P((x+D)/span+1,(y+D)/span+1,n)=P((x+D)/span,(y+D)/span,n)
                 end
            end
        end
    end
    Q=sum(P,3);
    Q=db(Q,'power');
figure;
colormap('gray');
[x,y]=meshgrid(-D:span:D,-D:span:D);
Q=Q';                   %画图时总出现x和y颠倒的问题，不能解决， 
mesh(x,y,Q);       % 故用mesh(y,x,Q)代替mesh(x,y,Q),或用mesh(x,y,Q')代替mesh(x,y,Q)
    xlabel('x');
    ylabel('y');
    zlabel('干扰温度:dB');
%     sum(size(Q))
    Limit=sum(sum(Q))/(size(Q,1)*size(Q,2));
%     Limit
% Limit=-span;
    Dt=zeros(size(Q));
    for i=1:size(Q,1)
        for j=1:size(Q,2)
            if Q(i,j)>=Limit+5
                Dt(i,j)=0;
            elseif Q(i,j)>=Limit
                  Dt(i,j)=1;  
                else Dt(i,j)=2;
            end
        end
    end
    figure;
%     Dt
colormap('gray');
 pcolor(x,y,Dt);
 shading interp;
  xlabel('x');
 ylabel('y');
 title('白色区域：白空，灰色区域：灰空，黑色区域：黑空');
% mesh(x,y,Dt)
% end