generateaeronauticalchannel.m

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

function [state,out] = GenerateAeronauticalChannel(state,NrOfSamples)

% 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.

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

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

if NrOfSamples <1
    error('NrOfSamples must be >0')
end%if

if ~isequal(NrOfSamples,round(NrOfSamples))
    error('NrOfSamples must be integer')
end %if

out.TimeVector=state.time+(0:1/state.Params.SampFreq:1/state.Params.SampFreq*NrOfSamples-1/state.Params.SampFreq);
state.time=out.TimeVector(end)+1/state.Params.SampFreq;

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

out.CurrentAltitude=max(state.Params.StartAltitude-out.TimeVector*state.Params.DescendingRate,state.AntennaHeight);
state.CurrentAltitude=out.CurrentAltitude(end);                                      

% ----------- Direct Path -------------------

out.Direct.Delay                    =   zeros(1,length(out.TimeVector));
out.Direct.Amplitude                =   ones(1,length(out.TimeVector));

% ----------- Fuselage refracted ------------


out.Fuselage.Refracted.Delay        =   zeros(1,length(out.TimeVector));
out.Fuselage.Refracted.Amplitude    =   10^((state.Fuselage.Refracted.ConstantPower+state.Fuselage.Refracted.Power)/20)*GenerateFadingSignal(state.Fuselage.Refracted.FreqVector,state.Fuselage.Refracted.Phases,out.TimeVector,0);
                                        
                                        

% ----------- Fuselage reflected ------------

out.Fuselage.Reflected.Delay        =   state.Fuselage.Reflected.Delay*ones(1,length(out.TimeVector));
out.Fuselage.Reflected.Amplitude    =   exp(j*state.Fuselage.Reflected.ConstantPhase)*10^(state.Fuselage.Reflected.Power/20)...
                                        + 10^((state.Fuselage.Reflected.ConstantPower+state.Fuselage.Reflected.Power)/20)*GenerateFadingSignal(state.Fuselage.Reflected.FreqVector,state.Fuselage.Reflected.Phases,out.TimeVector,0);

                                   
% --------- Ground Echo ----------

out.Ground.Delay = (out.CurrentAltitude)*2*sin(state.Params.Elevation)/Co; 

% --------- State Machine --------


for Sample=1:NrOfSamples
        
 %   disp(['Sample:',num2str(Sample)])

    % select Prob Matrix in Dependency from Altitude
    

        
    AltitudeWindow=find((out.CurrentAltitude(Sample)>state.Params.AltitudeVector));
    AltitudeWindow=AltitudeWindow(1);
  
    r=rand(1,1);
    
    NewState=find( ...
        r>=(state.Ground.TransitionMatrix.upper{AltitudeWindow}(state.Ground.LastState,:)) ...
        &...
        r <(state.Ground.TransitionMatrix.lower{AltitudeWindow}(state.Ground.LastState,:)) ...
            );
    out.Ground.state(Sample)=NewState;
    state.Ground.LastState=NewState;
    out.Ground.AltitudeWindow(Sample)=AltitudeWindow;

    out.Ground.AprioriAttenuation(Sample)=state.Params.MarkovAtten{AltitudeWindow}(NewState);
   
end %for

% Amplitude Generation

out.Ground.Amplitude    =   10.^(state.Ground.Power/20)*GenerateFadingSignal(state.Ground.FreqVector,state.Ground.Phases,out.TimeVector,state.Params.DescendingRate/state.Params.WaveLength);
out.Ground.Amplitude    =   10.^(out.Ground.AprioriAttenuation/20).*out.Ground.Amplitude;
    
% =================================================================================
%                               Private Functions
% =================================================================================


function out=GenerateFadingSignal(Freqs,Phases,Time,Doppler)

for dhv=1:length(Freqs)
    signal(dhv,1:length(Time))=exp(j*(2*pi*Time*(Freqs(dhv)+Doppler)+Phases(dhv)));
end % for

out=sum(signal)/sqrt(length(Freqs));