init_aeronauticalchannel.m

国外发布的航空信道matlab模型

function state=Init_AeronauticalChannel(SampFreq,StartAltitude,DescendingRate,AircraftType,MarcovProbs,MarkovAtten,Altitudevector,Elevation,Azimuth)

% Disclaimer: The channel simulation software "aeronautical channel model" has been
% developed by DLR under ESA-Contract.  The use of this software is free of charge. 
% Anyone using this software, will use it exclusively at his/her own risk and responsibility.
% Neither ESA nor DLR nor their partners, suppliers, affiliates nor any other party involved
% in creating, producing, or delivering the software is liable for any direct, indirect,
% incidental, special, consequential, punitive or other damages whatsoever including business
% interruption, loss of use, data, information, profits (regardless of the form of action,
% including but not limited to contract, negligence or other tortious act) arising out of or
% in connection with the use of the software. Without limiting the foregoing, everything in
% the software is provided "as it is" without representation, warranty or condition of any
% kind, either express, implied, or statutory including, but not limited to, implied
% representations, warranties or conditions of merchantabiltiy, fitness for a particular
% purpose, durability, title, non-infringement of intellectual property rights or
% inter-operability of products or services.
%
%
%
%Parameters:    SampFreq        = Sampling Frequency in Hz
%               StartAltitude   = Altitude where the decent begins in m 
%               DescendingRate  = Descending rate during landing in m/s
%               AircraftType    = A340 or VFW614
%               MarcovProbs     = cell array of nxn matrices of transition probabilities for the ground reflection
%               MarkovAtten     = cell array of 1xn matrix of attenuations attatched to the MarcovProbs matrix
%               Altitudevector  = Vector of Altitudes identifying the element of the cell array
%               Elevation       = Elevation of the Satellite [5 ... 90] 

%Natural constants
Co = 2.998e8; %speed of light in m/s


% ------------------------ Internal Parameters -------------------
% please note that the spectral parameters for the reflections is inside the functions at the end of this program

state.Params.CarrierFrequency = 1.57542e9;
state.Params.WaveLength = Co/state.Params.CarrierFrequency;

% ----- Fuselage Filter  ------

state.Fuselage.NrFuselageRealisations = 300;
state.Fuselage.FuselageEchoBW = 2;%Hz
state.Fuselage.Reflected.Power = -14.1; %dB
state.Fuselage.Refracted.Power = -14.1; %dB
state.Fuselage.Reflected.Delay = 1.5e-9; %s
state.Fuselage.Refracted.ConstantPhase=rand(1,1)*2*pi;
state.Fuselage.Reflected.ConstantPhase=rand(1,1)*2*pi;

% ----- Ground Reflection ------

state.Ground.NrGroundRealisations = 300;
state.Ground.GroundEchoBW = 20;%Hz
state.Ground.LastState = 1; % Initialise Start state of the Markov model
state.AntennaHeight = 5;%m
state.Ground.Power = 0; %dB

% ------------------------ security checks -----------------------

% -- SampFreq --
if SampFreq<20
    error ('Sampling Frequency must be >20 Hz')
end %if

state.Params.SampFreq=SampFreq;

% -- StartAltitude --

if (StartAltitude > 1500)| (StartAltitude<100) 
    error('StartAltitude must be between 100 and 1500m')
end %if

state.Params.StartAltitude=StartAltitude;

% -- DescendingRate --

if (DescendingRate>10) | (DescendingRate<0.5)
    error('DescendingRate must be between 0.5 and 10 m/s')
end %if

state.Params.DescendingRate=DescendingRate;

% -- AircraftType --

PossibleAircraftTypes={'A340','VFW614'};

if ~any(strcmp(AircraftType,PossibleAircraftTypes))
    error('Only A340 or VFW614 are allowed Aircraft Types')
end %if

state.Params.AircraftType=AircraftType;

% -- MarcovProbs --

if ~iscell(MarcovProbs)
    error('MarcovProbs must be a cell array')
end %if

hulp=size(MarcovProbs);
if hulp(1)~=1 & hulp(2)>=0 
    error ('Cell array must be at 1xn')
end %if

if hulp(2)~=length(Altitudevector)
    error('The MarcovProbs Cell must have as many elements as the altitude vector')
end%if

