% Skript inv_sinc_1.m, in dem ein inverses sinc-Filter
% untersucht wird

clear;
% ------ Frequenzgang des Halteglied-Nullter-Ordnung
fs = 1000;       Ts = 1/fs;
nfft = 1024;     df = fs/nfft;
f = -fs/2:df:fs/2-df;
H = exp(-j*pi*f*Ts).*sinc(f*Ts);    % Frequenzgang des Halteglieds nullter
            % Ordnung geteil durch Ts
% ------ Einfaches Kompensationsfilter der Ordnung 2
h = [-1, 18 -1]/16;
z = exp(j*2*pi*f*Ts);
Hk = polyval(h,z)./(z.^2);
figure(1);     clf;
plot(f/fs, abs(H));     hold on;
plot(f/fs, abs(Hk),'k');
plot(f/fs, (abs(H).*abs(Hk)),'k');
title(['Amplitudengang des Halteglieds nullter Ordnung,'...
    ' des Kompensationsfilters und das Ergebnis']);
xlabel('Relative Frequenz f/fs'); grid on;

% ------ Aufruf der Simulation
Tsim = 1;        % Simulationsdauer
delta = 1;       % Verspätung für die Ausrichtung der Signale
sim('inv_sin1', [0, Tsim]);    % Simulationsaufruf
t = y.time;
y1 = y.Data(:,1);    % Analoger Eingang
y2 = y.Data(:,2);    % Ausgang des Filters
y3 = y.Data(:,3);    % Abgetastetes Signal

figure(2);    clf;
plot(t, [y1,y3,y2]);
La = axis;    axis([0.5, 0.54, La(3:4)]);
title('Analoger Eingang, Abtastwerte und Ausgang des Filters');
xlabel('Zeit in s');     grid on;
legend('Analoger Eingang','Abgetastetes Signal','Ausgang des Filters');
% ------- Amplitudenspektrum FFT/N
nord = 2;
ys = y.Data(nord*100:end,:);   % Ohne Einschwingen
nY = length(ys(:,1));
Y = fft(ys,nY)/nY;       

figure(3);    clf;
subplot(311), plot((-(nY-1)/2:(nY-1)/2)*10*fs/nY,fftshift(abs(Y(:,1)))); 
La = axis;    axis([-2000, 2000, 0, 0.55])
title(['Amplitudenspektrum des analogen Eingangssignals']);
xlabel(['Hz fs = ',num2str(fs)]);     grid on;

subplot(312), plot((-(nY-1)/2:(nY-1)/2)*10*fs/nY,fftshift(abs(Y(:,3)))); 
Hn = sinc((-(nY-1)/2:(nY-1)/2)/(0.1*nY));     hold on
plot((-(nY-1)/2:(nY-1)/2)*fs/(0.1*nY), abs(Hn)*0.5,'r');
La = axis;    axis([-2000, 2000, La(3:4)]);
title(['Amplitudenspektrum des D/A gewandelten Signals ohne Kompensationsfilter']);
xlabel('Hz');     grid on;

subplot(313), plot((-(nY-1)/2:(nY-1)/2)*10*fs/nY,fftshift(abs(Y(:,2)))); 
La = axis;    axis([-2000, 2000, La(3:4)])
title(['Amplitudenspektrum des D/A gewandelten Signals mit einfachen Kompensationsfilter']);
xlabel('Hz');     grid on;


% ------- Inverses sinc-Filter der Ordnung 50
fp = 0.4*2;    fst = 0.45*2;
dp = 0.001;    dst = 40;
nord = 50;
d = fdesign.isinclp('N,Fp,Fst',nord,fp,fst);
hd = design(d,'SincFrequencyFactor',0.5,'SincPower',1);
hd = hd.numerator;
%hd = firceqrip(50,0.8,[0.01 0.01],'invsinc',[0.42 1.5]);

figure(4);   clf;
[Hd, w] = freqz(hd,1,'whole');
f = (-256:255)/512;
H = exp(-j*pi*f).*sinc(f);
plot((-256:255)/512, fftshift(abs(Hd)));    hold on;
plot((-256:255)/512, abs(H),'k');
plot((-256:255)/512, abs(H)'.*fftshift(abs(Hd)),'k');
title(['Amplitudengang des Halteglieds nullter Ordnung,'...
    ' des Kompensationsfilters und das Ergebnis']);
xlabel('Relative Frequenz f/fs');   grid on;

% ------ Aufruf der Simulation
Tsim = 1;       delta = nord/2;      h = hd;
sim('inv_sin1', [0, Tsim]);
t = y.time;
y1 = y.Data(:,1);    % Analoger Eingang 
y2 = y.Data(:,2);    % Ausgang des Filters 
y3 = y.Data(:,3);    % Abgetastetes Signal 

figure(5);    clf;
plot(t, y1, t,y3,'k', t, y2, 'k');
%plot(t, [y1, y3, y2]);
La = axis;    axis([0.5, 0.54, La(3:4)]);
title('Analoger Eingang, Abtastwerte und Ausgang des Kompensationsfilters');
xlabel('Zeit in s');     grid on;
legend('Analoger Eingang','Abgetastetes Signal','Ausgang des Filters');

% ------- Amplitudenspektrum FFT/N
ys = y.Data(nord*100:end,:);   % Ohne Einschwingen 
nY = length(ys(:,1));      Y = fft(ys,nY)/nY;       

figure(6);    clf;
subplot(311), plot((-(nY-1)/2:(nY-1)/2)*10*fs/nY,fftshift(abs(Y(:,1)))); 
La = axis;    axis([-2000,2000, 0, 0.55])
title(['Amplitudenspektrum des analogen Eingangssignals']);
xlabel(['Hz fs = ',num2str(fs)]);     grid on;

subplot(312), plot((-(nY-1)/2:(nY-1)/2)*10*fs/nY,fftshift(abs(Y(:,3)))); 
Hn = sinc((-(nY-1)/2:(nY-1)/2)/(0.1*nY));    hold on
plot((-(nY-1)/2:(nY-1)/2)*fs/(0.1*nY), abs(Hn)*0.5,'r');
La = axis;    axis([-2000,2000, La(3:4)])
title('Amplitudenspektrum des D/A gewandelten Signals ohne Kompensationsfilter');
xlabel('Hz');     grid on;

subplot(313), plot((-(nY-1)/2:(nY-1)/2)*10*fs/nY,fftshift(abs(Y(:,2)))); 
La = axis;    axis([-2000,2000, 0, 0.55])
title(['Amplitudenspektrum des D/A gewandelten Signals mit Kompensationsfilters']);
xlabel('Hz');     grid on;
