% Skript ideale_filter_2.m, in dem die Annäherung
% des idealen BP-Filters als FIR-Filter 
% untersucht wird

clear;
% ------- Parameter des Frequenzgang des Filters
fs = 1000;      % Abtastfrequenz
Ts = 1/fs;      % Abtastperiode
fp1 = 100;      % Durchlassfrequenz
fst1 = 300;     % Stopfrequenz 
N = 127;        % Ungerade Anzahl der Koeffizienten
                % das gerade Ordnung bedeutet      
n = 0:N-1;
n1 = round(fp1/(fs/N));    n2 = round(fst1/(fs/N)); 
                % Indizes im ersten Nyquist-Intervall 
n3 = N-n2;                 n4 = N-n1;
                % Indizes im gespiegelten Bereich
Hid_BP = (n1<=n)&(n<=n2) + (n3<=n)&(n<=n4);
                % Idealer Frequenzgang (abgetastet)
% ------- Ideale begrenzte Einheitspulsantwort aus der inverse DFT
hid_BP = ifft(Hid_BP);    % Einheitspulsantwort
figure(1);    clf;
subplot(211), stem(n/N, Hid_BP);
title('Idealer Frequenzgang');
xlabel('f/fs');    grid on;
subplot(212), stem(-(N-1)/2:(N-1)/2, fftshift(real(hid_BP)));
title(['Einheitspulsantwort des BP-Filters über die ifft ']);
grid on; axis tight;

% ------- DTFT aus der Einheitspulsantwort mit dichterem Gitter 
N1 = 512;
n1 = (-N1/2:N1/2-1);
f = n1*fs/N1;    nf = length(f);
hid = fftshift(real(hid_BP));    % Erzeugung der nicht kausalen
                                 % Einheitspulsantwort um Index 0
Hid = zeros(1,nf);   % DTFT mit dichterem Frequenzraster
for k = 1:nf
    Hid(k) = sum(real(hid).*exp(-j*2*pi*f(k)*(-N/2:N/2-1)*Ts));
end;
figure(2);    clf;
subplot(211), plot(n1/N1, abs(Hid));
hold on; stem((-(N-1)/2:(N-1)/2)/N, abs(fftshift(fft(hid))));
hold off;    xlabel('f/fs');    grid on;
title('DTFT der Einheitspulsantwort und die iDFT');
subplot(212), plot(n1/N1, abs(Hid));
hold on; stem((-(N-1)/2:(N-1)/2)/N, abs(fftshift(fft(hid))));
hold off;    xlabel('f/fs');    grid on;
title('DTFT der Einheitspulsantwort und die iDFT (Ausschnitt)');
La = axis;   axis([0.05, 0.35, La(3:4)]);
% ------- Komplexer frequenzgang über die Funktion freqz und fft
            % mit Nullerweiterung
[Hf, w] = freqz(hid,1,N1,'whole');
figure(3);    clf;
subplot(211), plot((w-3.14)/(2*pi), abs(fftshift(Hf)));
title('Die DTFT über die Funktion freqz mit Nullerweiterung');
xlabel('f/fs');    grid on;
subplot(212), plot(n1/N1, abs(fftshift(fft(hid, N1))));
title('Die DTFT über die fft mit Nullerweiterung');
xlabel('f/fs');    grid on;

