function [Tpoints codes problem]=SFRAC(m,comp_offset,opamp_error,cap_mismatch,margin)
%Sharif Fast Residual Based ADC Charactizer
%Finds the transition points of a pipelined A2D by backward derivation algorithm
%All Inputs And Outputs Must Be Normalized to Fullscale Input : (-1, 1) 
%%%%%%%%%%%Designed By Muhsin Hassanpour Ghadi%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%Date : 07, April, 2013

% Inputs:
%m:							Vector(1,number of stages): Bit per stage for eash stage of pipelined ADC
%comp_offset				Matrix(2^maximum(bit per stage)-2,number of stages): Comparator Offset for each stage of pipelined ADC
%opamp_error				Vector(1,number of stages): Error due to opamp non-idealities
%cap_mismatch				Vector(1,number of stages): Capacitor Mismatch for each stage
%margin 					Number(1): Possible offset error of ADC

% Outputs:
%tpoints                    Vector(1, number of Tpoints): Stores transition points(tpoints)
%codes                      Vector(1, number of Tpoints): The digital code after tpoint
%problem                    Number(1): Returns the stage that misses a code

% Description:
% Backwardly finds the Tpoints of each stage of pipelined ADC,

% Example:
%	m=[2 1 1 1 1 1 1 2];
%	error=[4/2000 2/2000 2/2000 2/2000  2/400 2/400 2/400 0]; 
%   comp_offset=.02*randn(n,2^(max(m)+1)-2)*1;
%   cap_mismatch=0.001*[.34 1 1 1 1.5 1.5 1.5 0].*randn(1,n)*1;	

% 10 bit ADC: first stage 2.5 bit/stage following with 1.5 bit/stage and the last stage is 2bit flash
% *sum(m) is equal to resolution of ADC
% *length(m) is equal to number of stages in Pipelined ADC
% Opamp errors with dc gains in stage 1-4 equal to 2000 and stage 5-7 equal to 400, the last stage has no opamp, settling error =0;
% *In this version opamp error is represented by 'e' in vout=(1-e)*(G*Vin-Vdac)
% *e= 1/(f*Gain) + settling error; in which, opamp is under the feedback equal to f and with the dc gain of Gain, settling error due to finite BandWidth 
% *error vector must be equal to 'm' vector in length
% Comparator offset (1* sigma of guassian distribution) equal to 20mv
% *The first stage is 2.5 bit/stage and MUST have 6 comparator offsets, 
% *The 2-7 stages are 1.5 bit/stage and MUST have 2 comparator offsets,
% *The last stage flash architecture uses 3 comparators.
% *Now * in this vector are considered as comparator offsets:
%	[ * * * * * * * *
%	  * * * * * * * 0
%	  * 0 0 0 0 0 0 *
%	  * 0 0 0 0 0 0 0
%     * 0 0 0 0 0 0 *
% 	  * 0 0 0 0 0 0 0] 
% Capacitor Mismatches are equal to (1*sigma of guassian distribuition) equal to 0.001*(.34 1 1 1 1.5 1.5 1.5 0)
% *The vector must be equal to 'm' vector in length

% INL, DNL Calculation Codes:
%       [Tpoints,codes,prob] = PABDS(m,comp_offset,error,cap_mismatch,90);
%		nbit=sum(m);
%		lengthcode=length(codes);
%       v1=Tpoints(1);
%       vlast=Tpoints(end);
%       LSB=(vlast-v1)/(2^nbit-2);
%       tpoint1=Tpoints(1);
%       codename=codes(1)+1;
%       DNL=zeros(1,2^nbit)-1;
        
%        for counter1=2:lengthcode
%            tpoint2=Tpoints(counter1);
%            DNL(codename)=DNL(codename)+(tpoint2-tpoint1)/LSB;
%            codename=codes(counter1)+1;
%            tpoint1=tpoint2;
%        end
%        INL=zeros(1,2^nbit);
%        previnl=0;
%        for counter2=2:2^nbit-1          
%            previnl=previnl+DNL(counter2);
%            INL(counter2)=previnl;
%        end

% Monte Carlo Simulation:
% for i=1:number_of_samples
%           tline =fgetl(fileid);
%           error=sscanf(tline,'%f')';
%           tline =fgetl(fileid);
%           cap_mismatch=sscanf(tline,'%f')';

%           comp_offset=zeros(n,widthm);
%           for j=1:widthm
%               tline =fgetl(fileid);
%               comp_offset(:,j)=sscanf(tline,'%f');
%           end
%           [tpoint,codes,prob] = real_pipelined_adc_code5(m,comp_offset,error,cap_mismatch,90);
% end


