0
$\begingroup$

I have a fixed point DSP library fft function that does not use a block exponent. I don't have prior knowledge of that implementation but whatever I could read online, my understanding is that its used in each butterfly stage to determine redundant sign bits and hence left shift by same amount.However, in my case, library function does not give me parameter for the block exponent.

1) Is it possible to do same precision improvement as happens with block exponent for FFT , but done on streamed/block of input data that goes into fft?
2) if yes, Can anyone show me how that can be done using small C code, assuming 24 bit fixed point data(4 int,20 frac) and explain the logic behind it ( with clear explanation at bit/binary level,as I couldn't find that anywhere) assuming input data divided into frames/blocks of size 256 for FFT size 256.

Thanks sedy

$\endgroup$
1
$\begingroup$

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
$\endgroup$
0
$\begingroup$

Block (floating point) has been implemented in a Vector Signal Processor developed by Zoran (sold by now) over 30 years ago.

You may get very close results to such implementation by employing a floating point (even double precision) and than scale the result so the maximal value of the FFT uses all 24 bits (usually you use 2's complement representation of the signal).

I do not think that you will find articles that really explain the power of Block Floating Point as means to greatly improve the processing of fixed point hardware.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.