function filter_analyse(h, fs, filter_type)
%
% Die Funktion erstellt das Bodediagramm der Impulsantwort des Filters h
% Die reale/en Grenzfrequenz/en fc1 und fc2 [Hz] werden grafisch bestimmt
% 
% Input:
% h           = Impulsantwort des Filters
% fs [Hz]     = Sampling frequency = Abtastfrequenz
% filter_type = 'low' = Tiefpass, 'high' = Hochpass,, 'bandpass' = Bandpass
%               und 'bandreject' = Bandsperre
%
% Output:
% Darstellung der realen Grenzfrequenz 1 und 2 (2 - wenn benutzt) im command window
% Genauigkeit 0.1 Hz
% -------------------------------------------------------------------------

% Anzahl Messwerte
N = 2^15; % 32768 Messwerte
% Frequenzspektrum mit FFT berechnen
H = fft(h, N);
% Berechnung des Amplitudengangs aus dem komplexen Frequenzvektor H:
mag = abs(H(2:(N/2)+1)); % Darstellung von df...fn [Hz]
mag = 20*log10(mag); % Umrechnung in dB
% Auflösung des Spektrums
df = fs/N;
% Nyquistfrequenz
fn = fs/2; 
% Frequenzvektor df...fn [Hz]
f = df : df : fn;
% Berechnung der Phase in Grad
phase = angledim(unwrap(angle(H(2:(N/2)+1))),'radians','degrees');

% Ausgabe des Bodediagramms
figure('Name','Bodediagram - Filter Frequency response','NumberTitle','off')
% Darstellung Betrag 
subplot(211); 
semilogx(f,mag,'b','LineWidth',2);

% 3dB Grenzfrequenz/-en grafisch bestimmen
% Grenzfrequenz = Abfall der Verstärkung um 1/sqrt(2) = - 3.0103 dB
dbwert = 20*log10(1/sqrt(2));

switch lower(filter_type)
    case 'high' % Hochpass
        % Grenzfrequenz suchen
        index1 = find(mag>=dbwert,1,'first');
        disp(sprintf('\n\t %s %8.2f \n',['real filter cutoff frequency fc1: '],f(index1-1)));
        hold on;
        plot(f(index1-1),mag(index1-1),'ro','LineWidth',2); % Grenzfrequenz fc1 markieren
        hold off;
        text_fc1 = sprintf(['real cutoff frequency fc1 = ' num2str(f(index1-1),'%8.2f'),' Hz']);
        if f(index1-1) > 0.5*fn
            legend('filter frequency response',text_fc1,'Location','NorthWest');
        else
            legend('filter frequency response',text_fc1,'Location','SouthEast');
        end
        
    case 'low' % Tiefpass
        % Grenzfrequenz suchen
        index1 = find(mag<=dbwert,1,'first');
        disp(sprintf('\n\t %s %8.2f \n',['real filter cutoff frequency fc1: '],f(index1-1)));
        hold on;
        plot(f(index1-1),mag(index1-1),'ro','LineWidth',2); % Grenzfrequenz fc1 markieren
        hold off;
        text_fc1 = sprintf(['real cutoff frequency fc1 = ' num2str(f(index1-1),'%8.2f'),' Hz']);
        if f(index1-1) > 0.5*fn
            legend('filter frequency response',text_fc1,'Location','SouthWest');
        else
            legend('filter frequency response',text_fc1,'Location','NorthEast');
        end
    
    case 'bandpass'
        % 1. Grenzfrequenz suchen
        index1 = find(mag>=dbwert,1,'first');
        disp(sprintf('\n\t %s %8.2f \n',['real filter cutoff frequency fc1: '],f(index1-1)));
        hold on;
        plot(f(index1-1),mag(index1-1),'ro','LineWidth',2); % Grenzfrequenz fc1 markieren
        if (strcmp(filter_type, 'bandpass')) || (strcmp(filter_type, 'bandreject'))
            % 2. Grenzfrequenz suchen       
            index2 = find(mag(index1:end)<=dbwert,1,'first');
            disp(sprintf('\n\t %s %8.2f \n',['real filter cutoff frequency fc2: '],f(index1+index2-2)));
            plot(f(index1+index2-2),mag(index1+index2-2),'go','LineWidth',2); % Grenzfrequenz fc1 markieren
            hold off;
            text_fc1 = sprintf(['real cutoff frequency fc1 = ' num2str(f(index1-1),'%8.2f'),' Hz']);
            text_fc2 = sprintf(['real cutoff frequency fc2 = ' num2str(f(index1+index2-2),'%8.2f'),' Hz']);        
            legend('filter frequency response',text_fc1, text_fc2,'Location','NorthWest');
        end
    
    case 'bandreject'
        % 1. Grenzfrequenz suchen
        index1 = find(mag<=dbwert,1,'first');
        disp(sprintf('\n\t %s %8.2f \n',['real filter cutoff frequency fc1: '],f(index1-1)));
        hold on;
        plot(f(index1-1),mag(index1-1),'ro','LineWidth',2); % Grenzfrequenz fc1 markieren
        % 2. Grenzfrequenz suchen
        index2 = find(mag(index1:end)>=dbwert,1,'first');
        disp(sprintf('\n\t %s %8.2f \n',['real filter cutoff frequency fc2: '],f(index1+index2-2)));
        plot(f(index1+index2-2),mag(index1+index2-2),'go','LineWidth',2); % Grenzfrequenz fc1 markieren
        hold off;
        text_fc1 = sprintf(['real cutoff frequency fc1 = ' num2str(f(index1-1),'%8.2f'),' Hz']);
        text_fc2 = sprintf(['real cutoff frequency fc2 = ' num2str(f(index1+index2-2),'%8.2f'),' Hz']);        
        legend('filter frequency response',text_fc1, text_fc2,'Location','SouthWest');
    otherwise 
        error(['Unknown filter type in bode diagram: ', filter_type]);    
