% Retrieve the power delay spectrum P(tau), the impulse response h(tau,t)% and the transfer function H(f,t) of a sample COST-207 mobile channel  % using the "Gaussian Stationary Uncorrelated Scattering" (GSUS) % stochastic channel model by Schulze/Hoeher.%% function [PDS_out, h_out, tau_out, C207par, H_out, f_Ts_out] = ...%   GSUS (mPDS, N_echo, fs, fDmax, tau_M, g_Ta, ...%         [t_fs, [tau_axis, [f_axis [, verbose[, ch_name[, rseed]]]]]])% ----------------------------------------------------------------------% INPUT:% ----------------------------------------------------------------------% mPDS     :  String designating the desired transmission environment,%             i.e. the physical channel, according to COST-207:%             - 'TU':  "Typical Urban"  (non-hilly urban environment),%             - 'RA':  "Rural Area"     (non-hilly rural environment),%             - 'BU':  "Bad Urban"      (hilly urban environment),%             - 'HT':  "Hilly Terrain"  (hilly rural environment).%      Note:  These environments correspond to different mean power %             delay spectra (PDS, "Verzoegerungsleistungsspektren").%      N.B.:  To obtain non-standardized environments, see "COST207.m".% N_echo   :  Desired number of echo paths.  Example: N_echo=100.% fs       :  Symbol frequency in baud (=> Symbol period:  Ts=1/fs).%             Ex.: fs=13e6/48=270833 (GSM);  fs=1228800 (QualComm-CDMA).% fDmax    :  Maximum Doppler frequency  fDmax := v*f0/c0, where%             - v       is the speed in m/s of the mobile unit%             - f0      is the carrier frequency in Hertz%             - c0=3e8  is the speed of light in m/s in empty space.% tau_M    :  Oversampling factor determining the sampling period on the %             tau axis:  Ta = Ts/tau_M, where Ts is the symbol period.%      N.B.:  The continuous echo delay times are also quantized to fit%             a Ta-spaced delay time grid.% g_Ta     :  Impulse response g(tau) sampled at tau=n*Ta of a band-%             limiting filter. This may either be the receive filter or%             -- for linear modulation schemes -- the overall filter com-  %             posed of transmit and receive filter (refer to "get_gTR.m").%      N.B.:  For linear modulation, fDmax << 1/L_gTx must be satisfied,  %             where L_gTx is the transmit filter's length, because other-%             wise, g_Ta is not simply the convolution of the transmit %             and receive filter impulse responses.% [t_fs]   :  Vector giving the time instants t in multiples of Ts, %             where the impulse response h(tau,t) will be evaluated.%             To retrieve a single slice of the impulse response, set %             t_fs to a scalar value. To cover a t-range of D Doppler %             periods 1/fDmax, choose e.g. t_fs=((0:D/39:D)/fDmax)*fs.%             Default value is:  t_fs=((0:1/40:1-1/40)/fDmax)*fs.% [tau_axis]: The tau output range is determined by these two relative   %             factors used to suppress outputs close to zero. The 1st %             [2nd] factor gives a percentage of how many symbol periods %             are to be appended to the tau axis BEFORE min{tau} [AFTER %             max{tau}]. Thus, setting tau_axis=[0 0] yields very tight %             bounds, while the max. tau range is obtained by [1.0 1.0].%             Default value is tau_axis=[0.3 0.2].% [f_axis]:   Optional two-element vector determining the range and re- %             solution of the frequency axis for the transfer function:%             - f_axis(1): norm. freq. (must be <= tau_M) used to deter-%               mine the frequency range: |f|<= f_axis(1)/(2Ts). The va-%               lue 1.0 [2.0] corresponds to the Nyquist [symbol] freq.%             - f_axis(2) is the desired number of FFT points in each %               frequency interval with length 1/Ts.%             Default value is f_axis = [min(1.5,tau_M), 128].% [verbose]:  Optional vector determining the kind of info given:%             verbose(1)==1:  print warnings only (default)%                       ==2:  print warnings and info %             verbose(2)==1:  draw random parameters%             verbose(3)==1:  draw power delay spectrum P(tau)%             verbose(4)==1:  draw magn. impulse response |h(tau,t)|%             verbose(5)==1:  draw magn. transfer function |H(f,t)|,%                             where  H(f,t) := F{h(tau,t)} w.r.t. tau% [ch_name]:  Optional string for the figure titles designating the ch.:%             Default value is 'COST207-<mPDS> channel' (if isstr(mPDS)).% [rseed]  :  If this parameter is given, rseed is used to set the %             seeds for the random numbers generator. In this way,%             the mobile channel's parameters can be reproduced.% ----------------------------------------------------------------------% OUTPUT:% ----------------------------------------------------------------------% PDS_out  :  Length N column vector containing a sample power delay  %             spectrum (PDS) P(tau) for tau = n*Ta. To obtain the symbol %             rate sampled PDS, use PDS_out(1:tau_M:length(PDS_out)).%      Note:  N = L*tau_M+1, where L is an integer depending on tau_axis.%      N.B.:  If tau_axis==[1.0 1.0], the following is exactly true:%             sum(PDS_out)/tau_M = sum(abs(g_Ta).^2)/tau_M, i.e. the ex-%             pected mean power amplification of the GSUS channel is one, %             if and only if the power amplification of g(tau) is unity.% h_out    :  (N x length(t_fs)) matrix with a sample impulse response%             h(tau,t) of the selected time-variant COST207 channel.%             Each column contains one slice h(tau,t0) where the tau %             axis is sampled at tau=n*Ta=n*Ts/tau_M and t0 is fixed. %             Conversely, each row represents one channel coefficient  %             h(tau0,t) sampled at t=t_fs*Ts for a fixed value tau0.%      Note:  Use h_out(1:tau_M:size(h_out,1),:) to obtain the%             corresponding symbol rate impulse response.%      N.B.:  For tau_axis==[1.0 1.0], the GSUS channel has an expected%             mean power amplification equal to sum(abs(g_Ta).^2)/tau_M.%             However, the actual mean power amplification of h_out%             will differ from this value as long as length(t_fs)~=inf.% tau_out  :  Length N column vector with the normalized time delay axis  %             tau_out=tau*fa=tau/Ta for use with graphics commands.% C207par  :  "COST207.m" outputs. The 1st column contains the quantized %             and shifted echo delay times in multiples of Ts, the 2nd col.%             the norm. Doppler frequencies and the 3rd col. the phases.% H_out    :  (NFFT x length(t_fs)) matrix with a sample transfer func-%             tion H(f,t) of the selected time-variant COST207 channel.%             Each column contains one slice H(f,t0) of the transfer%             function for a fixed value of t0. Conversely, each row re-%             presents the rate of change of each coefficient H(f0,t).% f_Ts_out :  Length NFFT column vector with the normalized freq. axis  %             -f_axis(1)/2Ts<=f<f_axis(1)/2Ts for use with graphics cmds.% ----------------------------------------------------------------------% PREREQUISITES:% ----------------------------------------------------------------------% Required functions:  "COST207.m"% ----------------------------------------------------------------------% REFERENCES:% ----------------------------------------------------------------------% - H. Schulze: "Stochastische Modelle und digitale Simulation von Mobil-%   funkkanaelen", Kleinheubacher Berichte, Vol.32, pp 473-483, 1989.% - P. Hoeher: "A Statistical Discrete-Time Model for the WSSUS Multipath%   Channel", IEEE Trans. on Veh. Technol, Vol.41, no.4, pp 461-468, 1992.% ----------------------------------------------------------------------% AUTHOR:  Marcus Benthin, 01.07.92  (m-files: "tap207*.m", "eng207*.m")%          Dieter Boss,    17.03.97% ----------------------------------------------------------------------% <-------------------------------  max. linewidth for atops  ------------------------------------->|function [PDS_out, h_out, tau_out, C207par, H_out, f_Ts_out] = ...  GSUS (mPDS, N_echo, fs, fDmax, tau_M, g_Ta, t_fs, tau_axis, f_axis, verbose, ch_name, rseed)  if nargin <  7,  t_fs = ((0:1/40:1-1/40)/fDmax)*fs; end;  % 40 times slices within one min.   if nargin <  8,  tau_axis = [0.3, 0.2]; end;              % Doppler period  if nargin <  9,  f_axis = [1.5, 128];   end;  if nargin < 10,  verbose = [1, 0, 0, 0, 0];  end;         % print warnings only  of_C207 = 'of a sample COST207';  if isstr(mPDS),  of_C207 = sprintf('%s-%s', of_C207, mPDS); end;  if nargin < 11,  ch_name = sprintf('%s channel', of_C207);  else,            ch_name = sprintf('%s  (%s)  channel', of_C207, ch_name); end;  t_fs = t_fs(:);                                           % force column vector  if length(tau_axis)==1, tau_axis(2)=tau_axis(1); end;  if length(f_axis)==1,   f_axis(2)=128; end;%--------------------------------------------------------------------------------------------------%  Step S0:  Input-Argumente N_echo, tau_M, g_Ta, t_fs ueberpruefen und ggf. Warnings ausgeben%--------------------------------------------------------------------------------------------------  if verbose(1) > 0    if (N_echo < 50)      disp(sprintf('Warning (GSUS.m):  N_echo=%d should be >= 50 to assure a good', N_echo));              disp('                   approximation of a GAUSSIAN stationary');              disp('                   uncorrelated scattering model.');    end;    if (tau_M < 1) | (tau_M > 64)      disp(sprintf('Warning (GSUS.m):  tau_M=Ts/Ta=%d should be in the range  1..64.', tau_M));    end;    energie_g = sum(abs(g_Ta).^2)/tau_M;    if abs(1-energie_g) > 1e-10              disp('Warning (GSUS.m):  h_out will NOT have an expected mean power ampli-');              disp('                   fication of one, because the energy of g_Ta is');      disp(sprintf('                   not one, but \\int |g(tau)|^2 dtau = %f.', energie_g));    end;    if length(t_fs) > 1      sl_per_TDmin = round((fs/fDmax)/((max(t_fs)-min(t_fs)+1)/(length(t_fs)-1)));      if sl_per_TDmin < 20        disp(sprintf('Warning (GSUS.m):  Calculating only %d<20 impulse response slices per', sl_per_TDmin));        disp(sprintf('                   min. Doppler period 1/fDmax = %dT may lead to', round(fs/fDmax)));                disp('                   an insufficient approx. of time variant behavior.');      end;    end;  end%--------------------------------------------------------------------------------------------------%  Step S1:  Auswuerfelung der Mobilfunkparameter (Rayleigh-Anteile):%              (1) Echo-Laufzeiten in Mikrosek.:   tauE_MHz := tauE*1MHz = tauE/1us%              (2) Normierte Dopplerfrequenzen :   fDnorm = fD/max{fD}  =>  |fDnorm| < 1 %              (3) Startphasen                 :   0 < THETA < 2pi. %            Danach werden die Echo-Laufzeiten auf das Abtastraster n*Ta=n*Ts/tau_M quantisiert %            [tauE_fa:=round(tauE/Ta)].%--------------------------------------------------------------------------------------------------  if nargin > 11    [tauE_MHz, fDnorm, THETA, tau_pdf] = COST207 (mPDS, N_echo, rseed);  else     [tauE_MHz, fDnorm, THETA, tau_pdf] = COST207 (mPDS, N_echo);  end;  N_echo = length(tauE_MHz);              % N_echo kann sich um +/-1 geaendert h. => Neubestimmung  Ts_MHz = 1e6/fs;                        % Symbolperiode in Mikrosekunden.  Ta_MHz = Ts_MHz/tau_M;                  % Abtastperiode in Mikrosekunden.  tauE_fa = round(tauE_MHz/Ta_MHz);       % Echo-Laufzeiten in ganzzahligen Vielfachen des                                           % Abtasttaktes (koennen auch null sein).  tauDcl_MHz = min(tau_pdf(:,3));         % Dauer des kuerzesten Echo-Clusters (z.B. 2us bei 'HT').  no_tauEcl_fa = ceil(tauDcl_MHz/Ta_MHz); % Anzahl der quantisierten Echo-Laufzeiten, die in das                                          % kuerzeste Echo-Cluster fallen. Annahme: das Ta-Raster  if (no_tauEcl_fa < 4) & (verbose(1) > 0)% "trifft" die Anfangszeit dieses Clusters.            disp('WARNING (GSUS.m):  The approximation of the tau-continuous channel im-');            disp('                   pulse response may be inexact, since there is/are');    disp(sprintf('                   only %d<4 possible value(s) for the echo delay times', no_tauEcl_fa));    disp(sprintf('                   within the shortest echo cluster (duration: %.1fus).', tauDcl_MHz));            disp('                   Increase tau_M or fs for a better approximation.');