function [output] = window_sinc_filter(order, fc1, fc2, fs, filter_type, window_type, analyse_plot)
% Die Funktion berechnet die Filterkoeffizienten eines FIR Window-Sinc Filters. 
% Hierzu werden die sinc Funktion sin(x)/x und ein wählbares Fenster verwendet.
%
% Inputparameter:
% order        = Filterordnung ->  Anzahl Filterkoeff. = order + 1
%                Bedingung: order = gerade ganze Zahl
%
% fc1 [Hz]     = cutoff frequency 1 = Grenzfrequenz der ersten Filterstufe 
%                -> für Hochpass, Tiefpass, Bandpass und Bandsperre
%                Bedingung: 0 < fc1 < fs/2
%
% fc2 [Hz]     = cutoff frequency 2 = Grenzfrequenz der zweiten Filterstufe 
%                -> neben fc1 zusätzlich nur für Bandpass und Bandsperre
%                Bedingung: 0 < fc1 < fc2 < fs/2
%
% fs [Hz]      = Sampling-frequency = Abtastfrequenz
%
% filter_type  = 'low' = Tiefpass, 'high' = Hochpass,, 'bandpass' = Bandpass
%                und 'bandreject' = Bandsperre
%
% window_type  = Fenstertyp: 'Hamming' oder 'Blackman'
%
% analyse_plot = Bodediagramm, Frequenzantwort linear und Sprungfunktion(nur Tiefpass)
%                des Filters darstellen: 'y' oder 'n'
%
% Outputparameter:
% output     = Filterkoeffizienten = Impulsantwort des Filters
%
% Quelle des Filter-Algorithmus:
% The Scientist and Engineer's Guide to Digital Signal Processing
% By Steven W. Smith, Ph.D.
% Chapter 16: Window-Sinc Filter
% http://www.dspguide.com/ch16.htm
% ------------------------------------------------------------------------- 

% Überprüfung der Inputparameter:
% Abtastfrequenz überprüfen
if (fs < 1) || (~isreal(fs)) 
    error(['Wrong sampling frequency fs [Hz]: ', sprintf('%8.1f',fs),' ...must be a real positiv number']);
end  
% Grenzfrequenz 1 überprüfen
Fc1 = fc1 / fs ; % Normierung auf fn (Fc1 = 0...0.5 fs)
if (Fc1 <= 0) || (Fc1 >= 0.5)
    error(['Wrong cutoff frequency fc1 [Hz]: ', sprintf('%8.1f',fc1),' ...choose: 0 < fc1 < ',sprintf('%8.1f',fs/2),' = fs/2']);
end
% Grenzfrequenz 2 überprüfen
% wird nur für Bandpass- und Bandsperrfilter benötigt
if (strcmp(filter_type, 'bandpass')) || (strcmp(filter_type, 'bandreject'))
    Fc2 = fc2 / fs ; % Normierung auf fn (Fc2 = 0...0.5 fs)
    if ((Fc2 <= 0) || (Fc2 >= 0.5) || (Fc2 <= Fc1))
        error(['Wrong cutoff frequency fc2 [Hz]: ', sprintf('%8.1f',fc2),' ...choose: ',sprintf('%8.1f',fc1) ' < fc2 < ',sprintf('%8.1f',fs/2),' = fs/2']);
    end
end
% Filterordnung überprüfen
M  = order;
if (mod(M, 2) ~= 0) % Filterordnung gerade?
    error(['Wrong Filterorder: ', sprintf('%8.1f',order),' ...must be even']);
end

% Window Typ überprüfen
if (~strcmp(window_type, 'hamming')) && (~strcmp(window_type, 'blackman'))
    error(['Unknown window type: ', window_type,' ...choose: hamming or blackman']);
end