end % end switch

axis([10^-1 fn -200 20]);
title('Bode Diagram - Filter frequency response','FontSize',18); 
xlabel('Frequency [Hz]','FontSize',14); 
ylabel('Magnitude [dB]','FontSize',14); 
grid on; 
% Darstellung Phase 
subplot(212);
semilogx(f,phase,'LineWidth',2);
hold on;
plot(f(index1-1),phase(index1-1),'ro','LineWidth',2); % Grenzfrequenz fc1 markieren
if (strcmp(filter_type, 'bandpass')) || (strcmp(filter_type, 'bandreject'))
    plot(f(index1+index2-2),phase(index1+index2-2),'go','LineWidth',2); % Grenzfrequenz fc2 markieren
end
hold off;
axis([10^-1 fn 1.01*min(phase) 1000]);
xlabel('Frequency [Hz]','FontSize',14); 
ylabel('Phase [deg]','FontSize',14); 
grid on;

% Lineare Darstellung
figure('Name','Filter Frequency response','NumberTitle','off')
f = 0 : df : fn;
plot(f,abs(H(1:(N/2)+1)),'LineWidth',2);
if (strcmp(filter_type, 'low')) || (strcmp(filter_type, 'high'))
    %nur eine Grenzfrequenz
    if 2*f(index1-1) > fn
        f_end = fn;
    else
        f_end = 2*f(index1-1);
    end
    axis([0 f_end 0 1.2])    
else
    if 2*f(index1+index2-2) > fn
        f_end = fn;
    else
        f_end = 2*f(index1+index2-2);
    end
    axis([0.5*f(index1-1) f_end 0 1.2])
end
title('Filter frequency response','FontSize',18); 
xlabel('Frequency [Hz]','FontSize',14); 
ylabel('Amplitude','FontSize',14); 
grid on;

% Sprungantwort für Tiefpass - Filter
if (strcmp(filter_type, 'low'))
    figure('Name','Filter step response','NumberTitle','off')
    % Länge der Sprungfunktion abhängig von der Länge des Filters 
    if (length(h) < 1024)
        N = 1024;
    else
        N = 2^nextpow2(length(h));
    end
    % Sprungfunktion
    stepfunction = [ones(N,1)]; 
    % Sprungantwort erstellen
    step_response = FFT_Faltung(stepfunction,h);     
    plot(step_response(1:N),'LineWidth',2);
    title('Filter step response','FontSize',18); 
    xlabel('Sample number','FontSize',14); 
    ylabel('Amplitude','FontSize',14);
    grid on;
end