% Z = plot3Dfigure(r,a,phi,tao,L)
% -------------------
% 功能：繪制功率－延遲－角度分布圖。
% 輸入：
% 1、歸一化多徑延遲：r
% 2、各多徑分量復幅度：a
% 3、波達方向：phi
% 4、LOS分量時延tao 
% 5、多徑數L
% 輸出：
% 信道總的沖激響應：Z

function Z = plot3Dfigure(r,a,phi,tao,L)

r = [tao*1e6,r];    % 將歸一化多徑延遲轉換為實際時延
r = round(r /0.005);    % 將時延r量化，令分辨率為5000ns
a = [1,a];
phi = [0,phi] * 180 / pi;   % 將波達方向由弧度換算成角度
phi = round(phi /1);    % 將波達方向量化，令分辨率為1度

h = Gauss3DPuls(0.05,10,0.02,4);    % 生成一合適的高斯脈沖，令時間寬度為500ns，角度寬度為10度
maxh = max(max(h)); % 記錄此高斯脈沖的最大值

[X,Y] = meshgrid(0-0.05:0.005:3.33+0.05, -180-25:1:180+25); % 生成TOA-DOA坐標“格點”矩陣
Z=zeros(361+50,677+10); 

% 用各徑功率對高斯脈沖進行加權并將其疊加，得到總的信道沖擊響應
for i=1:L
    temph = h * abs(a(i));
    tempZ = zeros(361+50,677+10);
    tempZ(phi(i)-10+180+25:phi(i)+10+180+25,r(i)-10+10:r(i)+10+10) = temph;  
    Z = Z + tempZ;
end

% 將閾值設定為比直射路徑分量功率低20dB，所有低于這個閾值到達的多徑分量都被忽略
TempZ = Z;
NotZeroZ = TempZ~=0;
TempZ(NotZeroZ) = 10*log10(Z(NotZeroZ).^2/max(max(TempZ))^2) + 20;
LessthanZeroZ = TempZ<=0;
TempZ(LessthanZeroZ) = 0;

figure;
mesh(X,Y,TempZ),colormap('jet');
axis([0 3.4 -195 195 0 max(max(TempZ))+1]) 
colormap(jet)
view(67,22)
xlabel('TOA(微秒)')
ylabel('DOA(度)')
zlabel('功率門限比(dB)')