%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Lundeby.m
%
%Esta funcao realiza as iteracoes do metodo de Lundeby.
%Aplicar metodo de Lundeby para determinar ponto de truncamento.
%
%[ponto,C,rms_dB]=lundeby(IR,fs,flag)
%
%Tem como entrada a resposta impulsiva e a taxa de amostragem
%Tem como saida o ponto de cruzamento da assintota de decaimento o nivel de
%ruido e se desejado, a constante C usada para compensar o truncamento no
%calculo da curva de Schroeder.
%Se nao for fornecida variavel de saida, a funcao imprime um grafico dos
%niveis encontrados.

function [ponto,C,rms_dB]=lundeby2(ir,fs,flag)

energia_ir = ir.^2;

rms_dB = 10*log10(mean(energia_ir(round(.9*length(energia_ir)):end))/max(energia_ir));

t = floor(length(energia_ir)/fs/0.01);
v = floor(length(energia_ir)/t);

for n=1:t, % ablakozott simitas
    media(n) = mean(energia_ir((((n-1)*v)+1):(n*v)));
    eixo_tempo(n) = ceil(v/2)+((n-1)*v);
end

plot(eixo_tempo,media);
mediadB = 10*log10(media/max(energia_ir));

r = find(mediadB > rms_dB+10, 1, 'last');
if any (mediadB(1:r) < rms_dB+10), r = find(mediadB(1:r) < rms_dB+10, 1); end

if isempty(r), r=10;
elseif r<10, r=10; end;

[A,B] = intlinear(eixo_tempo(1:r),mediadB(1:r));
cruzamento = (rms_dB-A)/B;
%display(cruzamento);

if rms_dB > -20
    %Relacao sinal ruido insuficiente
    ponto=length(energia_ir);
    if nargout==2
        C=0;
    end
else
    
    %%%%%%%%%%%%%%%%%%%%%%%%INICIA A PARTE ITERATIVA DO PROCESSO%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    erro=1;
    INTMAX=50;
    vezes=1;
    while (erro > 0.0001 && vezes <= INTMAX)

        %Calcula novos intervalos de tempo para media, com aproximadamente p passos por 10dB
        clear r t v n media eixo_tempo;

        p = 5;                          %numero de passos por decada
        delta = abs(10/B);              %numero de amostras para o a linha de tendencia decair 10dB
        v = floor(delta/p);             %intervalo para obtencao de media
        t = floor(length(energia_ir(1:round(cruzamento-delta)))/v);

%display(t);
%display(v);
%display(delta);

        if t < 2                        %numero de intervalos para obtencao da nova media no intervalo
            t=2;                        %que vai do inicio ate 10dB antes do ponto de cruzamento.
        elseif isempty(t)
            t=2;
        end
    
        for n=1:t
            media(n) = mean(energia_ir((((n-1)*v)+1):(n*v)));
            eixo_tempo(n) = ceil(v/2)+((n-1)*v);
        end

%display(media);
%display(eixo_tempo);

        mediadB = 10*log10(media/max(energia_ir));
    
        clear A B noise energia_ruido rms_dB;
        [A,B] = intlinear(eixo_tempo,mediadB);

%display(cruzamento);
        %nova media da energia do ruido, iniciando no ponto da linha de tendencia 10dB abaixo do cruzamento.
        noise = energia_ir(round(cruzamento+delta):end);
        if (length(noise) < round(.1*length(energia_ir)))
            noise = energia_ir(round(.9*length(energia_ir)):end); 
        end       

%display(length(noise));
        rms_dB = 10*log10(mean(noise)/max(energia_ir));
%display(rms_dB);
%display([A B]);

        %novo ponto de cruzamento.
        erro = abs(cruzamento - (rms_dB-A)/B)/cruzamento;
        cruzamento = round((rms_dB-A)/B);
        vezes = vezes + 1;

        %display(erro);
        %display(cruzamento);

    end
end

if nargout == 1
    if cruzamento > length(energia_ir)     %caso o sinal nao atinja o patamar de ruido
        ponto = length(energia_ir);        %nas amostras fornecidas, considera-se o ponto
    else                                        %de cruzamento a ultima amosta, o que equivale
        ponto = cruzamento;                     %a nao truncar o sinal.
    end
else
    if cruzamento > length(energia_ir)
        ponto = length(energia_ir);
    else
        ponto = cruzamento;
    end
    C=max(energia_ir)*10^(A/10)*exp(B/10/log10(exp(1))*cruzamento)/(-B/10/log10(exp(1)));
end

if (nargout == 0 || flag == 1)
    figure
    plot((1:length(energia_ir))/fs,10*log10(energia_ir/max(energia_ir)));
    hold
%    stairs(eixo_tempo/fs,mediadB,'y');
    plot((1:cruzamento+1000)/fs,A+(1:cruzamento+1000)*B,'k');
    line([cruzamento-1000,length(energia_ir)]/fs,[rms_dB,rms_dB],'Color','k');
    plot(cruzamento/fs,rms_dB,'o','MarkerFaceColor','w','MarkerSize',10);
    hold
end
