getPeakCoeff(fs,fc,Q,Gain_dB,mode)

function [b,a] = getPeakCoeff(fs,fc,Q,Gain_dB,mode)
% ******************************************************
% * Generates the filter coefficients for a Peak-Filter 
% * ----------------------------------------------------
% * Usage : getPeakCoeff(fs,fc,Q,GAin_dB,mode)
% * 
% *         fs      -   sampling frequency
% *         fc      -   center frequency
% *         Q       -   Q-factor
% *         gain_dB -   amplification
% *         mode    -   Udo Zoelzer ['Zoel'] or 
% *                     Robert Bristow Johnson ['RBJ']
% * ----------------------------------------------------
% * Copyright (c) 2004
% * Author  :   Tobias May
% * Date    :   11 Nov 2004
% * Version 1.0
% *
% * Fachhochschule OOW Standort Oldenburg
% * Studiengang Hrtechnik und Audiologie
% ******************************************************

if(fc>(fs/2))
    error('Cutoff frequency greater than Nyquistfrequency');
end

switch mode
    case 2
        
        K = tan(pi*fc/fs); 
        T = 1/fs;
        K = tan(2 * pi * fc* (T/2));                                     
        V0 = 10^(Gain_dB/20);     
        
        b0 = (1 + V0/Q * K + K^2)/(1 + 1/Q * K + K^2);                                            
        b1 = (2 * (K^2 - 1))/(1 + 1/Q * K + K^2); 
        b2 = (1 - V0/Q * K + K^2)/(1 + 1/Q * K + K^2);  
        a0 = 1;
        a1 = (2 * (K^2 - 1))/(1 + 1/Q * K + K^2); 
        a2 = (1 - 1/Q * K + K^2)/(1 + 1/Q * K + K^2);                                              
    
        b = [b0, b1, b2];
        a = [a0, a1, a2];
   
        
    case 1
        
        Ommega_rad = 2 * pi * fc / fs ;
        BandWidth_Oct = log((2*Q^2 +1)/(2*Q^2) + sqrt(((2*Q^2 + 1)/Q^2)^2 / 4 -1))/log(2);
       
        V0 = 10^(Gain_dB/20);
        
        Gamma = sinh((log(2)/2) * BandWidth_Oct * (Ommega_rad/sin(Ommega_rad))) * sin(Ommega_rad);
        
        b0 = (1 + Gamma * sqrt(V0)) / (1 + (Gamma/sqrt(V0))); 
        b1 = ( - 2 * cos(Ommega_rad))    / (1 + (Gamma/sqrt(V0)));
        b2 = ( 1 - Gamma * sqrt(V0)) / (1 + (Gamma/sqrt(V0)));
        a0 = 1;
        a1 = ( -2 * cos(Ommega_rad)) / (1 + (Gamma/sqrt(V0)));
        a2 = ( 1 - (Gamma/sqrt(V0)))/ (1 + (Gamma/sqrt(V0)));
        
        b = [b0, b1, b2];
        a = [a0, a1, a2];
        
         
    otherwise
        error('Wrong filter type');
end

%FilterAnalysis(b,a,fs,sprintf('%s',name));