% =======================================================================
% project722    (includes scenario editor)
% =======================================================================
% Initialize ============================================================
clear
close all
clc
warning off 
% basic inputs ==========================================================

fc=2000;         % MHz  Carrier frequency
F=16;            % sampling rate: fraction of wave length
V=10;            %  m/s MS1 speed 
NFFT=1024;       % Number of points in FFT
Nsamples=200;    % Number of samples 
avPower=-20;     % sigma^2  Raverage power

% geometry inputs ========================================================

dBS=1000;     
angleBS=130;
BSx=dBS*cosd(angleBS) % location of transmitter (BS) x-coordinate
BSy=dBS*sind(angleBS)  % location of transmitter (BS) y-coordinate

% locations of point scatterers =========================================

fig=figure;
plot(BSx,BSy,'k^'), hold on

% indirect parameters ===================================================

lambdac=300/fc;    % m wavelength
Dx=lambdac/F;      % m sampling spacing 
ts=Dx/V;           % s time sampling interval
fs=1/ts;           % Hz sampling frequency
kc=2*pi/lambdac;   % propagation constant
fm=V/lambdac       % max Doppler shift
cc=3e8;            % speed of light
timeaxis=ts.*[0:Nsamples-1];

% ========================================================================
MS0=-V*timeaxis(end)/2;   % initial location of receiver (MS) x-coordinate

MSx=MS0+V.*timeaxis;    % MS route along x-axis
MSy=zeros(Nsamples,1)';  % MS route along x-axis (y=0)
plot(MSx,MSy,'k','LineWidth',5)

MINx=min(min(BSx,MSx))-500;
MAXx=max(max(BSx,MSx))+500;
MINy=min(min(BSy,MSy))-500;
MAXy=max(max(BSy,MSy))+500;
axis([MINx MAXx MINy MAXy])
plot([0 0],[MINy MAXy], 'k:')
plot([MINx MAXx],[0 0], 'k:')

%=========================================================================
% SCENARIO EDITOR 
% ========================================================================
% placing point-scatterers in propagation scenario.

[SCx,SCy] = getpts(fig);
NSC=length(SCx);
plot(SCx,SCy,'k+');

% =======================================================================
% calculate distance matrix 
% =======================================================================

distBSSC=sqrt((BSx-SCx).^2+(BSy-SCy).^2);

distBSSCext=repmat(distBSSC,1,Nsamples);

distSCMS=zeros(NSC,Nsamples);
for ii=1:Nsamples
    distSCMS(:,ii)=sqrt((SCx-MSx(ii)).^2+SCy.^2);
end

distBSSCMS=distBSSCext+distSCMS;

% ======================================================================

distBSMS1aux=sqrt((BSx-MSx).^2+(BSy-MSy).^2);   
distBSMS1=min(min(distBSMS1aux));               % Ref distance is min BSMS dist 

% level of scattered contributions

a=(distBSMS1./distBSSC(:)).*(1./distSCMS(:,1));
DeltaPower=avPower-10*log10(sum(a.^2));
deltaa=10.^(DeltaPower/20);             % to achieve reference power
a=deltaa*a;

% =====================================================================
% Define time-varying complex magnitudes of point scatterer contributions 
% amplitudes remain constant while phases change

aa=zeros(Nsamples,NSC);     % create variable 

for k1=1:Nsamples           % scan route points
    for k2=1:NSC            % scan scatterers
        aa(k1,k2)=a(k2)*exp(-j*kc*distBSSCMS(k2,k1));  % time-varying phase
    end
end
 
% ======================================================================

distBSSCMS1=distBSSCMS-distBSMS1;     % set a new refernece for delays wrt to 
DelaysNormalized=distBSSCMS1/cc;      % arrival of direct ray, here assumed 
                                      % to be totally blocked 

DelaysNormalized=round(DelaysNormalized*1e9);  % convert quantify delays to 1 ns
%figure,mesh(DelaysNormalized)
auxx=size(DelaysNormalized);

auxx2=max(max(DelaysNormalized))+1;    % to include 0 ns delay
DelayProfile=zeros(auxx(2),auxx2);     % Create delay profile with step of 1 ns