num_stage=length(m);                        %Number of Stages
num_bit=sum(m);                             %Number of Bits

CTP_in_t=zeros(1,2*2^num_bit)+10;           %Main Tpoints
CTP_in_c=zeros(1,2*2^num_bit);              %Main Codes
CTP_out_t=zeros(1,2*2^num_bit)+10;          %Tpoints For Next Stage
CTP_out_c=zeros(1,2*2^num_bit);             %Codes For Next Stage


count_stage=num_stage-1;                    %Number of Stage that is analyzing
real_gain=2.^m+cap_mismatch;                %Real gain of stages
expected_gain=2.^m;                         %Ideal gain of stages
epsilon=1-opamp_error;                      %e=(1-errors)
problem=0;


% Defining codes and T points of the last stage
CTP_out_t(1:(2^m(end)-1))=linspace(-1+5/2^(m(end)+1),1-3/2^(m(end)+1),2^m(end)-1)+comp_offset(num_stage,2*(1:(2^m(end)-1)));
CTP_out_c(1:(2^m(end)-1))=1:(2^m(end)-1);
lastquantizercode=(2^m(end)-1);

powers=expected_gain(num_stage)/2;

% Meging All the stages of the pipelined ADC to get the equal quantizer
% count_stage is the number of stages, for the first merge it equals to n-1
while (count_stage>0)
    
    %Merge Function:
    %Merges an MDAC with its next quantizer
    lastcode=0;
    begincode=0;
    for i=0:(2^(m(count_stage)+1)-2)
        
        % line1 and line2 are two borders of each MDAC interval
        if i>0;line1=-1+3/2^(m(count_stage)+1)+(i-1)/2^m(count_stage)+comp_offset(count_stage,i);else line1=-1-margin;arraycode=0;end
        if i<(2^(m(count_stage)+1)-2);line2=-1+3/2^(m(count_stage)+1)+(i)/2^m(count_stage)+comp_offset(count_stage,i+1);else line2=1+margin;end
        
        % Each border of MDAC sub-doamin must be checked if they make a
        % Tpoint
        % begincode is the beginning code of the sub-domain
        % arraycode is the place of Tpoint in the arrayT
        % i is the code of the MDAC sub-domain
        if arraycode>0
            % Aka Vi+1=(G*Vi-VDAC)*E
            vi=(real_gain(count_stage)*line1-(1+cap_mismatch(count_stage))*(i+1-2^m(count_stage)))*epsilon(count_stage);
            
            % If its the first border 
            if arraycode==1
                
                CTP_in_t(1)=line1;
                CTP_in_c(1)=1;
                lastcode=begincode;
                arraycode=2;
            elseif arraycode>1
%                 begincode=sum(find((vi>CTP_out_t),1,'last'));
                begincode=0;    
                while vi>CTP_out_t(begincode+1)
                    begincode=begincode+1;
                    if begincode>length(CTP_out_t)
                        break;
                    end
                end
                
                if begincode==0
                    newcode=i*powers;
                else
                    
                    newcode=CTP_out_c(begincode)+i*powers;
                end
                if newcode~=lastcode;
                    CTP_in_t(arraycode)=line1;
                    CTP_in_c(arraycode)=newcode;
                    lastcode=newcode;
                    arraycode=arraycode+1;
                end
            end
        else
            arraycode=1;
        end
        
        % Each T point of next quantizer is checked here
        for j=(begincode+1):lastquantizercode
                newcode=CTP_out_c(j)+i*powers;
                vt=CTP_out_t(j)/real_gain(count_stage)/epsilon(count_stage)+(i+1-expected_gain(count_stage))*(1+cap_mismatch(count_stage))/real_gain(count_stage);
                if (vt<line2)
                    CTP_in_t(arraycode)=vt;
                    CTP_in_c(arraycode)=newcode;
                    lastcode=newcode;
                    arraycode=arraycode+1;
                else
                    break
                end
        end
    
    end
    
    lastquantizercode=arraycode-1;
    if count_stage>0;powers=expected_gain(count_stage)*powers;end
    if (problem==0 && lastquantizercode~=(2*powers-1));problem=count_stage;end
    count_stage=count_stage-1;
    CTP_out_t=CTP_in_t;
    CTP_out_c=CTP_in_c;
    CTP_in_t=CTP_in_t-CTP_in_t+10;
    CTP_in_c=CTP_in_c-CTP_in_c;

end

Tpoints=CTP_out_t(1:lastquantizercode);
codes=CTP_out_c(1:lastquantizercode);

end