Função para Ajuste Sinusoidal para Gretl

function [K,A,f,phi,rmq]=num_ajuste_sinusoidal(x,y)

    // [K,A,f,phi,rmq]=num_ajuste_sinusoidal(x,y)

    // Calcula os coeficientes de ajuste sinusoidal g(x)=A*sin(f*(x-phi))+K

    // para os pontos (x_i,y_i) e o resíduo do ajuste de mínimos quadrados

    // (C) Rudnei Dias da Cunha 2014

    m = max(size(x))

    K = sum(y)/m

    Y = y-K

    A = max(abs(Y))

    Y = asin(Y/A)

    H = [ m sum(Y) ; sum(Y) sum(Y.^2) ]

    b = [ sum(x) ; sum(x.*Y) ]

    a = H\b

    a_0 = a(1)

    a_1 = a(2)

    phi = a_0

    if abs(a_1)>0.0 then

        f = 1.0/a_1

        rmq = sqrt(sum((num_avalia_sinusoidal(x,K,A,f,phi)-y).^2))

    else

        f = []

        rmq = []       

    end

endfunction