function [x] = ir2112(ir, fs, delta)
%IR2112
% ir: impulzusválasz
% fs: mintavételi frekvencia
% delta: Ha delta>0, akkor addíg fut az illesztés, amíg az optimumtól 
% távolabb vagyunk mint delta, különben az intervallum b?vít?s módszert 
% használjuk 

    % Illesztési intervallumok, dB-ben kifejezve
    thld_low = -5;
    thld_mid = -10;
    thld_hig = -35;

    % Az IR hossza mintákban
    K = length(ir);
    
    % Id?vektor
    t = linspace(0, K-1, K)/fs;
    L = t(end);
    
    % EDC görbe
    ir_sqr = ir.^2;
    prev = 0;
    tmp = zeros(1, K);
    for i = 1:K,
        tmp(i) = ir_sqr(i)+prev;
        prev = tmp(i);
    end;
    ir_edc = 1-tmp/sum(ir_sqr);
    
    % Logarimikus EDC
    log_edc = 10*log10(ir_edc);
    
    % Megbecsüljük t_60-at thld_low és thld_mid intervallumból
    k_low = find(log_edc <= thld_low, 1);
    k_mid = find(log_edc <= thld_mid, 1);
    k_hig = find(log_edc <= thld_hig, 1);
    k_60 = (-60-thld_low)*(k_mid-k_low)/(thld_mid-thld_low);
    t_60 = k_60/fs;
    
    % Két m?ködési mód van, ha delta <= 0, akkor az intervallum b?vítget?s 
    % módszert használjuk, ha delta > 0, akkor addig keresünk, amíg az 
    % eltérés az optimumtól nagyobb mint delta 
    if delta > 0,
        
        % Az illesztést a teljes edc-re végezzük
        l = 1;
        m = K;
        
        % Kiindulási paraméterek (x3-ra kell még egy becslés)
        x = [1, 3*log(10)/t_60, 0];
        
        % Delta x
        dx = [Inf Inf Inf];

        % Iteratív algorimtus
        ok = 1;
        while ok && (norm(dx) > delta),

            % Modell
            Fx = x(1)*exp(-x(2)*t(l:m))+x(3)*(L-t(l:m));

            % Hiba
            E = Fx-ir_edc(l:m);

            % A modell deriváltjai az x1, x2, x3 paraméterek szerint
            Fdx1 = exp(-x(2)*t(l:m));
            Fdx2 = -x(1)*t(l:m).*Fdx1;
            Fdx3 = L - t(l:m);

            % Kiszámoljuk az A mátrixot
            a11 = sum(Fdx1.*Fdx1);
            a12 = sum(Fdx1.*Fdx2);
            a13 = sum(Fdx1.*Fdx3);
            a22 = sum(Fdx2.*Fdx2);
            a23 = sum(Fdx2.*Fdx3);
            a33 = sum(Fdx3.*Fdx3);
            A = [a11, a12, a13; a12, a22, a23; a13, a23, a33]; 

            % Kiszámoljuk a b vektort
            b = [sum(E.*Fdx1); sum(E.*Fdx2); sum(E.*Fdx3)];

            % Delta x
            dx_new = pinv(A)*b;

            % Ha delta x csökken, akkor konvergálunk
            if norm(dx_new) < norm(dx),

                % Elmozdulunk a paramétertérben
                x = x - dx_new';
                dx = dx_new;
                
            else

                % Nem konvergál az algoritmus, kilépünk
                ok = 0;

            end;

        end;        
        
    else
       
        % Az illesztést azon az intervallumon kezdjük, amin becsültünk, 
        % majd k_hig-ig b?vítjük
        l = k_low;
        m = k_mid;
        
        % A b?vítés mértéke
        dm = max([1 round((k_hig-k_mid)/10)]);
        
        % Kiindulási paraméterek (x3-ra kell még egy becslés)
        x = [1, 3*log(10)/t_60, 0];
        
        % Delta x
        dx = [Inf Inf Inf];

        % Iteratív algorimtus
        ok = 1;
        while ok && m < k_hig,

            % Modell
            Fx = x(1)*exp(-x(2)*t(l:m))+x(3)*(L-t(l:m));

            % Hiba
            E = Fx-ir_edc(l:m);

            % A modell deriváltjai az x1, x2, x3 paraméterek szerint
            Fdx1 = exp(-x(2)*t(l:m));
            Fdx2 = -x(1)*t(l:m).*Fdx1;
            Fdx3 = L - t(l:m);

            % Kiszámoljuk az A mátrixot
            a11 = sum(Fdx1.*Fdx1);
            a12 = sum(Fdx1.*Fdx2);
            a13 = sum(Fdx1.*Fdx3);
            a22 = sum(Fdx2.*Fdx2);
            a23 = sum(Fdx2.*Fdx3);
            a33 = sum(Fdx3.*Fdx3);
            A = [a11, a12, a13; a12, a22, a23; a13, a23, a33]; 

            % Kiszámoljuk a b vektort
            b = [sum(E.*Fdx1); sum(E.*Fdx2); sum(E.*Fdx3)];

            % Delta x
            dx_new = pinv(A)*b;

            % Ha delta x csökken, akkor konvergálunk
            if norm(dx_new) < norm(dx),

                % Elmozdulunk a paramétertérben
                x = x - dx_new';
                dx = dx_new;
                
                % B?vítjük az intervallumot
                m = m + dm;
                
            else

                % Nem konvergál az algoritmus, kilépünk
                ok = 0;

            end;

        end;
        
    end;
    
end

