function imp=channel_model;
% function imp=channel_model;
%
%  This function provides impulse responses based on the COST 273
%  MIMO channel model specifications. The impulse responses are
%  given as 4-way array:
%       [Snapshots x Imp_resp x BS_Antennas x MS_Antennas]
%           Snapshots .... the snapshots of the Imp_resp in time
%           Imp_resp ..... the delay index
%           BS_Antennas .. the number of BS antennas
%           MS_Antennas .. the number of MS antennas
%  A filtered version of the impulse response can be obtained by
%  uncommenting the last line of the code.
% The structured field "imp" contains the total impulse response, a
% filtered version of the impulse respons and impulse responses only
% covering the LOS component and the MPC in seperate matrizes.
%  Parameters marked with '**' are not supported up to now.
%
%  Copyright (C)2005  Helmut Hofstetter
% 
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
% 
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
% 
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
% 
%----------------------------------------------------------------------%
%                                                                      %
%              Author: Helmut Hofstetter <hofstetter@ftw.at>           %
%             Created: 25/05/2005 by HH                                %
%   Last modification: 08/07/2006 by HH                                %
%        Organization: Telecommunications Research Center Vienna,      %
%                      Tech Gate Vienna,                               %
%                      Donau-City Str. 1 / 3rd floor,                  %
%                      A-1220 Vienna, Austria.                         %
%             Project: FTW C9 "MIMO UMTS for future packet services"   %
%    Deliverable Type: This M-file is public available                 %
%                 URL: http://www.ftw.at/cost273/                      %
%         Source Code: MATLAB v7.0 SP2                                 %
%             History:                                                 %
%             Created: 1.00, 25/05/2005 by HH (ftw)                    %
%             Created: 1.10, 08/07/2006 by HH (Eurecom)                %
%                                                                      %
%----------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% External Parameter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

cell_type='macro';
scenario='large_urban_macro';

c0=3e8;             % speed of light
f_c=2e9;            % center frequency
h_BS=15;            % BS height
h_MS=1.5;           % MS height
r_BS=[0 0 h_BS];    % Position of BS
r_MS=[200 200 h_MS];    % Position of MS
MSAnr=4;            % Number of Antennas at MS
BSAnr=4;            % Number of Antennas at BS
d_Ant_BS=0.15*0.5;       % Antenna spacing at BS in [Lambda]
d_Ant_MS=0.15*0.5;       % Antenna spacing at MS in [Lambda]
r_Ant_BS=[];        % Antenna position at BS    %**
r_Ant_MS=[];        % Antenna position at MS    %**
theta_Ant_BS=0;     % Antenna orientation at BS    %**
theta_Ant_MS=0;     % Antenna orientation at MS    %**
P_L=[];             %**
h_rooftop=[];       % COST_Hata parameter: rooftop height    %**
w_road=[];          % COST_Hata parameter: width of street    %**
theta_road=[];      % COST_Hata parameter: road orientation    %**
No_of_imp_resp=10;  % Number of impulse responses to be calculated

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute BS-MS distance and related angles
[Theta_BS_MS,Phi_BS_MS,d_BS_MS]=cart2sph(r_MS(1)-r_BS(1),r_MS(2)-r_BS(2),r_MS(3)-r_BS(3));
tau_0=d_BS_MS/c0; % delay of the LOS

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Stochastic Parmeters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%% Visibility region %%%%%
R_C=100;                % not used upto now
L_C=20;                 % it is assumed that all the clusters are visible at the moment

%%%%% Cluster movement %%%%% 
mu_Cv=0;                % [m/s]  % not used up to now
sigma_Cv=0;             % [dB]

%%%%% Cluster power %%%%% 
k_tau=1;                % attenuation coefficient for clusters [dB/祍] (Eq. 4.8)
tau_B=10e-6;            % cut-off delay for attenuation (Eq. 4.8)

%%%%%  Line of sight %%%%% 
d_inf=500;              % [m] %** % not used up to now
R_L=30;                 % [m] %**
L_L=20;                 % [m] %**
EPL=0;                  % definition for the first moment
mu_K=(26-EPL)/6;        % [linear] %**
sigma_K=1;              % [dB] %**

