      (* ********************************************************** *)
      (*							    *)
      (*			  FFT				    *)
      (*		       APPLICATION			    *)
      (*							    *)
      (*  Fast Fourier math package.  Marcel Hendrix, 26 April 1987 *)
      (*       Joe Barnhart, Dr Dobb's Journal, September 1984	    *)
      (*  	   December 13 1989, made it 20% faster		    *)
      (*   This app shows I/O to FFT/IFFT and builds a synthesizer  *)
      (*							    *)
      (* ********************************************************** *)

	?DEF -graphics 0= IFTRUE LOAD \drawing\drivers\graphics.lib
				 TRUE TO ?AutoCr TRUE TO ?AutoWindow
				 3270  8 TO FontWidth
			  IFEND

			ONLY FORTH DEFINITIONS
				DECIMAL

	     ?DEF -fft 0= IFTRUE LOAD \digitize\fft
			  IFEND


?REM -synth


: .HELP
CR ." -----------------------------------------------------------------------"
CR ."		   Example Fast Fourier Transform Demo" 
CR ."		      Manipulates a complex SineWave" 
CR ." -----------------------------------------------------------------------"
CR
CR ." <spec> FFT		-- convert samples to spectrum."
CR ." <spec> IFFT		-- convert spectrum to samples."
CR ." 0Table	    ( Wave )	-- clears Wave table"
CR ." <n> !DC	    ( Wave )	-- sets DC level of spectrum to <n>."
CR ." <n> <x> !R#   ( Wave )	-- Real part of complex amplitude[x] <= <n>." 
CR ." <n> <x> !I#   ( Wave )	-- Imaginary part of complex amplitude[x] <= <n>." 
CR ." Wave>Data			-- copies Wave table to Data."
CR ." .FFT			-- do FFT on Data and show output."
CR ." .IFFT			-- do IFFT on Data and show output."
CR ." <spec> .FOURIER		-- show if data were spectrum (frequencies)"
CR ." <spec> .TIMEDOMAIN	-- show if data were samples (time domain)" 
CR ." <spec> .TABLE		-- lists values in either table." 
CR ." .DIFF%			-- check out relative difference after FFT+IFFT" CR
CR ." Note1: <spec> means any of the following XVECTOR's: Wave or Data"
CR ." Note2: !DC !R# and !I# work on the Wave table exclusively."
CR ." Note3: !SINE !SQUARE !PULSE !# !SAW : five pre-defined sessions." 
CR ."        These fill Wave (for user modification) then do Wave>Data etc."
CR ;

				CLS .HELP 

-- Tools for graphic output

: LARGESTMOD?				\ <descriptor> --- <max modulus>
		0 >S 
		DUP /ELEMENTS 1+
		1 ?DO 
			I OVER X@V |X| 
			S> UMAX >S
		 LOOP   DROP 		\ descriptor can go
		S> ;

: LARGESTREAL?				\ <descriptor> --- <|max real|>
		0 >S 
		DUP /ELEMENTS 1+
		1 ?DO 
			I OVER X@V NIP ABS 
			S> UMAX >S
		 LOOP   DROP 		\ descriptor can go
		S> ;

0. DVALUE xsteps
  0 VALUE rectwidth

: FRAME		Color DUP 8 XOR PosblColors CIRCULAR 
		?DUP IF TO Color ENDIF
		0 0 MOVE  Xmax      Ymax     #RECTANGLE 
		2 2 MOVE  Xmax 4 -  Ymax 4 - #RECTANGLE 
		TO Color ;


	-- Fourier spectrum always shown in dB (logarithmic Y-axis)

: .FOURIER	DUP LARGESTMOD? DUP	\ <descriptor> --- <>
		  0= IF 2DROP ." All zeroes." EXIT
		  ENDIF LOG(1+X) >S
		Xmax OVER /ELEMENTS   TO xsteps		\ dvalue.. trick..
		Xmax OVER /ELEMENTS / TO rectwidth
		GCLEAR FRAME
		Hres Vres !XY-RANGE  Or! TO Pen
		DUP /ELEMENTS 
		0 ?DO   I xsteps */  0  MOVE
		        I 1+ OVER X@V |X| LOG(1+X)
			 Ymax 8 - S */  
			rectwidth SWAP #RECTANGLE
		 LOOP DROP -S ;		\ descriptor and maximum


	-- Time domain data is assumed real, and can be negative.

: .TIMEDOMAIN	DUP LARGESTREAL? >S 	\ <descriptor> --- <>
		S 0= IF DROP -S ." All zeroes." EXIT
		  ENDIF 
		Xmax OVER /ELEMENTS   TO xsteps		\ dvalue.. trick..
		Xmax OVER /ELEMENTS / TO rectwidth
		GCLEAR FRAME
		Hres Vres !XY-RANGE  Or! TO Pen
		DUP /ELEMENTS 
		0 ?DO   I xsteps */  Ymax 2/  MOVE
		        I 1+ OVER X@V NIP
			 Ymax 2/ S */  
			rectwidth SWAP #RECTANGLE
		 LOOP DROP -S ;		\ descriptor and maximum
		
: /,		6 ASHR , ;		\ <n> --- <>  0 < entry < 255
: X,		, 0 , ;			\ <Im> --- <>

DECIMAL

256 16384 XVECTOR Data

256 16384 XVECTOR SineWave -256 2 CELLS * ALLOT ( vector initialized to a sine)

  8192 X,  8334 X,  8477 X,  8763 X,  8905 X,  9190 X,  9332 X,  9473 X,
  9755 X,  9895 X, 10173 X, 10312 X, 10450 X, 10723 X, 10859 X, 11127 X,
 11260 X, 11392 X, 11654 X, 11783 X, 12037 X, 12163 X, 12288 X, 12533 X,
 12653 X, 12890 X, 13007 X, 13122 X, 13347 X, 13457 X, 13673 X, 13778 X,
 13984 X, 14084 X, 14183 X, 14374 X, 14467 X, 14647 X, 14734 X, 14819 X,
 14983 X, 15062 X, 15213 X, 15286 X, 15356 X, 15491 X, 15554 X, 15675 X,
 15732 X, 15787 X, 15889 X, 15937 X, 16026 X, 16066 X, 16104 X, 16174 X,
 16204 X, 16259 X, 16283 X, 16304 X, 16339 X, 16352 X, 16372 X, 16379 X,
 16383 X, 16379 X, 16372 X, 16352 X, 16339 X, 16304 X, 16283 X, 16259 X,
 16204 X, 16174 X, 16104 X, 16066 X, 16026 X, 15937 X, 15889 X, 15787 X,
 15732 X, 15675 X, 15554 X, 15491 X, 15356 X, 15286 X, 15213 X, 15062 X,
 14983 X, 14819 X, 14734 X, 14647 X, 14467 X, 14374 X, 14183 X, 14084 X,
 13984 X, 13778 X, 13673 X, 13457 X, 13347 X, 13122 X, 13007 X, 12890 X,
 12653 X, 12533 X, 12288 X, 12163 X, 12037 X, 11783 X, 11654 X, 11392 X,
 11260 X, 11127 X, 10859 X, 10723 X, 10450 X, 10312 X, 10173 X,  9895 X,
  9755 X,  9473 X,  9332 X,  9190 X,  8905 X,  8763 X,  8477 X,  8334 X,
  8191 X,  8049 X,  7906 X,  7620 X,  7478 X,  7193 X,  7051 X,  6910 X,
  6628 X,  6488 X,  6210 X,  6071 X,  5933 X,  5660 X,  5524 X,  5256 X,
  5123 X,  4991 X,  4729 X,  4600 X,  4346 X,  4220 X,  4095 X,  3850 X,
  3730 X,  3493 X,  3376 X,  3261 X,  3036 X,  2926 X,  2710 X,  2605 X,
  2399 X,  2299 X,  2200 X,  2009 X,  1916 X,  1736 X,  1649 X,  1564 X,
  1400 X,  1321 X,  1170 X,  1097 X,  1027 X,   892 X,   829 X,   708 X,
   651 X,   596 X,   494 X,   446 X,   357 X,   317 X,   279 X,   209 X,
   179 X,   124 X,   100 X,    79 X,    44 X,    31 X,    11 X,     4 X,
     0 X,     4 X,    11 X,    31 X,    44 X,    79 X,   100 X,   124 X,
   179 X,   209 X,   279 X,   317 X,   357 X,   446 X,   494 X,   596 X,
   651 X,   708 X,   829 X,   892 X,  1027 X,  1097 X,  1170 X,  1321 X,
  1400 X,  1564 X,  1649 X,  1736 X,  1916 X,  2009 X,  2200 X,  2299 X,
  2399 X,  2605 X,  2710 X,  2926 X,  3036 X,  3261 X,  3376 X,  3493 X,
  3730 X,  3850 X,  4095 X,  4220 X,  4346 X,  4600 X,  4729 X,  4991 X,
  5123 X,  5256 X,  5524 X,  5660 X,  5933 X,  6071 X,  6210 X,  6488 X,
  6628 X,  6910 X,  7051 X,  7193 X,  7478 X,  7620 X,  7906 X,  8049 X,


: .FFT	Data FFT
	Data .FOURIER ;

: .IFFT Data IFFT
	Data .TIMEDOMAIN ;

: .ELAPSED					\ <d> ---
	?MS 2SWAP D- 2DUP
	<# # # #  '.' HOLD #S #>
	TYPE ."	 seconds used --> "
	D2/
	<# # # '.' HOLD #S #>
	TYPE ."	 ms / 256 point FFT." ;

: .SPEED
	." Testing..." CR ?MS
	#100 0 DO
		 Data FFT
		 Data IFFT
	     LOOP
	.ELAPSED ;

: PRINT#	S->D DUP 0< IF '-' 		\ <n> --- <>
			  ELSE BL 
			 ENDIF >S
		DABS <# # # # # # S> HOLD #> TYPE 
		2 SPACES ;

: .TABLE	DUP /ELEMENTS 1+		\ <descriptor> --- <>
		1 ?DO WAIT? 0=
		WHILE I OVER X@V PRINT# PRINT# 
		 LOOP DROP ;			\ prints 5 complex #s/line

#256 #16384 XVECTOR Wave

: 0Table	Wave /ELEMENTS 1+
		 1 DO 
		      0 0  I Wave X!V 
	 	 LOOP ;

-- The numbering of indices can be QUITE confusing:

-- 1  2  ....... 128   129   130 .......  256	<- element of XVECTOR
-- dc fund. ...  127th Nyq.  127'  ..     fund' <- frequencies
-- 0  1     ...  127   <------ automatic -----> <- used in R!# & I!#


: !DC		0 SWAP 1 Wave X!V ;		\ set dc component.

: !R#		#127 MIN >S			\ <a> <h#> --- <>  1 < h# < 127
		S 1+ Wave X@V DROP SWAP		\ throw away Re, keep Im
		XDUP S 1+ Wave X!V		\ cosine term => sine amp=0
		X' #257 S> - Wave X!V ;

: !I#		#127 MIN >S			\ <a> <h#> --- <>  1 < h# < 127
		S 1+ Wave X@V UNDER		\ throw away Im, keep Re
		XDUP S 1+ Wave X!V  		\ sine => complex conjugate
		X' #257 S> - Wave X!V ;

: Wave>Data	Wave Data  
		Data &SIZE  4 CELLS +  CMOVE ;

: !Pt		Data /ELEMENTS MIN 0 MAX 1+	\ <a> <'time'> --- <>
		0 -ROT Wave X!V ;		\ 0 < time < 255

  #128 VALUE Zero
#16000 VALUE Amp
  #127 VALUE MaxH

: +PAGE		Page2See 1+ >GPAGE ;
: =PAGE		Page2Write GPAGE> ;

: DO-IT		Wave>Data 
		+PAGE Data .TIMEDOMAIN  =PAGE KEY DROP
		+PAGE .FFT   =PAGE KEY DROP
		+PAGE .IFFT  =PAGE HOME ;

: >DO-IT<	Wave>Data 
		+PAGE Data .FOURIER  =PAGE KEY DROP
		+PAGE .IFFT =PAGE KEY DROP
		+PAGE .FFT  =PAGE HOME ;

: !#		0Table
		MaxH 1+ 1 DO   Amp I /  I !I#
		     2 +LOOP
		>DO-IT< ;

: !PULSE	0Table 
		Zero 0 DO Amp I !Pt 
		     LOOP 
		DO-IT ;

: !SQUARE	0Table 
		Data /ELEMENTS 
		0 DO Amp
		     I Zero < IF NEGATE 
			   ENDIF 
		     I !Pt 
		LOOP 
		DO-IT ;

: !SAW		0Table 
		Data /ELEMENTS DUP 2/ >S 
		0 DO I S - Amp S */  I !Pt 
		LOOP -S
		DO-IT ;

: Sine>Wave	SineWave Wave  
		Wave &SIZE  4 CELLS +  CMOVE ;

: !SINE		Sine>Wave DO-IT ;

: .DIFF%	Data &SCALE >S
		Data /ELEMENTS 2/ 
		1 ?DO #100
		      I Wave X@V  |X| 
		      I Data X@V  |X| 
		      - S */ 4 .R
		 LOOP -S ;

		Sine>Wave

		        (* End of information *)