% Auswahl Filtertyp
switch lower(filter_type)
    
    case 'high' % Hochpass - Filter
                
        output = zeros(M+1, 1); % Init
        % Fensterfunktion erstellen
        window = Fenster(M+1, window_type);
        % Berechnung der Filterkoeffizienten
        for i = 0:M 
            if 2 * i == M 
                B(i+1) = 2*pi*Fc1; 
            else
                B(i+1) = sin(2*pi*Fc1 * (i-(M/2))) / (i-(M/2)); 
            end
            B(i+1) = B(i+1) * window(i+1); 
        end                

        % Verstärkungsfaktor des Filters auf 1 normieren
        B = B./sum(B);
        % Filterkoeef. übergeben
        % Tiefpass in Hochpass durch Inversion des Spektrums wandeln
        output      = - B;
        output((M/2)+1) = output((M/2)+1) + 1;
        
    case 'low' % Tiefpass - Filter

        output = zeros(M+1, 1); % Init
        % Fensterfunktion erstellen
        window = Fenster(M+1, window_type);

        % Berechnung der Filterkoeffizienten
        for i = 0:M 
            if 2 * i == M   % Multiplikation ist schneller als Division 
                B(i+1) = 2*pi*Fc1; 
            else 
                B(i+1) = sin(2*pi*Fc1 * (i-(M/2))) / (i-(M/2)); 
            end
            B(i+1) = B(i+1) * window(i+1); 
        end

        % Verstärkungsfaktor des Filters auf 1 normieren
        B      = B./sum(B);
        % Filterkoeef. übergeben
        output = B;

    case 'bandpass' % Bandpass - Filter

        output = zeros(M+1, 1); % Init
        % Fensterfunktion erstellen
        window = Fenster(M+1, window_type);
        
        % 1. Filterstufe: Berechnung der Filterkoeffizienten
        for i = 0:M 
            if 2 * i == M   % Multiplikation ist schneller als Division 
                A(i+1) = 2*pi*Fc1; 
            else 
                A(i+1) = sin(2*pi*Fc1 * (i-(M/2))) / (i-(M/2)); 
            end
            A(i+1) = A(i+1) * window(i+1); 
        end

        % 2. Filterstufe: Berechnung der Filterkoeffizienten
        for i = 0:M 
            if 2 * i == M 
                B(i+1) = 2*pi*Fc2; 
            else 
                B(i+1) = sin(2*pi*Fc2 * (i-(M/2))) / (i-(M/2)); 
            end
            B(i+1) = B(i+1) * window(i+1); 
        end

        % Verstärkungsfaktor des Filters auf 1 normieren
        A = A./sum(A);
        B = B./sum(B);
        % Tiefpass in Hochpass durch Inversion des Spektrums wandeln
        B          = - B;
        B((M/2)+1) = B((M/2)+1) + 1;
        output     = A + B;
        % Filterkoeef. übergeben
        % Bandsperre in Bandpass durch Inversion des Spektrums wandeln
        output          = - output;
        output((M/2)+1) = output((M/2)+1) + 1;

    case 'bandreject' % Bandsperr - Filter

        output = zeros(M+1, 1); % Init
        % Fensterfunktion erstellen
        window = Fenster(M+1, window_type);

        % 1. Filterstufe: Berechnung der Filterkoeffizienten
        for i = 0:M 
            if 2 * i == M   % Multiplikation ist schneller als Division 
                A(i+1) = 2*pi*Fc1; 
            else 
                A(i+1) = sin(2*pi*Fc1 * (i-(M/2))) / (i-(M/2)); 
            end
            A(i+1) = A(i+1) * window(i+1); 
        end

        % 2. Filterstufe: Berechnung der Filterkoeffizienten
        for i = 0:M 
            if 2 * i == M   
                B(i+1) = 2*pi*Fc2; 
            else 
                B(i+1) = sin(2*pi*Fc2 * (i-(M/2))) / (i-(M/2)); 
            end
            B(i+1) = B(i+1) * window(i+1); 
        end

        % Verstärkungsfaktor des Filters auf 1 normieren
        A = A./sum(A);
        B = B./sum(B);
        % Tiefpass in Hochpass durch Inversion des Spektrums wandeln
        B          = - B;
        B((M/2)+1) = B((M/2)+1) + 1;
        % Filterkoeef. übergeben
        output = A + B;
        
    otherwise % ungültiger Filtertyp
        error(['Unknown filter type: ', filter_type ' ...choose: high, low, bandpass or bandreject']);
end % switch

% Filter Frequenzantwort darstellen?
% Eingabewert überprüfen
if (~strncmpi(analyse_plot, 'y', 1)) && (~strncmpi(analyse_plot, 'n', 1))
    error(['Unknown analyse_plot command: ', analyse_plot ,' ...choose: y or n']);
elseif (strncmpi(analyse_plot, 'y', 1))
    % Plots erstellen
     filter_analyse(output, fs, filter_type);
end