SDeconv(hin)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function SDeconv is self deconvolution.
% gout = SDeconv(hin)
% input: hin length of hin must be 2*L +1, where L is the length of g
% output: gout the output with length L
% 
% By Xiaozhou Huang
% March 9th, 2004
% Note: give h=g*g, find out g
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function gout = SDeconv(hin)

if rem(length(hin),2) ~= 1
    error('length of the input must be odd');
end
L = (length(hin)-1)/2;    % L+1 is the length of channel

gout=zeros(1,L+1);

if (abs(hin(1))>abs(hin(length(hin))))   % forward elimination

gout(1) = sqrt(hin(1));

for ii=1:L
    temp=0;
    for jj=1:ii-1
        temp=temp+gout(jj+1)*gout(ii-jj+1);
    end
    gout(ii+1) = (hin(ii+1)-temp)/(2*gout(1));
end
else                                    % backward substitution
    gout(L+1)=sqrt(hin(length(hin)));
    for ii=L-1:-1:0
        temp=0;
        for jj=ii+1:L-1
            temp=temp+gout(jj+1)*gout(ii+L-jj+1);
        end
        
        gout(ii+1)=(hin(ii+1+L)-temp)/(2*gout(L+1));
    end
    
end