for dhv=1:length(MarcovProbs)
    hulp=size(MarcovProbs{dhv});
    if hulp(1)~=hulp(2) | hulp(1)<2
        error ('MarcovProbs must be a square matrix at least 2x2')
    end %if
    if ~all(abs((sum(MarcovProbs{dhv}'))-ones(1,hulp(1)))<eps)
        error (['The sum of the transition probabilities of MarcovProbs must be one! The Line sum is:',num2str(sum(MarcovProbs{dhv}'))])    
    end %if
end %for

state.Params.MarcovProbs=MarcovProbs;

% -- MarkovAtten --

if ~iscell(MarkovAtten)
    error('MarkovAtten must be a cell array')
end %if

dhv=size(MarkovAtten);


if dhv(2)~=length(Altitudevector)
    error('The MarcovAtten Cell must have as many elements as the altitude vector')
end%if

if dhv(1)~=1 | dhv(2) <1    
    error(['MarkovAtten must be a 1xn cell array'])
end %if

for sdhv=1:dhv(2)
    dhv=size(MarkovAtten{sdhv});
    if dhv(1)~=1 | dhv(2)~=hulp(1)
        error(['MarkovAtten must be a 1x',num2str(hulp(1)),' matrix'])
    end %if
end %for 

state.Params.MarkovAtten=MarkovAtten;


% -- AltitudeVector --

if max(Altitudevector)<StartAltitude
    error('The StartAltitude must be covered by the Altitudevector and must therefore be smaller than the largest element of the AltitudeVector')
end %if



if ~isequal(Altitudevector,fliplr(flipud(sort(Altitudevector))))
    error('Altitudevector must be a descending line vector')
end %if

state.Params.AltitudeVector=[Altitudevector(2:end),0];

% -- Elevation --

if Elevation < 10 | Elevation > 70
    error('Elevation bust be between 10 and 70 deg')
end %if

state.Params.Elevation=Elevation/180*pi;

% -- Azimuth --

if Azimuth <10 | Azimuth >=350 | (Azimuth>170 & Azimuth <190)
    error('Azimuth must be between 10...170 or 190 ... 350')
end%if

state.Params.Azimuth=Azimuth/180*pi;

disp ('Security checks passed')


% =============================================
% =========   Initialisations   ===============
% =============================================

% ----- time -----

state.time=0;

% ----- Altitude ------

state.CurrentAltitude=state.Params.StartAltitude;



% ----- Fuselage Filter Refracted ------

f=linspace(-state.Fuselage.FuselageEchoBW,state.Fuselage.FuselageEchoBW,1e5+1);
[y,state.Fuselage.Reflected.ConstantPower]=FuselageFcn(f,state.Params.Elevation,state.Params.Azimuth,state.Params.AircraftType);

state.Fuselage.Refracted.ConstantPower=state.Fuselage.Reflected.ConstantPower;

Y=cumsum(y);
Y=Y/Y(end);

r=rand(1,state.Fuselage.NrFuselageRealisations);

random=r;                                                          % standard implementation

state.Fuselage.Refracted.FreqVector=interp1(Y,f,random,'spline');
state.Fuselage.Refracted.Phases=2*pi*rand(1,state.Fuselage.NrFuselageRealisations);


% ----- Fuselage Filter Reflected ------

r=rand(1,state.Fuselage.NrFuselageRealisations);

random=r;                                                          % standard implementation

state.Fuselage.Reflected.FreqVector=interp1(Y,f,random,'spline');

state.Fuselage.Reflected.Phases=2*pi*rand(1,state.Fuselage.NrFuselageRealisations);

% ----- Ground Reflection ------
clear Y y f;

f=linspace(-state.Ground.GroundEchoBW,state.Ground.GroundEchoBW,1e5+1);
y=GroundFcn(f);

Y=cumsum(y);
Y=Y/Y(end);

r=rand(1,state.Ground.NrGroundRealisations);

random=r;                                                          % standard implementation

state.Ground.FreqVector=interp1(Y,f,random,'spline');
state.Ground.Phases=2*pi*rand(1,state.Ground.NrGroundRealisations);


disp('System initialised')

% Initialize state machine probability vectors

for dhv = 1:length(MarcovProbs)
    sz=size(MarcovProbs{dhv});
    for hhh=1:sz(2)
        l=length(MarcovProbs{dhv}(hhh,:));
        hulp=cumsum(MarcovProbs{dhv}(hhh,:));
        state.Ground.TransitionMatrix.upper{dhv}(hhh,:)=[0,hulp(1:end-1)];
        state.Ground.TransitionMatrix.lower{dhv}(hhh,:)=hulp;

    end %for hhh
end %for

% =================================================================
%                        Private functions
% =================================================================


function [out,ConstantPower] = FuselageFcn(f,Elevation,Azimuth,AircraftType)

switch AircraftType
case 'A340'
    load FuselageSpec_A340.mat;
case 'VFW614'
    load FuselageSpec_ATTAS.mat;
otherwise
    error('This error must never happen')
end

Az = abs(asin(sin(Azimuth))*180/pi);
El = Elevation*180/pi;

C2 = [El^4 El^3 El^2 El 1] * Ab2 * [Az^4 Az^3 Az^2 Az 1]';              %dB
C3 = [El^4 El^3 El^2 El 1] * Ab3 * [Az^4 Az^3 Az^2 Az 1]';              %1/Hz
 
ConstantPower = - [El^4 El^3 El^2 El 1] * Amean * [Az^4 Az^3 Az^2 Az 1]'; %dB

out = 10.^((-(C2+ConstantPower)+C2.*exp(C3.*abs(f)))/20);               %amplitude normalized


function out = GroundFcn(f)

C1 = 2.92;                                                              %deviation in Hz
out = exp(-(f.^2)/(2*C1^2));                                            %amplitude normalized