/*
	This C language subroutine computes the
	autocoherence sequence of its floating
	point input.

	The Run Time Library for the C Language.

	Gordon A. Sterling	(617) 461 - 3076
	DSP Development Tools Engineering
    
	Created on 5/27/90	

	#include <stats.h>			<Header>
	float autocoh(float out[], float in[],\ <Prototype>
		      int samples, int lags);	<Prototype>

	
*/

#include "lib_glob.h"
#include "sts_glob.h"

.SEGMENT/CODE	Code_Space_Name;

.FILE 		RTL_FILENAME;

.GLOBAL	    	_autocoh;

_autocoh:	R1=PASS R8, FETCH_RETURN
		R9=PASS R12, put(R9);
		F12=FLOAT R9, put(R7);

		F7=1.0;
		F8=F7+F7, readparam4(R2); /*F8=2.0, rd 1st stk param	*/

		F0=RECIPS F12, R2=dma_i; /*Get 4 bit seed R0=1/D	*/		
		F12=F0*F12, put(R2); 	/*D'=D*R0, push dma_i		*/

		F7=F0*F7, F0=F8-F12;	/*F0=R1=2-D', F7=N*R0		*/
		F12=F0*F12, dm_lnt=dm_0; /*F12=D'=D'*R1, zero out dm_lnt*/
		F7=F0*F7, F0=F8-F12;	/*F7=N*R0*R1, F0=R2=2-D'	*/
		F12=F0*F12, R2=dmb_i; 	/*F12=D'=D'*R2			*/
		F7=F0*F7, F0=F8-F12;	/*F7=N*R0*R1*R2, F0=R3=2-D'	*/
		F7=F0*F7, R8=R1;	/*F7=N*R0*R1*R2*R3, restore r8	*/

		IF LE JUMP (PC, error_return) (DB); /*Must have > 1 point*/
		R1=R1-R1, dm_ptr=R8;    /* pt to slot			*/
		R12=R9-1, F0=DM(dm_ptr,1);
		
		LCNTR=R12, DO sum UNTIL LCE;
sum:		    F1=F1+F0, F0=DM(dm_ptr,1);
		F1=F1+F0, put(R2);	/*push dmb_i			*/
		F1=F1*F7, dm_bse=R4;	/*dm_ptr implicitly gets r4 also*/
		F1=F1*F1, R4=R8;

		LCNTR=R2, DO lag_loop UNTIL LCE;
		    R9=R9-1, dma_i=R4;
		    R0=R0-R0, dmb_i=R8;
		    F2=DM(dmb_i,1);
		    F3=DM(dma_i,1);
	    	    LCNTR=R9, DO element_loop UNTIL LCE;
			F3=F3*F2, F2=DM(dmb_i,1);
element_loop:  		F0=F0+F3, F3=DM(dma_i,1);
		    F3=F3*F2;
		    F0=F0+F3;
		    F0=F0*F7;
		    F0=F0-F1, R3=dm_bse; /*Restore orig r4 for return	*/
lag_loop:	    R4=R4+1, DM(dm_ptr,1)=F0;
		
restore_state:	get(dmb_i,1);
		get(dma_i,2);
		R0=PASS R3, get(R7,3);
		get(R9,4);
		RETURN (DB);
		RESTORE_STACK
		RESTORE_FRAME

error_return:	JUMP (PC, restore_state) (DB);
		R0=0;
		nop;

.ENDSEG;