% Skript hilbert_filter_2.m, in dem ein Hilbert-Filter
% untersucht wird. Arbeitet mit Modell hilbert_filter2.slx

clear;
% -------- Parameter des Systems
fs = 1000;     Ts = 1/fs;   % Abtastfrequenz und Periode
% -------- Entwicklung des Hilbert-Filters
nord = 64;           % Gerade Ordnung
Tw = 0.2;            % 'Transition' Band
d = fdesign.hilbert(nord, Tw);
hd = design(d, 'equiripple');       
nfft = 512;
[Hd, w] = freqz(hd.numerator, 1, nfft);

figure(1);    clf;
subplot(121), stem(0:nord, hd.numerator);
title('Einheitspulsantwort des Filters'); xlabel('n');
grid on;   axis tight;
subplot(222), plot((0:nfft-1)/(2*nfft), 20*log10(abs(Hd)));
title('Amplitudengang des Hilbert-Filters');
xlabel('Relative Frequenz f/fs');     grid on;
phase = unwrap(angle(Hd)) + ((0:nfft-1)*nord*pi/(2*nfft))';
subplot(224), plot((0:nfft-1)/(2*nfft), phase);
title('Phasengang ohne den linearen Teil');
xlabel('Relative Frequenz f/fs');     grid on;
La = axis;    axis([La(1:2), -2, 0]);

% ------ Aufruf der Simulation
Tsim = 10;         % Simulationsdauer
%Tsim = 100;        % Simulationsdauer
delta = nord/2;       % Verspätung für die Ausrichtung der Signale
k = 1;           % Zwei Eingangssignale
sim('hilbert_filter2', [0, Tsim]);    % Simulationsaufruf
t = y.time;
y1 = y.Data(:,1);    % Verspätetes Signal
y2 = y.Data(:,2);    % Ausgang des Hilbert-Filters

% -------- Spektrum des analytischen Signals
ys = y.Data(nord*10:end,:);     % Daten ohne Einschwingen
nfft = 256;
Yanal = pwelch((ys(:,1)+j*ys(:,2)), hamming(nfft),50,nfft,fs);
Y1 = pwelch((ys(:,1)+j*eps), hamming(nfft),50,nfft,fs);       
Y2 = pwelch((ys(:,2)+j*eps), hamming(nfft),50,nfft,fs);

figure(2);    clf;
subplot(311), plot((0:nfft-1)/nfft, 10*log10(Y1));
axis tight;  La = axis;    axis([La(1:2), -80, 0.8*La(4)]);
title('Spektrale Leistungsdichte des verspäteten Signals');
xlabel('Relative Frequenz f/fs');    grid on;

subplot(312), plot((0:nfft-1)/nfft, 10*log10(Y2));
axis tight;  La = axis;    axis([La(1:2), -80, 0.8*La(4)]);
title('Spektrale Leistungsdichte am Ausgang des Hilbert-Filters');
xlabel('Relative Frequenz f/fs');    grid on;

subplot(313), plot((0:nfft-1)/nfft, 10*log10(Yanal));
title('Spektrale Leistungsdichte des analytischen Signals');
xlabel('Relative Frequenz f/fs');    grid on;
axis tight;  La = axis;    axis([La(1:2), -100, 0.8*La(4)]);

figure(3);    clf;
subplot(311), plot((-nfft/2:nfft/2-1)/nfft, 10*log10(fftshift(Y1)));
axis tight;  La = axis;    axis([La(1:2), -80, 0.8*La(4)]);
title('Spektrale Leistungsdichte des verspäteten Signals');
xlabel('Relative Frequenz f/fs');    grid on;

subplot(312), plot((-nfft/2:nfft/2-1)/nfft, 10*log10(fftshift(Y2)));
axis tight;  La = axis;    axis([La(1:2), -80, 0.8*La(4)]);
title('Spektrale Leistungsdichte am Ausgang des Hilbert-Filters');
xlabel('Relative Frequenz f/fs');    grid on;

subplot(313), plot((-nfft/2:nfft/2-1)/nfft,...
                     10*log10(fftshift(abs(Yanal))));
title('Spektrale Leistungsdichte des analytischen Signals');
xlabel('Relative Frequenz f/fs');    grid on;
axis tight;  La = axis;    axis([La(1:2), -100, 0.8*La(4)]);

figure(4);    clf;
subplot(311), plot((-nfft/2:nfft/2-1)/nfft, fftshift(Y1));
title('Spektrale Leistungsdichte des verspäteten Signals');
xlabel('Relative Frequenz f/fs');    grid on;

subplot(312), plot((-nfft/2:nfft/2-1)/nfft, fftshift(Y2));
title('Spektrale Leistungsdichte am Ausgang des Hilbert-Filters');
xlabel('Relative Frequenz f/fs');    grid on;

subplot(313), plot((-nfft/2:nfft/2-1)/nfft,...
                     fftshift(abs(Yanal)));
title('Spektrale Leistungsdichte des analytischen Signals');
xlabel('Relative Frequenz f/fs');    grid on;

disp('Mittlere Leistung aus dem Zeitverlauf'), std(y1)^2,
disp('Mittlere Leistung aus der spektralen Leistungsdichte Y1'),...
    sum(Y1)*fs/256,
disp('Mittlere Leistung aus der spektralen Leistungsdichte Y2'),...
    sum(Y2)*fs/256,
disp('Mittlere Leistung aus der spektralen Leistungsdichte Yanal'),...
    sum(Yanal)*fs/256,
