% ##################################################################################
% ##  Loesung:  Schaetzung des Leistungsdichtespektrums nach Blackman-Tuckey      ##
% ##  --------------------------------------------------------------------------  ##
% ##  Benoetigte(s) m-File(s):  lblack.m und darin: lrader.m                      ##
% ##################################################################################

N  = 2^12;            
N_vor=500;            %Anzahl der Abtastwerte um das System
                      %in den eingeschw. Zustand zu setzen
n  = randn(1,N);      % Modellrauschprozesse
h1 = [1 -4 6 -4 +1];  x1_all = filter(h1,1,n); x1=x1_all(N_vor+1:length(x1_all));
h2 = [1 0.8];         x2_all = filter(1,h2,n); x2=x2_all(N_vor+1:length(x2_all));

% "Wahres" Leistungsdichtespektrum zum Vergleich berechnen
NFFT = 2^12;  
Sx1x1 = abs(fft(h1,NFFT)).^2;                axis_Sx1=[0 1 0 1.35*max(Sx1x1)];
Sx2x2 = ones(1,NFFT)./abs(fft(h2,NFFT)).^2;  axis_Sx2=[0 1 0 1.35*max(Sx2x2)];
fT = 0:1/NFFT:1-1/NFFT;

% Schaetzung des Leistungsdichtespektrums fuer M=8, 32, 128
win = 'triang';
   
Sxx=lblack(x1,   8,win,NFFT);
figure;   plot(fT,Sx1x1,'--'); xlabel('Omega/2pi = f T'); 
hold on;  plot(fT,Sxx); ylabel('Sxx(exp(j*Omega))');
title('Bl.-Tukey-Sch. von Sx1x1(exp(j*Omega)), M=8 (Dreieck)');axis(axis_Sx1);  

Sxx=lblack(x1,   32,win,NFFT); 
figure;   plot(fT,Sx1x1,'--'); xlabel('Omega/2pi = f T');
hold on;  plot(fT,Sxx); ylabel('Sxx(exp(j*Omega))');
title('Bl.-Tukey-Sch. von Sx1x1(exp(j*Omega)), M=32 (Dreieck)'); axis(axis_Sx1); 

Sxx=lblack(x1,   128,win,NFFT); 
figure;   plot(fT,Sx1x1,'--'); xlabel('Omega/2pi = f T');
hold on;  plot(fT,Sxx); ylabel('Sxx(exp(j*Omega))');
title('Bl.-Tukey-Sch. von Sx1x1(exp(j*Omega)), M=128 (Dreieck)'); axis(axis_Sx1); 
input('....Zum Fortfahren: RETURN druecken');

Sxx=lblack(x2,   8,win,NFFT);
figure;   plot(fT,Sx2x2,'--'); xlabel('Omega/2pi = f T'); 
hold on;  plot(fT,Sxx); ylabel('Sxx(exp(j*Omega))');
title('Bl.-Tukey-Sch. von Sx2x2(exp(j*Omega)), M=8 (Dreieck)');axis(axis_Sx2);  

Sxx=lblack(x2,   32,win,NFFT); 
figure;   plot(fT,Sx2x2,'--'); xlabel('Omega/2pi = f T');
hold on;  plot(fT,Sxx); ylabel('Sxx(exp(j*Omega))');
title('Bl.-Tukey-Sch. von Sx2x2(exp(j*Omega)), M=32 (Dreieck)'); axis(axis_Sx2); 

Sxx=lblack(x2,   128,win,NFFT); 
figure;   plot(fT,Sx2x2,'--'); xlabel('Omega/2pi = f T');
hold on;  plot(fT,Sxx); ylabel('Sxx(exp(j*Omega))');
title('Bl.-Tukey-Sch. von Sx1x1(exp(j*Omega)), M=128 (Dreieck)'); axis(axis_Sx2); 
% #####  EOF  #####