	include "sincos.asm"
	include "fftr2a_x.asm"
	org lh:
COEFFS	dsm	ANALYSIS_LEN
	sincos 	ANALYSIS_LEN,COEFFS
	org p:
ProcessData
	ejsr
;	move	#67,r0
;	move	r0,r1
;	jsr	SCI_Snd
;	erts
;First transfer the data from the CDC_BUF to the analysis buffer
;At the same time we calculate the mean
	clr	a x:CDC_RCV_PTR,r0
	clr 	b r0,r4
	move	#(CDC_BUF_LEN-1),m0
	move	m0,m4
	move	#ANALYSIS_BUF,r1
	move	#-1,m1
	move	r1,r5
	move	m1,m5
  	move	x:(r0)+,x0 y:(r4)+,y0
	do	#ANALYSIS_LEN,PrDat_loop1
	  add	  x0,a x0,x:(r1)+ y0,y:(r5)+
	  add	  y0,b x:(r0)+,x0 y:(r4)+,y0
PrDat_loop1
	move	(r0)- ; Restore address of the last sample
	move	r0,x:CDC_RCV_PTR
Test1	nop
;We calculate the mean
	rep	#LOG2_AN_LEN
	asr	b
	rep	#LOG2_AN_LEN
	asr	a
;We subtract the mean
	move	#ANALYSIS_BUF,r0
	move	#-1,m0
	move	a,x1
	move	b,y1
	do	#ANALYSIS_LEN,PrDat_loop2
	  move	  x:(r0),a
	  sub	  x1,a y:(r0),b
	  sub	  y1,b a,x:(r0)
	  move	  b,y:(r0)+
PrDat_loop2
;Now we calculate the FFT
        fftr2a	ANALYSIS_LEN,ANALYSIS_BUF,COEFFS
;Then we calculate the fmax using the simplest 'percentile' algorithm
;The results of the FFT are now sorted in the reverse bits order
;To simplify the further analysis we descramble them
;using the temporary buffer
	move	#0,m0
	move	#(ANALYSIS_LEN/2),n0
	move	#ANALYSIS_BUF,r0
	move	#-1,m1
	move	#TMP_AN_BUF,r1
	do	#ANALYSIS_LEN,Unscramble_loop1
	  move	  l:(r0)+n0,x
	  move	  x,l:(r1)+
Unscramble_loop1
;Because access to the analysis buffer is faster, we copy the unscrambled
;data back to the analysis buffer 
	move	#-1,m0
	move	#ANALYSIS_BUF,r0
	move	#-1,m1
	move	#TMP_AN_BUF,r1
	do	#ANALYSIS_LEN,Unscramble_loop2
	  move	  l:(r1)+,x
	  move	  x,l:(r0)+
Unscramble_loop2
;We send the marker of new data record
	move	#$80,r1
	jsr	SCI_Snd
;Now we analyse the flow toward the transducer (f>0)
;First we calculate the power of the signal
	move 	#-1,m0
	move	#ANALYSIS_BUF+1,r0
	move 	m0,m4
	move	r0,r4
	clr	a
	clr	b
	do	#(ANALYSIS_LEN/2)-1,PrDat_loop3
	  clr	  b x:(r0),x0 y:(r4)+,y0
	  mac	  x0,x0,b
	  mac	  y0,y0,b
	  add	  b,a b,l:(r0)+
PrDat_loop3
;In the analysis buffer we have now square of the spectrum
;Accumulator "a" contains the sum of the power (Ptot)	
;Now we calculate the 3/4 of the power and place it in the "b"
	move	a,b
	rep	#2
	asr	a
	sub	a,b
;Now we calculate the sum of the power once again, until the
;sum exceedes 3/4 Ptot
	move	#ANALYSIS_BUF+1,r0
	move	#-1,m0
	clr	a
	do	#(ANALYSIS_LEN/2)-1,PrDat_loop4
	  move	  l:(r0)+,x
	  move	  lc,r3
	  add	  x,a
	  cmp	  a,b
	  jge	  PrDat_l4_cont
	  enddo
	  nop
PrDat_l4_cont
	nop
PrDat_loop4
;R3 contains the maximum frequency (in fact ANALYSIS_LEN/2-Fmax)
	move 	r3,r1
	jsr	SCI_Snd ; Send the Fmax through the SCI
	move	r3,a
	rep	#$0f
	asl	a
	move	a1,x:SSI_TX_BUF ; Send the Fmax to the analog output
