% Skript welch_1.m, in dem die Ermittlung der spektralen
% Leistungsdichte mit Hilfe der DFT untersucht wird
% Arbeitet mit Modell welch1.slx

clear;
% --------- Parameter
fs = 1000;         % Abtastfrequenz
Ts =1/fs;          % Abtastperiode
fd = 200;          % Durchlassfrequenz des FIR-Tiefpassfilters (fd < fs/2)
nfft = 256;        % Anzahl Stützstellen für die FFT
pw = sum(hamming(nfft).^2)/nfft; % Mittlere Leistung der Fensterfunktion

% --------- Aufruf der Simulation
Tsim = 10;         % Dauer der Simulation
options = simset('solver', 'ode45');
sim('welch1',[0, Tsim], options);
t = x.time;                    % Simulationszeit
xf = x.signals.values(:,1);    % Gefiltertes Rauschen 
xu = x.signals.values(:,2);    % Weißes Rauschen Varianz 1

Sxx = Sxx.signals.values(:,:,end);
Sxxw = Sxxw.signals.values(:,:,end);

figure(1);     clf;
subplot(211), stem(t, xu);
hold on, plot(t, xu);   hold off;
La = axis;     axis([500*Ts, 700*Ts, min(xu), max(xu)]);
title('Weißes Rauschen der Varianz 1');
xlabel('Zeit in s');     grid on;

subplot(212), stem(t, xf);
hold on, plot(t, xf);   hold off;
La = axis;     axis([500*Ts, 700*Ts, min(xf), max(xf)]);
title('Gefiltertes weißes Rauschen');
xlabel('Zeit in s');     grid on;

figure(2);     clf;
subplot(211), plot((0:nfft-1)*fs/nfft, 10*log10(Sxx));
title('Spektrale Leistungsdichte in dBW/Hz im Bereich 0 bis fs');
xlabel('Frequenz in Hz');     grid on;      
subplot(212), plot((-nfft/2:nfft/2-1)*fs/nfft, 10*log10(fftshift(Sxx)));
title('Spektrale Leistungsdichte in dBW/Hz im Bereich -fs/2 bis fs/2');
xlabel('Frequenz in Hz');     grid on;      

figure(3);     clf;
subplot(211), plot((0:nfft-1)*fs/nfft, 10*log10(Sxxw));
title(['Spektrale Leistungsdichte in dBW/Hz mit Welch-Block',...
    ' im Bereich 0 bis fs']);
xlabel('Frequenz in Hz');     grid on;      

subplot(212), plot((-nfft/2:nfft/2-1)*fs/nfft, 10*log10(fftshift(Sxxw)));
title(['Spektrale Leistungsdichte in dBW/Hz mit Welch-Block',...
    ' im Bereich -fs/2 bis fs/2']);
xlabel('Frequenz in Hz');     grid on;      
