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

% M = melbank(fs, n, p, eq)
% 
% Funkcia ráti banku filtrov rovnomerne rozložených v Mel 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:0) 
% M     - výstupná matica, každý riadok sa viaže na filter

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

%% Výpočet centra frekvencií
c_m = zeros(3, p);
c_m(2, :) = (1:p).*(hz2mel(fix(fs/2))/(p + 1));
c_m(1, :) = c_m(2, :) - c_m(2, 1);
c_m(3, :) = c_m(2, :) + c_m(2, 1);

c_f = mel2hz(c_m);

%% Filter banky
M = zeros(p, n);
f = linspace(0, fix(fs/2), n);

for i = 1:p
    ind1 = sum(f < c_f(1, i)) + 1;
    ind2 = sum(f <= c_f(2, i));
    
    M(i, ind1:ind2) = (f(ind1:ind2) - c_f(1, i)) ./ ...
        (c_f(2, i) - c_f(1, i));
    
    ind1 = sum(f <= c_f(2, i));
    ind2 = sum(f <= c_f(3, i));
    
    M(i, ind1:ind2) = (f(ind1:ind2) - c_f(3, i)) ./ ...
        (c_f(2, i) - c_f(3, i));
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