% Skript lattice_11.m, in dem eine Lattice-Struktur
% für ein 'prediction' FIR-Filter untersucht wird

clear;
% ------- FIR-Filter
nord = 32;
hn = fir1(nord, 0.4);   % FIR-Filter zur Generierung eines Signal
    % mit einer Bekannten Autokorrelationsfunktion
rxx = conv(hn, hn);    % Autokorrelation am Ausgang des Filters
    % das mit unabhängiger Sequenz angeregt wird
figure(1);   clf;
subplot(211), stem(0:nord, hn);
title(['Einheitspulsantwort des FIR-Filters zur',...
             'Erzeugung der korrelierten Sequenz']);
xlabel('n');    grid on;   axis tight;       
subplot(212), stem(-nord:nord, rxx);
title('Autokorrelationsfunktion der korrelierten Sequenz');
xlabel('n');    grid on;   axis tight;         

% ------- Normale Gleichung zur Identifikation 
    % des Prediction-Filters    
RXX = toeplitz(rxx(33:end));        % Matrix der normalen Gleichung
[a, e, kl] = levinson(rxx(33:end)); % Koeffizienten des voraussage 
    % FIR-Filters für x[n]
b = fliplr(a);                     % Koeffizienten des rückwärts 
    % voraussage FIR-Filters für x[n-M]
figure(2);     clf;
na = nord;
subplot(311), stem(0:na, a);
title('Koeffizienten des Voraussage-Filters für x[n] inklussive a(0) = 1')
xlabel('k');   grid on;
subplot(312), stem(0:na, b);
title(['Koeffizienten des rückwärts Voraussage-Filters für '...
    'x[n-M] inklussive b(M) = 1'])
xlabel('k');   grid on;

k = tf2latc(a);   % Berechnung der Koeffizienten der Gitter-Struktur
subplot(313), stem(0:na-1, k);
title('Reflection-Koeffizienten der Lattice-Struktur');
xlabel('k');   grid on;
% ------- Filterung eines unabhängigen Zufallsignals der Varianz 1
nx = 200;
randn('seed',157893);
xn = randn(1,nx);
x = filter(hn,1,xn);  % Gefilterte Sequenz mit Autokorrelation Rx
var_x = var(x),
% ------- Prediction und revers prediction Filterung
f = filter(a,1,x);               % Voraussagefehler x[n]
var_f = var(f),                  % Varianz des Schätzungsfehlers
g = filter(b,1,x);               % Voraussagefehler x[n-M]
ypred = filter(-a(2:end),1,x);   % Vorausgesagte Sequenz x[n]
yrev_pred = filter(-b(1:end-1),1,x); % Vorausgesagte Sequenz x[n-M]

figure(3);   clf;
subplot(311), plot(0:nx-1, x);
title('Zufallssequenz x mit Autokorrelation rxx')
xlabel('n');  grid on;
subplot(312), plot(0:nx-1, [0,ypred(1:end-1)]); % Mit einem Takt verspätet
title('Vorausgesagte Sequenz x[n] aus ypred (verspätet mit einem Schritt)')
xlabel('n');  grid on;
subplot(313), plot(0:nx-1, f, 0:nx-1, x-[0,ypred(1:end-1)]);
title('Fehler der Voraussage von x[n]');
xlabel('n');  grid on;

figure(4); 
subplot(311), plot(0:nx-1, [zeros(1,nord),x(1:end-nord)]);
title('Zufallssequenz x mit Autokorrelation rxx (verspätet mit M Schritte)')
xlabel('n');  grid on;
subplot(312), plot(0:nx-1, yrev_pred(1:end)); 
title('Rückwärts vorausgesagte Sequenz x[n-M] aus yrev-pred')
xlabel('n');  grid on;
subplot(313), plot(0:nx-1, g);
title('Fehler der rückwärts Voraussage von x[n-M]');
xlabel('n');  grid on;

% ------- Die direkte Filterung mit den Koeffizienten K
[fd,gd] = latcfilt(k,x);
figure(5);
subplot(211), plot(0:nx-1, fd)
title('Fehler der Voraussage von x[n]');
xlabel('n');   grid on;
subplot(212), plot(0:nx-1, gd)
title('Fehler der Voraussage von x[n-M]');
xlabel('n');   grid on;