% Skript PSD_korrelation_1.m, in dem die PSD
% über die Autokorrelationsfunktion ermittelt wird
% Arbeitet auch mit dem Modell PSD_korrelation1.slx

clear;
% -------- Gefilterte Sequenz
nx = 10000;         % Größe der Sequenz
varianz = 1;        % Varianz des weißen Rauschens
randn('seed', 17937);     % Startwert für die Rauschgenerierung
x = sqrt(varianz)*randn(1,nx);   % Rauschsequenz

h = fir1(256, 0.4); % Impulsantwort (Koeffizienten) des FIR-Filters 
y = filter(h,1,x);  % Gefiltertes Signal
ny = nx;
Ts = 1/1000;      fs = 1000;  % Abtastperiode und Abtastfrequenz

% -------- Autokorrelationsfunktion
max_lag = 50;   % Maximalwert für die Verspätung der Autokorrelation
%[Ryy,lag] = xcorr(y,max_lag, 'biased');  % Autokorrelationsfunktion
     % mit Fensterfunktion Bartlett gewichtet
[Ryy,lag] = xcorr(y,max_lag, 'unbiased');  % Autokorrelationsfunktion
     % ohen Fensterfunktion
nyy = length(Ryy);
Ryyw = Ryy.*bartlett(nyy)';  % Gewichtung der Autokorrelation
%Ryyw = Ryy.*hamming(nyy)';

figure(1);    clf;
subplot(211), plot((0:ny-1)*Ts, x,(0:ny-1)*Ts, y,'r');
title('Unabhängige Zufallssequenz und gefilterte Sequenz');
xlabel('Zeit in s');     grid on;    axis tight;
La = axis;     axis([300*Ts, 500*Ts, -4, 4])

subplot(212), plot(lag, Ryyw);
title('Mit Bartlett-Fenster gewichtete Autokorrelationsfunktion');
xlabel('Verspätung');     grid on;

% -------- DFT der Autokorrelationsfunktion;
nfft = 2^(nextpow2(nyy));    % Anzahl Stützstellen (Bins) der FFT
Syy = abs(fft(Ryyw,nfft));   % Betrag der FFT der Autokorrelation 

figure(2);    clf;
subplot(211), plot((0:nfft-1)/nfft, 10*log10(Syy));
title('Spektrale Leistungsdichte in dBW/Relative-Frequenz');
xlabel('Relative Frequenz f/fs');   grid on;

subplot(212), plot((-nfft/2:(nfft/2-1))/nfft, 10*log10(fftshift(Syy)));
title('Spektrale Leistungsdichte in dBW/Relative-Frequenz');
xlabel('Relative Frequenz f/fs');   grid on;
axis tight;

figure(3);    clf;
subplot(211),plot((0:nfft-1)/(nfft*Ts), 10*log10(Ts*Syy));
title('Spektrale Leistungsdichte in dBW/Hz');
xlabel('Frequenz Hz');   grid on;

subplot(212), plot((-nfft/2:(nfft/2-1))/(nfft*Ts), 10*log10(fftshift(Syy)));
title('Spektrale Leistungsdichte in dBW/Hz');
xlabel('Frequenz Hz');   grid on;
axis tight;

% -------- Parseval-Theorem
Pzeit = sum(Syy)/nfft,      % Leistung aus der spektralen Leistungsdichte
Pfreq = sum(y.^2)/ny,       % Leistung aus der Zeitfunktion

% -------- Untersuchung mit Simulink-Modell PSD_korrelation1.slx
lag_max = 64;
nfft = 2^(nextpow2(2*lag_max));
Tsim = 10;

options = simset('Solver', 'FixedStep');
sim('PSD_korrelation1',[0, Tsim],options);

Rxx_sim = Rxx_sim.signals.values(:,:,end);    % Autokorrelationsfunktion
nxx = length(Rxx_sim);
Sxx = Sxx.signals.values(:,:,end);    % Spektrale Leistungsdichte

figure(4);     clf;
subplot(211),   plot(-(nxx-1)/2:(nxx-1)/2, Rxx_sim);
axis tight
title('Gemittelte Autokorrelationsfunktion aus der Simulink-Simulation');
xlabel('Verspätung');     grid on;
subplot(212), 
plot((-nfft/2:(nfft/2-1))/(nfft*Ts), 10*log10(fftshift(Sxx)));
title('Spektrale Leistungsdichte in dBW/Hz aus der Simulink-Simulation')
xlabel('Frequenz in Hz');     grid on;

