gwssus.m

这是关于移动无线信道COS207模型的MATLAB源程序

  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');