%By default, the example selects an ITU-R P.681-11 LMS channel and configures the channel to an urban scenario with 3.8 GHz carrier frequency having a mobile terminal moving at a speed of 2 m/s.
% Create an ITU-R P.681-11 channel or Lutz LMS channel
clear
set(0,'DefaultAxesTitleFontWeight','normal');
chan = p681LMSChannel; %p681LMSChannel;
% Set channel properties
if isa(chan,"p681LMSChannel")
    % For ITU-R P.681 LMS channel
    % Environment type
    environment = input('Environment(0:urban,1:subur.,2:RuralWooded,3:Village,4:Residential)');
    switch environment
        case 0
            chan.Environment="Urban";
        case 1
            chan.Environment="Suburban";
        case 2 
             chan.Environment="RuralWooded";
        case 3
             chan.Environment="Village";
        case 4
             chan.Environment="Residential";
        otherwise
            disp('error');
            return;
    end
    % Carrier frequency (in Hz)
    chan.CarrierFrequency = 2e9;
    % Elevation angle with respect to ground plane (in degrees)
    chan.ElevationAngle = input('Elevation angles (degrees): ');
    % Speed of movement of ground terminal (in m/s)
    chan.MobileSpeed = 2;
    % Direction of movement of ground terminal (in degrees)
    chan.AzimuthOrientation = 0;
else
    % For Lutz LMS channel
    % Rician K-factor (in dB)
    chan.KFactor = 5.5;
    % Lognormal fading parameters (in dB)
    chan.LogNormalFading = [-13.6 3.8];
    % State duration distribution
    chan.StateDurationDistribution = "Exponential";
    % Mean state duration (in seconds)
    chan.MeanStateDuration = [21 24.5];
    % Maximum Doppler shift (in Hz)
    chan.MaximumDopplerShift = 2.8538;
end
% Sampling rate (in Hz)
chan.SampleRate = 400;
%Assign a suitable initial state for the channel.
chan.InitialState = "Good";

%Set the fading technique used to realize the Doppler spectrum. The fading technique is either "Filtered Gaussian noise" or "Sum of sinusoids". When FadingTechnique property is set to "Sum of sinusoids", you can also set the number of sinusoids through NumSinusoids property.
chan.FadingTechnique = "Filtered Gaussian noise"; 

%Initialize random number generator with seed. Vary the seed to obtain different channel realizations. The default value 73 is an arbitrary value.
seed = 73;
chan.RandomStream = "mt19937ar with seed";
chan.Seed = seed;
%Display the properties of the channel.

disp(chan)
%Use get to show all properties
%Channel Model
%Generate the channel for a duration of 100 seconds. Use random samples as input waveform.

% Set random number generator with seed
rng(seed);
% Channel duration (in seconds)
chanDur = 100;
% Random input waveform
numSamples = floor(chan.SampleRate*chanDur)+1;
in = complex(randn(numSamples,1),randn(numSamples,1));
% Pass the input signal through channel
[fadedWave,channelCoefficients,sampleTimes,stateSeries] = step(chan,in);
%Channel Visualization
%Visualize the power profile, the space series, and the state series generated as part of channel modeling.

%Plot the power profile of input waveform and the faded waveform.
figure(1)
plot(sampleTimes,20*log10(abs(in)),sampleTimes,20*log10(abs(fadedWave)))
title(['Power Profile of Waveform for Duration ' num2str(chanDur) ' seconds'])
legend('Input Waveform', 'Faded Waveform')
xlabel('Time (in s)')
ylabel('Power (in dB)')


%Plot the space series to show how the instantaneous power of the channel envelope varies with time.

figure(2)
plot(sampleTimes,20*log10(abs(channelCoefficients)))
title(['Space Series of Channel for Duration ' num2str(chanDur) ' seconds'])
xlabel('Time (in s)')
ylabel('Path Gain (in dB)')


%Plot the state series to show how the channel state varies with time.

figure(3)
plot(sampleTimes,stateSeries)
title(['State Series of Channel for Duration ' num2str(chanDur) ' seconds'])
axis([0 sampleTimes(end) -0.5 1.5])
xlabel('Time (in s)')
ylabel('State')
figure(4)
x=10.^((-30:1:5)/20);
[CDF,PDF]=ecdf2(x,abs(channelCoefficients));
%hold off
plot(20*log10(x),CDF,'r-.');
xlabel('x0 [dB]')
ylabel('P[x < x_0]')
hold on
titlestr=sprintf('environment %s, elevation angle=%d',chan.Environment,chan.ElevationAngle);
title(titlestr)