end  tauD_MHz = sum(tau_pdf(size(tau_pdf,1),2:3))-tau_pdf(1,2); % Max. theo. Laufzeit-Differenz in us.  tauD_fa = round(tauD_MHz/Ta_MHz);       % Maximale theo. Laufzeit-Differenz in Vielfachen von Ta.%--------------------------------------------------------------------------------------------------%  Step S2:  Leistungsnormierung des Mobilfunkkanals durch Normierung der Impulsantwort g_Ta%--------------------------------------------------------------------------------------------------  g_Ta = g_Ta(:)/sqrt(N_echo);            % Hiermit wird sichergestellt, dass der Mobilfunkkanal    L_g_fa = length(g_Ta);                  % die gleiche mittlere Leistungsverstaerkung hat wie das                                            % Filter, d.h. mean_t{\int |h(tau,t)|^2 dtau} =                                          % \int P(tau) dtau != \int |g(tau)|^2 dtau%--------------------------------------------------------------------------------------------------%  Step S3:  Bestimmung des Verzoegerungs-Leistungsspektrums (Power Delay Spectrum) P(tau) und der%            zeitvarianten Gesamt-Impulsantwort h(tau,t). Die tau-Achse wird dabei im Abtasttakt%            tau=n*Ta aufgeloest, die t-Achse im Raster t=m*Ts, wobei m Werte aus t_fs annimmt.%--------------------------------------------------------------------------------------------------%  PDS_Ta := P(n*Ta)    = g^2(n*Ta) * (\sum_{i=1}^{N_echo} \delta((n-tauE_fa(i))*Ta))%                       = \sum_{i=1}^{N_echo} g^2((n-tauE_fa(i))*Ta)%  h_Ta := h(n*Ta,m*Ts) = \sum_{i=1}^{N_echo} exp(j*phase(i,m*Ts)) g((n-tauE_fa(i))*Ta), %                         wobei phase(i,t) := (omegaD(i)*t+THETA(i))%--------------------------------------------------------------------------------------------------%  Bem: - Das schnellste Echo liegt im Indexbereich (1:L_g_fa) von PDS_Ta bzw. h_Ta;  %         das langsamste bei (1:L_g_fa)+max(tauE_fa)-min(tauE_fa).%       - Die Echozeit tau=0 liegt dort, wo das schnellste Echo sein Betragsmaximum hat (g_mpos)%--------------------------------------------------------------------------------------------------  tauEmin_fa = min(tauE_fa);  g2_Ta = abs(g_Ta).^2;                   % Energiedichte des Sende-*Empfangsfilters (fuer PDS_Ta)  omegaD = 2*pi*fDmax*fDnorm;             % Dopplerkreisfrequenzen in rad/s (Spaltenvektor)  tt = t_fs.'/fs;                         % Absolute Beobachtungszeit in Sekunden (Zeilenvektor)  no_slices = length(tt);                 % Number of time slices t to be calculated (fuer h_Ta).  PDS_Ta = zeros(tauD_fa+L_g_fa,1);       % Power Delay Spectrum im Abtasttakt  h_Ta   = zeros(tauD_fa+L_g_fa,no_slices); % Gesamt-Impulsantwort h(tau,t)  for echo_cnt = 1:N_echo,    pos = tauE_fa(echo_cnt)-tauEmin_fa;    ind_Ta = pos+1:pos+L_g_fa;            % Relevanter Index-Bereich von PDS_Ta und h_Ta                                          % PDS_Ta: Ueberlagerung von durch die Echos gewichteten    PDS_Ta(ind_Ta) = PDS_Ta(ind_Ta) + g2_Ta;    % und verschobenen Energiedichten                                          % h_Ta: Ueberlagerung von durch exp-Faktoren gewichteten                                          %   und verschobenen Sende*Empfangsfilter-Impulsantworten                                          % M. Benthin, Diss., TUHH, Gl. (3.65)    h_Ta(ind_Ta,:) = h_Ta(ind_Ta,:) + g_Ta*exp(j*(omegaD(echo_cnt)*tt+THETA(echo_cnt)));  end;  clear g2_Ta omegaD tt;  if verbose(1) > 1    tauQ_MHz = tauE_fa*Ta_MHz;            % Quantisierte Echo-Laufzeiten in us    mean_delay_MHz = mean(tauQ_MHz);      % Mittlere Echo-Laufzeit in us     del_spread_MHz = std(tauQ_MHz);       % "Delay spread" in us     mmh_fa = mean(mean(h_Ta.'));    EsPDS_fa = sum(PDS_Ta)/tau_M;         % Mittlere Gesamtenergie des zeitkont. PDS fuer alle t    Esh_fa   = sum(abs(h_Ta).^2)/tau_M;   % Gesamtenergie pro Kanalzeitscheibe    disp(sprintf('Info (GSUS.m): params %s', ch_name));            disp('  Quantized echo delay times:');    disp(sprintf('    - range :%6.2fus <= tau_q <= %6.2fus', min(tauQ_MHz), max(tauQ_MHz)));    disp(sprintf('    - mean  :%6.2fus =%6.2fT', mean_delay_MHz, mean_delay_MHz*Ta_MHz));    disp(sprintf('    - stddev:%6.2fus =%6.2fT', del_spread_MHz, del_spread_MHz*Ta_MHz));            disp('  Norm. Doppler frequencies:');    disp(sprintf('    - range :  %g <= |f_d/f_{d,max}| <= %g', min(abs(fDnorm)), max(abs(fDnorm))));            disp('  Impulse response:');    disp(sprintf('    - mean{h(tau,tk)}}:  %s', vec_str(mmh_fa)));    disp(sprintf('    - inst. energy:  En(tk) = \\int |h(tau,tk)|^2 dtau:  %f...%f', min(Esh_fa), max(Esh_fa)));    disp(sprintf('    - mean energy :  mean_tk{En(tk)}                 :  %f', mean(Esh_fa)));    disp(sprintf('    - exp. energy :  E{En(tk)} = \\int P(tau) dtau    :  %f', EsPDS_fa));  end;%--------------------------------------------------------------------------------------------------%  Step S4:  Verschiebung der tau-Achse, so dass tau=0 auf das Maximum des mittleren Verzoegerungs-%            Leistungsspektrums PDS_Ta faellt und anschliessende Bestimmung der outputs PDS_out und%            h_out aus PDS_Ta und h_Ta. Als erste Werte von PDS_out bzw. h_out wird derjenigen Wert%            von PDS_Ta bzw. h_Ta genommen, der tau_del(1) Symbolintervalle links des Maximums von%            PDS_Ta liegt. Dies stellt sicher, dass max(PDS_out) auf dem Symbolraster tau=n*Ts%            liegt. Fuer die Impulsantwort (d.h. max(h_out)) gilt dies im allgemeinen nicht!%--------------------------------------------------------------------------------------------------  [dummy, g_mpos] = max(abs(g_Ta));  L12_g_fs = [(g_mpos-1)/tau_M;           % "krumme" Anzahl der Symbolintervalle VOR  max(g_Ta)              (L_g_fa-g_mpos)/tau_M];     % "krumme" Anzahl der Symbolintervalle NACH max(g_Ta)  tau_axis = [min(abs(tau_axis(1)), 1.0); % force range 0 <= tau_axis(i) <= 1.0              min(abs(tau_axis(2)), 1.0)];  tau_del = fix(tau_axis.*L12_g_fs);      % Der fuer die Output-Variablen beruecksichtigte Bereich                                          % der tau/Ts-Achse faengt tau_del(1) Symbolintervalle VOR                                          % max(PDS_Ta) [dort ist tau=0] an und endet ca. tau_del(2)                                          % Symbolintervalle NACH der max. Echo-Laufzeit-Differenz.  [PDSmax PDSmpos] = max(PDS_Ta);         % Maximum des Verzoegerungs-Leistungsspektrums  ind_min_Ta = PDSmpos - tau_M*tau_del(1);% ist immer > 0, selbst fuer tau_axis(1)=1.0.  ind_max_Ta = PDSmpos + tau_M*fix((tauD_fa+tau_M*tau_del(2))/tau_M);  ind_Ta = ind_min_Ta : ind_max_Ta;       % Indices von PDS_Ta fuer die outputs. Fuer tau_axis=[1 1]:                                          % length(ind_Ta) = tau_M*fix(tauD_fa/tau_M)+L_g_fa.  PDS_out = PDS_Ta(ind_Ta);               % Build function output  h_out   = h_Ta (ind_Ta,:);              % Build function output  tau_out = (ind_Ta.'-PDSmpos)/tau_M;     % Build function output                                          % Zeitlich ausgerichtete Echo-Laufzeiten ermitteln  if tau_M == 1                           % wobei tauShft_fa die Verschiebung fuer tauE_fa ist.    tauShft_fa = PDSmpos-g_mpos;    else                                    % Es erfolgt eine Synchronisation auf das st"arkste Echo.     tauShft_fa = rem(PDSmpos-1,tau_M);    % Der Ursprung der Zeitachse wird dann so festgelegt,    end;                                    % dass er im Anfangssektor des Echoprofils liegt und das                                           % st"arkste Echo genau auf eine Symboltaktzeit f"allt.  tauES_fa = tauE_fa - tauShft_fa;        % Verschiebung der Echo-Laufzeiten, so dass die Laufzeit                                           % null auf das Max. des Verz.-Leistungsspektrums faellt.  tauES_fs = tauES_fa.'/tau_M;            % Verschobene quantisierte Echo-Laufzeiten in "krummen"                                           % Vielfachen von Ts (Zeilenvektor fuer Grafikkommandos)  if nargout > 3                          % Build function output    C207par = [tauES_fs(:), fDnorm, THETA];  end;  TDmin_fs = round(fs/fDmax);             % Minimum Doppler period in multiples of the symbol period                                          % (fuer Bilder)  clear h_Ta PDS_Ta;%--------------------------------------------------------------------------------------------------%  Step S5:  Ggf. Bestimmung der Uebertragungsfunktion H(f,t) = F{h(tau,t)} und der %            Scattering Function SF(tau,fD) [Verzoegerungs-Doppler-Leistungsspektrum]:%            - H_out  := H(f,m*Ts),   wobei  |f|  < fmax_2Ts/2Ts   mit  1<=fmax_2Ts<=tau_M%            - SF_out := SF(n*Ta,fD), wobei  |fD| < max{fD}.%--------------------------------------------------------------------------------------------------  if (nargout > 4) | (verbose(5)==1)                       % #####  Uebertragungsfunktion    NFFT_fs = round(f_axis(2));                            % In jedes Freq.intervall der Breite fs=1/Ts    NFFT = 2.^(nextpow2(NFFT_fs*tau_M));                   % fallen so immer NFFT_fs Stuetzstellen!    H_Ta = fft(h_out, NFFT);                               % NFFT-Punkte DFT in tau-Richtung.    H_Ta = [H_Ta(NFFT/2+1:NFFT,:); H_Ta(1:NFFT/2,:)];      % "fftshift" fuer eine Matrix    f_Ts = tau_M*((0:NFFT-1)/NFFT - 0.5).';                % Zugehoerige norm. Frequenzwerte                                                           % - tau_M/2 <= f*Ts < tau_M/2    fmax_2Ts = min(tau_M, f_axis(1));                      % fmax_2Ts = fmax*2Ts !<= tau_M    NFFT_out = fmax_2Ts*NFFT/tau_M;                        % FFT samples required for the freq.    ind_mitte = NFFT/2+1;                                  % range -fmax <= f < fmax    indices = ind_mitte-ceil(NFFT_out/2):ind_mitte+ceil(NFFT_out/2)-1;    H_out    = H_Ta(indices,:);                            % Build function output    f_Ts_out = f_Ts(indices);                              % Build function output    clear H_Ta f_Ts;  end;  if 0% if ((nargout > 5) | (verbose(6)==1)) & (no_slices>10)    % #####  Scattering Function    NFFT = 2048;    SF = abs(fft(h_out.', NFFT).').^2;                     % NFFT-Punkte DFT in t-Richtung.    SF = [SF(:,NFFT/2+1:NFFT)  SF(:,1:NFFT/2)];            % "fftshift" fuer eine Matrix    fD_Ts = ((0:NFFT-1)/NFFT - 0.5)/(t_fs(2)-t_fs(1)).';    fDmax_Ts = fDmax/fs;                                   % fDmax*Ts; auf fs norm. max. D.-frequenz    NFFT_out = fDmax_Ts*(t_fs(2)-t_fs(1))*NFFT;            % FFT samples required for the frequency    ind_mitte = NFFT/2+1;                                  % range -fmax <= f < fmax    indices = ind_mitte-ceil(NFFT_out/2):ind_mitte+ceil(NFFT_out/2)-1;    SF_out = SF(:,indices);                                % Build function output    fD_Ts_out = fD_Ts(indices);                            % Build function output    SF_out = SF_out/max(max(SF_out));                      % Normierung    fD_TDm_out = fD_Ts_out/fDmax_Ts;                       % Normierung  end;%--------------------------------------------------------------------------------------------------%  Step S6:  Bilder zeichnen (je nach den Werten von verbose und no_slices).%--------------------------------------------------------------------------------------------------%  verbose(2)==1:  Verteilung  Haeufigkeiten der Zufallsvariablen \tau_\nu, f_{d,\nu} und \Theta_\nu%  verbose(3)==1:  P(tau)      Power Delay Spectrum (PDS) inkl. Echo-Laufzeiten (2D plot)%  verbose(4)==1:  |h(tau,t)|  Betrag der Impulsantwort  inkl.   " - "  (3D plot, falls no_slices>1)%  verbose(5)==1:  |H(f,t)|    Betrag der Uebertragungsfunktion (3D plot, falls no_slices>1)%  verbose(6)==1:  SF(tau,fD)  Scattering Function [Verzoegerungs-Doppler-Leistungsspektrum] inkl. %                              Echo-Laufzeiten und Dopplerfrequenzen (falls no_slices>10)  %--------------------------------------------------------------------------------------------------   if verbose(2)==1                                         % #####  Plot random parameters    tauQ_MHz = tauE_fa*Ta_MHz;                             % Quantisierte Echo-Laufzeiten in us    fD = fDmax*fDnorm;                                     % Dopplerfrequenzen in Hertz    th = THETA/(2*pi);                                     % Phasenwinkel in rad/2pi    no_bins = N_echo/10;    figure;    subplot(321);    xLines = ones(2,1)*tauQ_MHz.';                         % Striche fuer die quant. Verzoegerungszeiten    yLines = [0; 1]*ones(size(tauQ_MHz.'));    plot(xLines,yLines); axis([0 tauD_MHz 0 1.2]);    title('Quantized echo delay times in us');     subplot(323);    xLines = ones(2,1)*fD.';                               % Striche fuer die Dopplerfrequenzen    yLines = [0; 1]*ones(size(fD.'));    plot(xLines,yLines); axis([-fDmax fDmax 0 1.2]);    title('Doppler frequencies in Hz');     subplot(325);    xLines = ones(2,1)*th.';                               % Striche fuer die Phasenwinkel    yLines = [0; 1]*ones(size(th.'));    plot(xLines,yLines); axis([0 1 0 1.2]);    title('Phase angles in rad/2pi');     subplot(322);    [h_tmp,x_tmp]=hist(tauQ_MHz, no_bins); bar(x_tmp,h_tmp);% Histogramm der quant. Verzoegerungszeiten    axis([0 tauD_MHz 0 1.1*max(h_tmp)]);    title('Histogram of quantized echo delay times');     subplot(324);    [h_tmp,x_tmp]=hist(fD, no_bins); bar(x_tmp,h_tmp);     % Histogramm der Dopplerfrequenzen    axis([-fDmax fDmax 0 1.1*max(h_tmp)]);    title('Histogram of Doppler frequencies');