function M = barkbank(fs, n, p, eq)

% M = barkbank(fs, n, p, eq)
% 
% Funkcia, ktorá vráti banku filtrov rovnomerne rozloženú v Barkovej
% stupnici
% fs    - vzorkovacia frekvencia
% n     - počet vzorkov na filter
% p     - počet filtrov v banke(default: depending on fs)
% eq    - ak je nastavený na 1 budú použité krivky rovnakej hlasitosti (default:1)
% M     - výstupná matica, každý riadok sa viaže na filter

%% Cesty a premenné
if ((nargin < 3) || (isempty(p)))
    if (fs <= 8e3)
        p = 17;
    elseif (fs <= 11e3)
        p = 19;
    elseif (fs <= 16e3)
        p = 21;
    elseif (fs <= 22e3)
        p = 23;
    else
        p = 27;
    end
end
if ((nargin < 4) || (isempty(eq)))
    eq = 1;
end

%% Výpočet centra frekvencií
c_b = zeros(5, p);
c_b(3, :) = (0:(p - 1)).*(hz2bark(ceil(fs/2))/(p - 1));
c_b(1, :) = c_b(3, :) - 2.5;
c_b(2, :) = c_b(3, :) - 0.5;
c_b(4, :) = c_b(3, :) + 0.5;
c_b(5, :) = c_b(3, :) + 1.3;

c_f = bark2hz(c_b);

%% Príprava frekvenčného merítka
f = linspace(0, ceil(fs/2), n);

%% Príprava banky
M = zeros(p, n);

for i = 1:p
    ind1 = sum(f <= c_f(1, i));
    ind2 = sum(f <= c_f(2, i));

    if (ind2 > 0)
        if (ind1 == 0)
            ind1 = 1;
        end
        
        M(i, ind1:ind2) = ...
            10.^(linspace(-2.5, -0.5, (ind2 - ind1 + 1)) + 0.5);
    end
    
    ind1 = sum(f <= c_f(2, i));
    ind2 = sum(f <= c_f(4, i));
        
    if (c_f(4, i) > 0)
        if (ind1 == 0)
            ind1 = 1;
        end
        
        M(i, ind1:ind2) = 1;
    end
    
    ind1 = sum(f <= c_f(4, i));
    ind2 = sum(f <= c_f(5, i));

    M(i, ind1:ind2) = ...
        10.^(-2.5*(linspace(0.5, 1.3, (ind2 - ind1 + 1)) - 0.5));
end

%% Ak je potreba, aplikovať krivky hlasitosti
if (eq)
    E = eql_curves(fs, n);
    
    tmp = c_f(3, :)';
    tmp = f(ones(p, 1), :) <= tmp(:, ones(1, n));
    
    ind = sum(tmp.');    
    m   = E(ind).';
    m   = m(:, ones(1, n));
    
    M = M.*m;
end