read data

% data format:
%   1) link(TX position, line, RX position, frequency band), where frequency
% band 1: 6.75 GHz, 2: 30 GHz, 3: 60 GHz
%   2) valid_idx 1: this measurement can be used, 0: this measurement is
% corrupted
%
% Details:
% link
%     Value (1x4053x2x7x4x1x1): CIR (snapshot,delay,polarization,tx az.
%     scan, tx el. scan, rx az. scan, rx el. scan)
%     txAz 1x7: angles in the TX scans in azimuth
%     txEl 1x4: angles in the TX scans in elevation
%     rxAz 1x1: angles in the RX scans in azimuth
%     rxEl 1x1: angles in the RX scans in elevation
%     Dim 1x7: dimensions of each 
%     TX_id: TX label
%     TX_xyz: TX position (x,y,z) (see map)
%     RX_id: RX label
%     RX_xyz: RX position (x,y,z) (see map)
%     Visibility: LOS, NLOS

%% read data
close all; clear all; clc;

load('Bosch')

frequency_label = {'6.75 GHz','30 GHz','60 GHz'};

%% calculate TX RX distance
clc; close all;
for tt = 1:size(link,1)    
    for ll = 1:size(link,2)
        for rr = 1:size(link,3)
            for bb = 1:size(link,4)            
                if valid_idx(tt,ll,rr,bb) == 1
                    TXRX_distance(tt,ll,rr,bb) = sqrt(sum(abs(link(tt,ll,rr,bb).TX_xyz - link(tt,ll,rr,bb).RX_xyz).^2));                    
                else
                    TXRX_distance(tt,ll,rr,bb) = NaN;
                end
            end
        end
    end
end

%% Multi-band PDP per position

close all; clc;
clear PDP PAP

%define dR for visualization
DR_dB = 30;
% axis are from a valid measurement
tau = link(1,1,10,1).Dim(2).Value(1:4053).*1E-9;

c0 = 299792458;

for tt = 1:size(link,1)    
    for ll = 1:size(link,2)
        for rr = 1:size(link,3)
            for bb = 1:size(link,4)
                        
                if valid_idx(tt,ll,rr,bb)
                    
                    PDP(tt,ll,rr,bb,:) = squeeze(10*log10( ...
                        mean(abs(link(tt,ll,rr,bb).Value(:,:,1,:)).^2,4)));                                       
                    RX_id{tt,ll,rr,bb} = link(tt,ll,rr,bb).RX_id;                    
                else                    
                    PDP(tt,ll,rr,bb,:) = nan(length(tau),1);                    
                    RX_id{tt,ll,rr,bb} = NaN;
                end
            end
            
            if sum(squeeze(valid_idx(tt,ll,rr,:)))
                
                figure
                aux = squeeze(PDP(tt,ll,rr,1,:));
                aux = aux-max(aux(:));
                aux(aux < -DR_dB) = -DR_dB;
                plot(tau.*c0,aux)
                hold on
                aux = squeeze(PDP(tt,ll,rr,2,:));
                aux = aux-max(aux(:));
                aux(aux < -DR_dB) = -DR_dB;                
                plot(tau.*c0,aux)
                aux = squeeze(PDP(tt,ll,rr,3,:));
                aux = aux-max(aux(:));
                aux(aux < -DR_dB) = -DR_dB;
                plot(tau.*c0,aux)
                axis([0 120 -DR_dB 0]);
                hold off
                grid on
                lH = legend( ...
                    [frequency_label{1} ' - ' RX_id{tt,ll,rr,1}], ...
                    [frequency_label{2} ' - ' RX_id{tt,ll,rr,2}], ...
                    [frequency_label{3} ' - ' RX_id{tt,ll,rr,3}])

                pause(1)
            end
        end
    end
end

%% PDP per line

close all; clc;
clear PDP 

%define dR for visualization
DR_dB = 30;
% axis are from a valid measurement
tau = link(1,1,10,1).Dim(2).Value(1:4053).*1E-9;

c0 = 299792458;

for tt = 1:size(link,1)    
    for ll = 1:size(link,2)
        for rr = 1:size(link,3)
            for bb = 1:size(link,4)
                        
                if valid_idx(tt,ll,rr,bb)
                    
                    PDP(tt,ll,rr,bb,:) = squeeze(10*log10( ...
                        mean(abs(link(tt,ll,rr,bb).Value(:,:,1,:)).^2,4)));
                    RX_id{tt,ll,rr,bb} = link(tt,ll,rr,bb).RX_id;
                    TXRX_d(tt,ll,rr,bb) = TXRX_distance(tt,ll,rr,bb);                                       
                    
                else
                    PDP(tt,ll,rr,bb,:) = nan(length(tau),1);
                    RX_id{tt,ll,rr,bb} = NaN;
                    TXRX_d(tt,ll,rr,bb) = NaN;
                end
            end                        
        end
                                
        %we limit the samples in delay domain to 120 m propagation
        %distance

        tau_idx = 1:find(tau.*c0 < 120,1,'last');

        figure            
        bb_idx = 1;
        aux = squeeze(PDP(tt,ll,:,bb_idx,tau_idx));
        aux = aux-max(aux(:));
        aux(aux < -DR_dB) = -DR_dB;
        surf(tau(tau_idx).*c0,squeeze(TXRX_d(tt,ll,:,bb_idx)),aux)
        shading interp
        axis([0 120 0 max(squeeze(TXRX_d(tt,ll,:,bb_idx))) -DR_dB 0]);
        grid on
        xlabel('Propagation distance (m)');
        ylabel('TX-RX distance (m)');
        zlabel('Norm. RX Power (dB)');
        view(25,75)
        title([frequency_label{bb_idx} ' - Line ' num2str(ll)]);

        figure
        bb_idx = 2;
        aux = squeeze(PDP(tt,ll,:,bb_idx,tau_idx));
        aux = aux-max(aux(:));
        aux(aux < -DR_dB) = -DR_dB;
        surf(tau(tau_idx).*c0,squeeze(TXRX_d(tt,ll,:,bb_idx)),aux)
        shading interp
        axis([0 120 0 max(squeeze(TXRX_d(tt,ll,:,bb_idx))) -DR_dB 0]);
        grid on
        xlabel('Propagation distance (m)');
        ylabel('TX-RX distance (m)');
        zlabel('Norm. RX Power (dB)');
        view(25,75)
        title([frequency_label{bb_idx} ' - Line ' num2str(ll)]);

        figure
        bb_idx = 3;
        aux = squeeze(PDP(tt,ll,:,bb_idx,tau_idx));
        aux = aux-max(aux(:));
        aux(aux < -DR_dB) = -DR_dB;
        surf(tau(tau_idx).*c0,squeeze(TXRX_d(tt,ll,:,bb_idx)),aux)
        shading interp
        axis([0 120 0 max(squeeze(TXRX_d(tt,ll,:,bb_idx))) -DR_dB 0]);
        grid on
        xlabel('Propagation distance (m)');
        ylabel('TX-RX distance (m)');
        zlabel('Norm. RX Power (dB)');
        view(25,75)
        title([frequency_label{bb_idx} ' - Line ' num2str(ll)]);

        pause        
    end
end