% Selection parameter
K_sel=0;        % choose between the twin clusters and single interacting clusters.
                % K_sel=0 enables single interacting clusters 
                % macro cells: K_sel=0;
                % micro cells: K_sel=0.5;
                % pico cells:  K_sel=1;
                % each float between 0 and 1 is possible

N_C_local=1;            % if N_local==1 only a local cluster around the MS exists
                        % if N_local==2 also a local cluster around the BS exists
                        % there is no possible scenario with just a local
                        % cluster around the BS but no local cluster around
                        % the MS
N_C_add_mean=0.18;      % mean number of add. clusters (twin-clusters + single-interaction clusters)
N_C_add=poissrnd(N_C_add_mean);
if N_C_add==0
    N_C_add=1;          % at least one add. cluster is required
end
N_C_twin=ceil(K_sel*N_C_add);   % twin clusters are slightly preferd:-)
N_C_single=N_C_add-N_C_twin;
N_C=N_C_local+N_C_add;  % total number of clusters
    
% single interacting clusters
r_min=50;              % [m]   minimum distance of a single interacting cluster from the BS
sigma_r=50;            % [m]   variance of the distance BS-MS
sigma_Theta_C_BS=60/180*pi;        % [rad] angular distribution of single interacting clusters

N_MPC=20;               % Number of MPCs per cluster
K_MPC=0;                % Rice factor of clusters   
 
% Diffuse radiation    %**
mu_diff=0.05;           %Std of diffuse radiation
sigma_diff=3.4;         %Var of diffuse radiation

%%%%% Delay spread %%%%% 
mu_tau=0.4e-6;          % [祍]  mean value
sigma_tau=3;            % [dB]  std

