% Skript interpol_1.m, in dem eine Interpolation mit Faktor
% L untersucht wird.

clear;
% -------- Eingangssignal
fp = 900;       % fmax Signal in Hz 
nellip = 6;     % Ordnung des Filters für das Nutzsignal
dp = 1;         % Welligkeit im Durchlassbereich und 
ds = 60;        % Dämpfung im Sperrbereich  
[b,a] = ellip(nellip, dp, ds, 0.9*2*pi*fp, 's');
[H, w] = freqs(b,a,512);
figure(1);     clf; % ----------------------
subplot(211), semilogx(w/(2*pi), 20*log10(abs(H)));
title(['Amplitudengang des analogen Tiefpassfilters für das'...
    ' Eingangssignal']);
xlabel('Frequenz in Hz');     grid on;    axis tight;
subplot(212), semilogx(w/(2*pi), angle(H));
title('Phasengang des analogen Tiefpassfilters für das Eingangssignal');
xlabel('Frequenz in Hz');     grid on;    axis tight;
axis tight;

% -------- Parameter des Interpolationssystems
L = 4;
fs = 2000;      % Abtastfrequenz
Ts = 1/fs;      % Abtastperiode
fs_prim = fs*L; % Abtastfrequenz des interpolierten Signals
nord = 256;
hint = fir1(nord, 0.9*1/L);     % FIR-Interpolationsfilter
[Ha,w] = freqz(hint,1,512);

figure(2);   clf; % ----------------------
subplot(211), plot(fs_prim*w/(2*pi), 20*log10(abs(Ha)));
grid on;   La = axis;   axis([La(1:2), -100, La(4)]);
title('Amplitudengang des FIR-Interpolationsfilter');
xlabel('Frequenz in Hz');     grid on;    axis tight;
subplot(212), plot(fs_prim*w/(2*pi), unwrap(angle(Ha)));
title('Phasengang des FIR-Interpolationsfilter');
xlabel('Frequenz in Hz');     grid on;    axis tight;


% -------- Aufruf der Simulation
Tsim = 1;
sim('interpol1',[0:Ts/L:Tsim]);   % Feste Schrittweite für 
                                  % die Ergebnisse
y1 = y.Data(:,1);        % Mit Ts abgetastetes Eingangssignal
y2 = y.Data(:,2);        % Ausgang des Aufwärtstaster (Pulse)
y3 = y.Data(:,3);        % Mit Faktor L interpoliertes Signal 
y4 = y.Data(:,4);        % Analoges bandbegrenztes Eingangssignal
t = y.Time;
figure(3);    clf; % ----------------------
subplot(311), plot(t, y4); hold on;
stairs(t, y1);   hold off;
y_max = max(y4);
La = axis;   axis([0.2, 0.22, -y_max,y_max]);    
title('Analoges und mit Ts abgetastetes Eingangssignal');
xlabel('s');    grid on;
subplot(312), stairs(t, y2); hold on;
stairs(t, y1);   hold off;
La = axis;   axis([0.2, 0.22,  -y_max,y_max])
title(['Mit Ts abgetastetes Eingangssignal und Ausgangssignal'...
    ' des Aufwärtstasters']);
xlabel('s');    grid on;
subplot(313), plot(t, y4);  hold on;
stairs(t, y3);   hold off;
La = axis;   axis([0.2, 0.22, -y_max,y_max])
title(['Analoges und interpoliertes Signal']);
xlabel('s');    grid on;

% -------- Spektrale Leistungsdichten
nfft = 256;
[Y2,w] = pwelch((y2+j*eps)*L,hamming(nfft),nfft/2,nfft,fs*L);
[Y3,w] = pwelch(y3+j*eps,hamming(nfft),nfft/2,nfft,fs*L);
figure(4);    clf; % ----------------------
subplot(211);    plot(w, 10*log10(Y2));
title(['Spektrum des Signals nach dem Aufwärtstaster']);
xlabel([' Hz (fs*L = ',num2str(fs*L),' Hz)']);  
ylabel('dBWatt/Hz');    grid on;

subplot(212);    plot(w, 10*log10(Y3));
title(['Spektrum des interpolierten Signals']);
xlabel([' Hz (fs*L = ',num2str(fs*L),' Hz)']);  
ylabel('dBWatt/Hz');    grid on;

% -------- Spektrum des mit Ts abgetastete Signal
y1Ts = y1(1:L:end); % Dezimierung, um die Abtastfrequenz fs zu erhalten
[Y1Ts,w1] = pwelch((y1Ts+j*eps),hamming(nfft),nfft/2,nfft,fs);
[Y1,w2] = pwelch((y1+j*eps),hamming(nfft),nfft/2,nfft,fs*L); % Spektrum
           % mit dem Einfluss des Halteglieds bullter Ordnung
figure(5);   clf; % ----------------------
subplot(211), plot(w1, 10*log10(Y1Ts)); 
title(['Spektrum des mit Ts abgetasteten Signals']);
xlabel([' Hz (fs = ',num2str(fs),' Hz)']);  
ylabel('dBWatt/Hz');    grid on;

subplot(212), plot(w2, 10*log10(Y1));
title(['Spektrum des Signals mit dem Einflus des Halteglieds',...
    ' nullter Ordnung']);
xlabel([' Hz (fs*L = ',num2str(fs*L),' Hz)']);  
ylabel('dBWatt/Hz');    grid on;


% -------- Parseval Theorem
Pzeit = var(y3),       % Varianz
Pfreq = sum(Y3)*fs*L/nfft,  % Leistung aus der Spektralen Leistungsdichte