for jj=1:auxx(2)           % scan route locations
    for ii=1:auxx(1)       % scan scatterers
        indexx=DelaysNormalized(ii,jj)+1;               %<--------- ????
        DelayProfile(jj,indexx)=DelayProfile(jj,indexx)+aa(jj,ii);     
                                    % put in corresponding delay bin 
                                    % complex amplitude of delta
    end
end

axisdelayprofile=[0:auxx2-1];   % axis in ns
figure, hold
for ii=1:Nsamples
 stem(axisdelayprofile,abs(DelayProfile(ii,:)))  % accumulate deltas with time and delay on same plot
%  pause
end

% ======================================================================= 
% Parameters of ideal PDP

% PDP parameters: D, S 
%========================================================================
D=sum(abs(DelayProfile(1,:)).*axisdelayprofile)/sum(abs(DelayProfile(1,:)))  
% 1st route point
S=sqrt(sum(abs(DelayProfile(1,:)).*(axisdelayprofile-D).^2)/sum(abs(DelayProfile(1,:))))
%=======================================================================
% Simulate channel sounding
%=======================================================================
% build sounding pulse

SequenceLength=511;      % length of sounding sequence
tau_sampl_spacing=1;     % in ns
tau_chip=100;            % in ns   1/Chip rate of sounding sequence
SamplesInImpulse=round(tau_chip/tau_sampl_spacing);
tau_chip_mod=SamplesInImpulse*tau_sampl_spacing;
% SlopePulse=(SequenceLength+1)/(SamplesInImpulse+1);
SlopePulse=SequenceLength/(SamplesInImpulse+1);
% SoundImpulse=[0:SamplesInImpulse]*SlopePulse-1;
SoundImpulse=[0:SamplesInImpulse]*SlopePulse;
%SoundImpulseAux=[SamplesInImpulse-1:-1:0]*SlopePulse-1;
SoundImpulseAux=[SamplesInImpulse-1:-1:0]*SlopePulse;
SoundImpulse=[SoundImpulse SoundImpulseAux];

figure, plot([0:length(SoundImpulse)-1],SoundImpulse)

% Normilize sound impulse to unit energy
ESoundImpulse=sum(SoundImpulse.^2);
SoundImpulse=SoundImpulse/sqrt(ESoundImpulse);
figure, plot([0:length(SoundImpulse)-1],SoundImpulse)


figure,hold
pdp=[];
for ii=1:Nsamples,
    pdpaux=conv(SoundImpulse,DelayProfile(ii,:));
    auxx3=length(pdpaux);
    plot([0:auxx3-1]-tau_chip,abs(pdpaux)) 
    pdp=[pdp, pdpaux'];
    %pause
end

figure, plot([0:auxx3-1]-tau_chip,10*log10(abs(pdpaux)))

figure,surf(timeaxis, [0:auxx3-1]-tau_chip,abs(pdp))    
shading interp, colormap(jet) 
view(2)
figure,surf(timeaxis, [0:auxx3-1]-tau_chip,abs(pdp))    
shading interp, colormap(jet) 

figure,surf(timeaxis, [0:auxx3-1]-tau_chip,10*log10(abs(pdp)))
shading interp, colormap(jet) 
view(2)
figure,surf(timeaxis, [0:auxx3-1]-tau_chip,10*log10(abs(pdp)))
shading interp, colormap(jet) 

% =======================================================================
% Compute averaged power delay profile APDP from instantaneous PDP
% =======================================================================
apdp=sum(abs(pdp'));    % trapose pdp since sum operates row-wise
apdp=apdp/Nsamples;      % 

figure, plot([0:auxx3-1]-tau_chip,apdp)

figure, plot([0:auxx3-1]-tau_chip,10*log10(apdp))

% =======================================================================
% Compute parameters of APDP
% =======================================================================
D=sum(apdp.*[0:auxx3-1])/sum(apdp)
S=sqrt(sum(apdp.*([0:auxx3-1]-D).^2)/sum(apdp))