%%%%% Angular spread %%%%% 
mu_theta_BS=0.81/180*pi;    % [rad] angular spread in azimuth, BS side, mean
sigma_theta_BS=0.34;        % [dB]  angular spread in azimuth, BS side, std
mu_phi_BS=0.5/180*pi;       % [rad] angular spread in elevation, BS side, mean
sigma_phi_BS=3;             % [dB]  angular spread in elevation, BS side, std
mu_theta_MS=35/180*pi;      % [rad] angular spread in azimuth, MS side, mean
sigma_theta_MS=0;           % [dB]  angular spread in azimuth, MS side, std
mu_phi_MS=0.81/180*pi;      % [rad] angular spread in elevation, MS side, mean
                            % the angular spread in elevation is defined according 
                            % to the COST259 recommmendations and differs from the rest:
                            % pdf(mu_phi_MS)=[0,45癩
sigma_phi_MS=0.34;          % [dB]  angular spread in elevation, MS side, std

%%%%% Shadow fading %%%%% 
sigma_S=6;                  % [dB]

%%%%% Autocorrelation distances %%%%% 
L_S=100;                    % [m] %** not yet used
L_tau=100;                  % [m] %**
L_thetaBS=100;              % [m] %**
L_phiBS=100;                % [m] %**
L_thetaMS=100;              % [m] %**
L_phiMS=100;                % [m] %**

%%%%% Cross correlations %%%%% 
rho=[1 -0.6 -0.6 0   0   0;...        % shadow fading
    -0.6 1   0.5 0   0   0;...        % delay spread
    -0.6 0.5 1   0   0   0;...        % theta BS
     0   0   0   1   0   0;...        % phi BS
     0   0   0   0   1   0;...        % theta MS
     0   0   0   0   0   1];          % phi MS       

corr_mat=chol(rho);    % the correlation matrix is caculated using the Cholesky factorization
 
%%%%% Polarisation %%%%% % not yet used
mu_XPD=-6;              %** 
sigma_XPD=2;            %**
mu_VVHH=0;              %**
sigma_VVHH=-inf;        %**
mu_VHHV=0;              %**
sigma_VHHV=-inf;        %**

[ms_vx,ms_vy]=pol2cart(rand(1)*2*pi,50);
v_MS=[ms_vx ms_vy 0];         % velocity of MS [m/s]

%%%%% Sampling interval for impulse response %%%%%
t_sample=1/(3.84e6);          % [sec] delay resolution of the impulse response
delta_t=500*t_sample;   % [sec] time between two impulse responses
t=0;                    % [sec] absoulte time for simulation

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Filter coefficients
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% channel is 4 times oversampled              
	cosf.Delay = 6; cosf.R = .2; cosf.Fs = 4; cosf.Fd = 1;   
	% init cosine filter 
	cosf.PropD = cosf.Delay * cosf.Fd;
	cosf.ys = rcosine(cosf.Fd, cosf.Fs, 'fir', cosf.R, cosf.Delay);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% angular spreads of all the clusters
corr_randn=randn(N_C,6)*corr_mat;       % compute correlated random variables
theta_C_BS=mu_theta_BS*10.^(0.05*sigma_theta_BS*corr_randn(:,3));
theta_C_MS=mu_theta_MS*10.^(0.05*sigma_theta_MS*corr_randn(:,5));
phi_C_BS=mu_phi_BS*10.^(0.05*sigma_phi_BS*corr_randn(:,4));
phi_C_MS=mu_phi_MS*10.^(0.05*sigma_phi_MS*corr_randn(:,6));

if findstr(scenario,'large_urban_macro')        % there is a special case for the 
                                                % anguler spread definition at the BS 
                                                % for the large urban macro cell scenario
                                                % for the local MS cluster
    phi_MS(1)=rand(1,1)*45/180*pi;     % since there is no cross correlation given for phi_MS this case is omitted.    
end % if findstr(scenario,'large_urban_macro')

% delay spread of the clusters
tau_C=mu_tau.*10.^(0.05*sigma_tau*corr_randn(:,2));    % holds true for all types of clusters

% precompute the cluster size; it is the size a cluster needs to result in
% the required delay spread.
d_tau=0.5*tau_C*c0; % radius of the clusters

% init variable
d_C_MS=zeros(N_C,1);    % distance of cluster from MS
d_C_BS=zeros(N_C,1);    % distance of cluster from BS
a_C_MS=zeros(N_C,1);    % main axis of the cluster ellipsoid in azimuth as seen from the MS
a_C_BS=zeros(N_C,1);    % main axis of the cluster ellipsoid in azimuth as seen from the BS
b_C_MS=zeros(N_C,1);    % axis of the cluster ellipsoid in the delay direction as seen from the MS
b_C_BS=zeros(N_C,1);    % axis of the cluster ellipsoid in the delay direction as seen from the BS
h_C_MS=zeros(N_C,1);    % height of the cluster ellipsoid as seen from the MS
h_C_BS=zeros(N_C,1);    % height of the cluster ellipsoid as seen from the BS
r_C_BS=zeros(N_C,3);    % 3-D positioning vector of BS
r_C_MS=zeros(N_C,3);    % 3-D positioning vector of MS
Phi_C_MS=zeros(N_C,1);    % direction of the cluster in azimuth as seen from the MS
Phi_C_BS=zeros(N_C,1);    % direction of the cluster in azimuth as seen from the BS
Theta_C_MS=zeros(N_C,1);      % direction of the cluster in elevation as seen from the MS
Theta_C_BS=zeros(N_C,1);      % direction of the cluster in elevation as seen from the BS
impn=zeros(No_of_imp_resp,40,MSAnr,BSAnr);          % initialize
imp_MPC=zeros(No_of_imp_resp,40,MSAnr,BSAnr);         % initialize
imp_LOS=zeros(No_of_imp_resp,40,MSAnr,BSAnr);          % initialize
imp_cosine=zeros(No_of_imp_resp,60,MSAnr,BSAnr);  % initialize

% The cluster directions are needed for rotating the cluster and nothing
% else. Therefore it only denotes the actual cluster direction for the twin
% clusters. For single interacting clusters and the local clusters it
% denotes the angle the cluster has to be rotated. This allows for a common
% processing of all the clusters for all further steps.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%