function [f,eta,fmap]=rec_fiaa_oct_c1(x,K,q_i,vt,f_in,eta_in)
%
% Recursive fast IAA
%
% x => Input data vector
% K => Number of frequency points
% q_i => Number of IAA iterations
% vt => noise weighting factor (1 is default)
% f_in => IAA reconstruction (a(m)) of previous scanline
% eta_in => noise of previous scanline
% f <= IAA reconstruction (a(m))
% eta <= noise vector

N=length(x);
eta= eta_in;
diaaf=f_in;


for k=1:q_i
    q=(ifft(diaaf))*K;   
    c=q(1:N)+vt*eta*[1;zeros(N-1,1)]; 
    [A,af]=levinson(c);
    A=A.'/sqrt(af);
    y=tvec_gs_i(A,x);
    diaa_num=fft([y;zeros(K-N,1)]);
    fa1=coeff_gs(A);
    Fa=fft( fa1,K);
    diaa_den=[Fa(1);Fa(end:-1:2)];
    dff=diaa_num./diaa_den;
    diaaf=abs(dff).^2;  
    aa=A; bb=conj(A(end:-1:1)); bb=[0;bb(1:end-1)];
    diag_a=cumsum(abs(aa).^2)-cumsum(abs(bb).^2);
    eta= mean( abs(y./diag_a).^2 );
 end
   
% f=diaaf;
f=dff;
   
% MAP
 q=(ifft(diaaf))*K;   
    c=q(1:N)+vt*eta*[1;zeros(N-1,1)]; 
    [A,af]=levinson(c);
    A=A.'/sqrt(af);
    y=tvec_gs_i(A,x);
    diaa_num=fft([y;zeros(K-N,1)]);
    fmap=diaaf.*diaa_num;


function [y,q]=tvec_gs_i(a,x)
%
% inverse Toeplitz times vector
% using GS form  
%
m=length(a);

z=a;
t1=[z;z(1);complex(zeros(m-1,1))];
t1t=[conj(z(1));complex(zeros(m-1,1));conj(z(1));conj(z(end:-1:2))];

w=[x;complex(zeros(m,1))];
tmp=ifft(fft(t1t).*fft(w));
y11=tmp(1:m);
w=[y11;complex(zeros(m,1))];
tmp=ifft(fft(t1).*fft(w));
y1=tmp(1:m); 
z=[complex(0);conj(a(end:-1:2))];
t1=[z;z(1);complex(zeros(m-1,1))];
t1t=[conj(z(1));complex(zeros(m-1,1));conj(z(1));conj(z(end:-1:2))];
w=[x;complex(zeros(m,1))];
tmp=ifft(fft(t1t).*fft(w));
y12=tmp(1:m);
w=[y12;complex(zeros(m,1))];
tmp=ifft(fft(t1).*fft(w));
y2=tmp(1:m);
y=y1-y2;


function phi=coeff_gs(a)
%
% Estimation of the coefficients of the trigonomertic polynomial
%
m=length(a);
M=[1:m]';
t1=a;
t2=[complex(0);conj(a(end:-1:2))];
s1=conj(t1(end:-1:1)).*M;
s2=conj(t2(end:-1:1)).*M;
w1=[s1;complex(zeros(m,1))];
t1x=[t1;t1(1);complex(zeros(m-1,1))];
tmp=ifft(fft(t1x).*fft(w1));
f1=tmp(1:m); 
w2=[s2;complex(zeros(m,1))];
t2x=[t2;t2(1);complex(zeros(m-1,1))];
tmp=ifft(fft(t2x).*fft(w2));
f2=tmp(1:m); 
f=f1-f2;
phi=[conj(f);f(end-1:-1:1)];
    