% ##################################################################################
% ##  Funktion: lywex.m;   AR-(p)-Approximation eines ARMA-Modelles nach Yule-Walker
% ##################################################################################
%
% function [ARcoeff] = lywex(MA,AR,p[,NFFT])
%
% Berechnung der Koeffizienten eines AR-Modelles aus den exakten
% Autokorrelationskoeffizienten
% Input-Argumente:                       
%   MA,AR    :   Parameter des ARMA-Modells
%   p        :   Grad des zu bestimmenden AR-Modells
%   [NFFT]   :   FFT-Laenge
% Output-Argument:
%   ARcoeff  :   Koeffizienten des Yule-Walker-AR-Modells (Zeilenvektor)
%   Sxx_ar   :   geschaetztes AR-Betragspektrum
%   Sxx_arma :   wahres LDS des ARMA-Modells

function [Sxx_ar, Sxx_arma, ARcoeff] = lywex(MA,AR,p,NFFT)

if nargin<4, NFFT = 2^10; end;
delta    = [1 zeros(1,NFFT-1)];
h_arma   = filter(MA,AR,delta);
H_arma   = fft(h_arma);
Sxx_arma = abs(H_arma).^2;           % LDS des ARMA-Modells
                                     % Wahre AKF des ARMA-Modells:
r0_arma = h_arma*h_arma';            % rxx(0)
for ii = 1:p                         % rxx(1)..rxx(p)
  r_arma(ii) = [zeros(1,ii) h_arma]*[h_arma zeros(1,ii)]';
end

R_arma  = toeplitz([r0_arma r_arma(1:length(r_arma)-1)]); 
ARcoeff = -inv(R_arma)*r_arma';
ARcoeff = [1 ARcoeff.'];
h_ar = filter(1, ARcoeff, delta);
H_ar = fft(h_ar);
if p == 0
  sigmak2 = 1;
else
  sigmak2 = R_arma(1,1)-r_arma*inv(R_arma)*r_arma';
end
Sxx_ar = sigmak2 * abs(H_ar).^2;     %geschaetztes AR-Betragspektrum
% #####  EOF  #####