i don't have a C version of fixed-point FFT, but i have an old 56K version. i used to have a 68K version but that's in a Mac that has died decades ago.
if you can get a good fixed-point version of the FFT in C that you like, i can show you where to place the magnitude tests and where to kick it into "divide-by-two" mode.
your magnitude test should confirm that at least the 3 MSBs are identical (that the magnitude does not exceed 1/4 full scale) for every real and imaginary result of an FFT pass for the following pass to be unscaled. if the magnitude of any value exceeds 1/4 full scale, then a sticky bit should be set and the following pass should be in divide-by-two mode.
page 255
opt cex,mex
lstcol 8,8,10,13,13
;
; copyright (C) 1994 by Robert Bristow-Johnson
; tel: 201/429-1509
;
;
;
; This subroutine is identical to the Motorola benchmark fftr2a
; except that scaling is optimized using block floating point.
; This requires a DSP56002, DSP56004, DSP56007 or DSP56001A.
; Some FFT passes have the divide by 2 scaling mode set (depending
; on the result of the S bit in the CCR). This insures
; against, overflow for any signal. N must be a power of 2.
;
; To be completely safe, the input should be scaled so that
; the magnitude of every real or imaginary part does not
; exceed 0.31066 ($27C3B5) .
;
;
; N-1
; X[n] = 1/M SUM{ x[k]exp(-j2 PI nk/N) }
; k=0
;
; where
; x[k] = input in normal order
; X[n] = output in bit reverse order
;
; N = 2^p
;
; M = 2^q (q = integer, number of passes that
; right shifting was necessary.
; 0 <= q <= p )
;
;
; input: r1 -> data; x:(r1) is real part, y:(r1) is imag part
; r2 -> twiddle coeffients; x:(r2) is cos, y:(r2) is sin
; n1 = log2(number of points) = p
; r4 = block floating point exponent input
;
; output: r4 = block floating point exponent output
; = q + r4 (input)
;
; uses: a, b, x0, x1, y0, y1, mr, ccr
; r1, r2, r3, r4, r5, r6
; n1, n2, n3, n5, n6
; m1, m2, m3, m4, m5, m6
; x:<data
; x:<coef
; x:<n_grp
;
; Leaves mr in no scaling mode.
;
; Requires 78 words of internal p: memory
; and p*( 3*N + 20 ) + 8*N + 15 instruction cycles.
;
; To execute optimally, the fft routine sould be in
; internal p: space.
;
;
org x:$0
data ds $1
n_grp ds $1
coef ds $1
org p:$40
fft
clr a r1,x:<data ;1
move #<1,a1 ;1
move a1,x:<n_grp ;1 init groups per pass
do n1,_shift ;3
asl a ;p
_shift
lsr a r2,x:<coef ;1
lsr a a1,n5 ;1 init bfly per group (and offset between bfly top and bfly bottom addresses)
move a1,n3 ;1 init twiddle pointer offset
move #>-1,m5 ;2 init bfly top and bfly bottom address modifiers
move m5,m6 ;1 for linear addressing
move m5,m4 ;1
move m5,m1 ;1
move m5,m2 ;1
move #<0,m3 ;1 init twiddle address modifier for bit-reversed addressing
andi #$7F,ccr ;1 clear sticky bit
;
; Do FFT passes with triple nested DO loop.
;
do n1,end_pass ;3
move x:<data,r5 ;p init bfly top input pointer
move r5,r1 ;p init bfly top output pointer
lua (r5)+n5,r6 ;p init bfly bottom input pointer
move x:<coef,r3 ;p init twiddle pointer
lua (r6)-,r2 ;2p init bfly bottom output pointer
move n5,n6 ;p init bfly pointer offsets
move n5,n1 ;p
move n5,n2 ;p
jset #7,sr,scale_down ;2p test scaling bit
andi #$F3,mr ;p reset scaling to normal scaling mode
do x:<n_grp,_end_group1 ;3p
move x:(r6),x1 y:(r3),y0 ;N-1 preload x1, lookup -sin value
move x:(r2),a y:(r5),b ;N-1 preload a and b
move x:(r3)+n3,x0 ;N-1 lookup -cos value, update twiddle pointer
do n5,_end_bfly1 ;3(N-1)
mac x1,y0,b y:(r6)+,y1 ;pN/2 Im[top] - sin*Re[bot]
macr -x0,y1,b a,x:(r2)+ y:(r5),a ;pN/2 Im[top] - sin*Re[bot] + cos*Im[bot] = Im[ top + twiddle*bot ]
subl b,a x:(r5),b b,y:(r1) ;pN/2 2*Im[top] - (Im[top] - sin*Re[bot] + cos*Im[bot]) = Im[top] + sin*Re[bot] - cos*Im[bot] = Im[ top - twiddle*bot ]
mac -x1,x0,b x:(r5)+,a a,y:(r2) ;pN/2 Re[top] + cos*Re[bot]
macr -y1,y0,b x:(r6),x1 ;pN/2 Re[top] + cos*Re[bot] + sin*Im[bot] = Re[ top + twiddle*bot ]
subl b,a b,x:(r1)+ y:(r5),b ;pN/2 2*Re[top] - (Re[top] + cos*Re[bot] + sin*Im[bot]) = Re[top] - cos*Re[bot] - sin*Im[bot] = Re[ top - twiddle*bot ]
_end_bfly1
move a,x:(r2)+n2 y:(r6)+n6,y1 ;N-1 update bfly top and bfly bottom pointers
move x:(r5)+n5,x1 y:(r1)+n1,y1 ;N-1
_end_group1
jmp <pass_common ;2p
scale_down
move (r4)+ ;(p) increment r4 to indicate one pass of scaling down
andi #$7F,ccr ;(p) clear sticky bit
ori #$04,mr ;(p) set scaling mode to divide by 2
do x:<n_grp,_end_group2 ;(3p)
move x:(r6),x1 y:(r3),y0 ;(N-1) preload x1, lookup -sin value
move x:(r2),a y:(r5),b ;(N-1) preload a and b
asl a x:(r3)+n3,x0 ;(N-1) double a to compensate for div by 2 scaling, lookup -cos value, update twiddle pointer
do n5,_end_bfly2 ;(3(N-1))
mac x1,y0,b y:(r6)+,y1 ;(pN/2)
macr -x0,y1,b a,x:(r2)+ y:(r5),a ;(pN/2)
subl b,a x:(r5),b b,y:(r1) ;(pN/2)
mac -x1,x0,b x:(r5)+,a a,y:(r2) ;(pN/2)
macr -y1,y0,b x:(r6),x1 ;(pN/2)
subl b,a b,x:(r1)+ y:(r5),b ;(pN/2)
_end_bfly2
move a,x:(r2)+n2 y:(r6)+n6,y1 ;(N-1) update bfly top and bfly bottom pointers
move x:(r5)+n5,x1 y:(r1)+n1,y1 ;(N-1)
_end_group2
pass_common
move n5,b1 ;p divide bfly per group by 2
lsr b x:
;
;pi equ 3.141592654
;freq equ 2.0*pi/@cvf(points)
;
; org x:coef
;count set 0
; dup points/2
; dc -@cos(@cvf(count)*freq)
;count set count+1
; endm
;
; org y:coef
;count set 0
; dup points/2
; dc -@sin(@cvf(count)*freq)
;count set count+1
; endm
;
; endm ;end of sincos macro
;
;
; this generates the FFT twiddle coefficients in normal order.
; fft expects that and uses bit reversed addressing to fetch the twiddle coefs.
;
; the values are negated because -1.0 can be exactly represented but +1.0 cannot.
; fft expects that and negates the multiplications in the butterflies.
;
sincos
andi #$F3,mr
move #>-1,m5
move m5,m2
clr a #delta_table,r5
move #<2,n5
move #<1,a1
do n1,_shift
asl a (r5)+n5
_shift
lsr a #<-1.0,x0 ;-cos
move a,n5 ;N/2
move p:(r5)+,x1 ;-cosdel
move p:(r5)+,y1 ;-sindel
mpyr -x0,x1,a #<0.0,y0 ;-sin
do n5,_sincos_loop
mpy -x0,y1,b a,x0
macr -y0,x1,b a,x:(r2) ; b = new sin
mpy -x0,x1,a b,y0
macr y0,y1,a b,y:(r2)+ ; a = new cos
_sincos_loop
rts
twopi equ 6.28318530717959
delta_table
dc -@cos(twopi/1.0)
dc -@sin(twopi/1.0)
dc -@cos(twopi/2.0)
dc -@sin(twopi/2.0)
dc -@cos(twopi/4.0)
dc -@sin(twopi/4.0)
dc -@cos(twopi/8.0)
dc -@sin(twopi/8.0)
dc -@cos(twopi/16.0)
dc -@sin(twopi/16.0)
dc -@cos(twopi/32.0)
dc -@sin(twopi/32.0)
dc -@cos(twopi/64.0)
dc -@sin(twopi/64.0)
dc -@cos(twopi/128.0)
dc -@sin(twopi/128.0)
dc -@cos(twopi/256.0)
dc -@sin(twopi/256.0)
dc -@cos(twopi/512.0)
dc -@sin(twopi/512.0)
dc -@cos(twopi/1024.0)
dc -@sin(twopi/1024.0)
dc -@cos(twopi/2048.0)
dc -@sin(twopi/2048.0)
dc -@cos(twopi/4096.0)
dc -@sin(twopi/4096.0)
dc -@cos(twopi/8192.0)
dc -@sin(twopi/8192.0)
dc -@cos(twopi/16384.0)
dc -@sin(twopi/16384.0)
dc -@cos(twopi/32768.0)
dc -@sin(twopi/32768.0)
END
