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

clear;
% -------- Eingangssignal
fp = 900;       % fmax Signal in Hz 
naord = 8;      % Ordnung des Filters für das Nutzsignal
dp = 1;         % Welligkeit im Durchlassbereich und 
ds = 60;        % Dämpfung im Sperrbereich  
[b,a] = cheby1(naord, dp, 0.85*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 in zwei Stufen
% mit je einem Interpolationsfaktor 10
L = 10;
fs = 2000;      % Abtastfrequenz des Eingangssignals
Ts = 1/fs;      % Abtastperiode
fs_prim = fs*L; % Abtastfrequenz des interpolierten Signals der 
% ersten Stufe
fs_sec = fs_prim*L; % der zweiten Stufe
nord = 128;     % Ordnung des Filters
hint = fir1(nord, 1/L);     % FIR-Interpolationsfilter
[Ha,w] = freqz(hint,1,512);

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

% -------- Aufruf der Simulation
Tsim = 0.5;
sim('interpol2',[0:Ts/(L*L):Tsim]);   % Feste Schrittweite für 
                                  % die Ergebnisse
y1 = y.Data(:,1);        % Mit Faktor L*L interpoliertes Signal 
y2 = y.Data(:,2);        % Mit Ts abgetastetes Eingangssignal
y3 = y.Data(:,3);        % Analoges bandbegrenztes Eingangssignal
t = y.Time;
figure(3);    clf; % ----------------------
subplot(211), plot(t, y3); hold on;
stairs(t, y2);   hold off;
y_max = max(y3);
La = axis;   axis([0.2, 0.22, -y_max,y_max]);    
title('Analoges und mit Ts abgetastetes Eingangssignal');
xlabel('s');    grid on;
subplot(212), stairs(t, y1); hold on;
plot(t, y3);   hold off;
La = axis;   axis([0.2, 0.22,  -y_max,y_max])
title(['Interpoliertes und zeitkontinuierliches Eingangssignal']);
xlabel('s');    grid on;

% -------- Spektrum des interpolierten Signals
nfft = 1024*4;
y1 = y1(2000:end);            % Ohne Einschwingen
[Y1,w] = pwelch(y1,hamming(nfft),nfft/2,nfft,fs*(L*L));% Spektrale
      % Leistungsdichte nur für das Nyquist Intervall 0 bis fs*L*L/2
figure(4);   clf; % ----------------------
plot(w, 10*log10(Y1));
title(['Spektrale Leistungsdichte des interpolierten Signals (Ausschnitt)']);
xlabel([' Hz (fs*L*L = ',num2str(fs*L*L),' Hz)']);  
ylabel('dBWatt/Hz');    grid on;
La = axis;   axis([0, 2000, -100, La(4)]);

% --------- Parseval Theorem
Pzeit = var(y1), 
Pfreq = sum(Y1)*fs*(L*L)/nfft,   