;Now we calculate the mean frequency...
;In the accumulator "a" we calculate the sum of the power (P)
;In the accumulator "b" we calculate the sum of P*f
	clr	a
	clr	b
	move	#1,r3 ; Frequency
	move	#-1,m3
	move	#(ANALYSIS_BUF+1),r0
	move	m3,m0
	do	#(ANALYSIS_LEN/2)-1,EndLoopPrDat5
	  move	  r3,x1
	  move	  x:(r0)+,y1
	  add	  y1,a
	  move	  (r3)+
	  mac	  x1,y1,b
EndLoopPrDat5	  
	jsr	CalcMeanFreq
	move	a1,r1
	jsr	SCI_Snd ;We send the mean frequency
;We send the logarithm of the power
	move	#ANALYSIS_BUF+1,r2
	move	#-1,m2
	do	#(ANALYSIS_LEN/2)-1,SendPower_l1
	  move	  l:(r2)+,a
	  move	  #0,r1
	  rep	  #$2f
	  norm	  r1,a
	  jsr	  SCI_Snd
SendPower_l1
	nop		  
;Now we analyse the flow outward the transducer (f<0)
;The algorithm is identical with the one used for analysis
;of flow toward the transducer. Only the loop limits are different.
;First we calculate the power of the signal
	move 	#-1,m0
	move	#ANALYSIS_BUF+ANALYSIS_LEN-1,r0
	move 	m0,m4
	move	r0,r4
	clr	a
	clr	b
	do	#(ANALYSIS_LEN/2)-1,PrDat_loop3b
	  clr	  b x:(r0),x0 y:(r4)-,y0
	  mac	  x0,x0,b
	  mac	  y0,y0,b
	  add	  b,a b,l:(r0)-
PrDat_loop3b
;In the analysis buffer we have now square of the spectrum
;Accumulator "a" contains the sum of the power (Ptot)
;Now we calculate the 3/4 of the power and place it in the "b"
	move	a,b
	rep	#2
	asr	a
	sub	a,b
;Now we calculate the sum of the power once again, until the
;sum exceedes 3/4 Ptot
 	move	#ANALYSIS_BUF+ANALYSIS_LEN-1,r0
	move	#-1,m0
	clr	a
	do	#(ANALYSIS_LEN/2)-1,PrDat_loop4b
	move	l:(r0)-,x
	move	lc,r3
	add	x,a
	cmp	a,b
	jge	PrDat_l4b_cont
	enddo
	nop
PrDat_l4b_cont
	nop
PrDat_loop4b
;R3 contains the maximum frequency (in fact ANALYSIS_LEN/2-Fmax)
	move 	r3,r1
	jsr	SCI_Snd ; Send the Fmax through the SCI
        move    r3,a
        rep     #$0f
        asl     a
        move    a1,x:SSI_TX_BUF+1 ; Send the Fmax to the analog output
;Now we calculate the mean frequency...
;In the accumulator "a" we calculate the sum of the power (P)
;In the accumulator "b" we calculate the sum of P*f	
	clr	a
	clr	b
	move	#1,r3 ; Frequency
	move	#-1,m3
	move	#(ANALYSIS_BUF+ANALYSIS_LEN-1),r0
	move	m3,m0
	do	#(ANALYSIS_LEN/2)-1,EndLoopPrDat5b
	  move	  r3,x1
	  move	  x:(r0)-,y1
	  add	  y1,a
	  move	  (r3)+
	  mac	  x1,y1,b
EndLoopPrDat5b
	jsr	CalcMeanFreq
	move	a1,r1
	jsr	SCI_Snd
;We send the logarithm of the power
	move	#ANALYSIS_BUF+ANALYSIS_LEN-1,r2
	move	#-1,m2
	do	#(ANALYSIS_LEN/2),SendPower_l1b
	  move	  l:(r2)-,a
	  move	  #0,r1
	  rep	  #$2f
	  norm	  r1,a
	  jsr	  SCI_Snd
SendPower_l1b
	nop		  
	erts
	
CalcMeanFreq
;This subroutine performs the b/a division, assuring the proper scaling
;INPUT: a,b  >0 !!!
;USES: r0,r3,m0,m3
;OUTPUT: A1.A0 - result
;First we normalize the values
	ejsr
	move	#0,r0
	move	#0,r3
	rep	#$2f
	norm	r0,a
	rep	#$2f
	norm	r3,b
	asr	b	;To ensure a > b
;Finally we divide the sum(f*p)/sum(p)
	move	a,x1
	and	#$fe,CCR
	rep	#$18
	div	x1,b
	add	x1,b
;B0 contains the result of division
	clr	a
	move	b0,a1
;Now we consider the results of normalisation
	move	r0,b
	move	r3,x1
	sub	x1,b
	dec	b
	rep	b1
	asr	a
	erts