'From Cuis 2.5 of 9 August 2010 [latest update: #534] on 10 August 2010 at 11:20:42 am'!

!FFT methodsFor: 'initialization' stamp: 'jmv 8/10/2010 10:51'!
basicRealData: real
	"Avoids conversion to FloatArray.
	This disables the fast plugin, but allows for 64 bit Double precision:"
	self
		basicRealData: real
		imagData: (real collect: [:i | 0.0])  "imaginary component all zero"! !

!FFT methodsFor: 'initialization' stamp: 'jmv 8/10/2010 10:59'!
basicRealData: real imagData: imag
	"Avoids conversion to FloatArray.
	If data is not in FloatArrays,  the fast plugin is disabled, but this allows for 64 bit Double precision:
	In any case, make the sinTable match the data, to allow either the plugin or the full precision"
	(real is: #FloatArray)
		ifTrue: [
			(sinTable is: #FloatArray)
				ifFalse: [ sinTable _ sinTable asFloatArray ]]
		ifFalse: [
			(sinTable is: #FloatArray)
				ifTrue: [ sinTable _ (0 to: n/4) collect: [:i | (i asFloat / (n//4) * Float pi / 2.0) sin ]]].
	realData _ real.
	imagData _ imag! !

!FFT methodsFor: 'initialization' stamp: 'jmv 8/10/2010 10:43'!
nu: order
	"Initialize variables and tables for transforming 2^nu points"
	|  j perms k |
	nu _ order.
	n _ 2 bitShift: nu-1.

	"Initialize permutation table (bit-reversed indices)"
	j_0.
	perms _ WriteStream on: (Array new: n).
	0 to: n-2 do:
		[:i |
		i < j ifTrue: [perms nextPut: i+1; nextPut: j+1].
		k _ n // 2.
		[k <= j] whileTrue: [j _ j-k.  k _ k//2].
		j _ j + k].
	permTable _ perms contents asWordArray.

	"Initialize sin table 0..pi/2 in n/4 steps."
	sinTable _ (0 to: n/4) collect: [:i | (i asFloat / (n//4) * Float pi / 2.0) sin].
	
	sinTable _ sinTable asFloatArray.
	realData _ FloatArray new: n.
	imagData _ FloatArray new: n.

	self initializeHammingWindow: 0.54.  "0.54 for Hamming, 0.5 for hanning"! !

!FFT methodsFor: 'initialization' stamp: 'jmv 8/10/2010 10:49'!
realData: real
	"By defaults use FloatArrays to allow for the fast primitive.
	This means precision is limeted to 32bit Float.
	For 64 bit Double precision, use #basicRealData:"
	self basicRealData: real asFloatArray! !

!FFT methodsFor: 'initialization' stamp: 'jmv 8/10/2010 10:49'!
realData: real imagData: imag
	"By defaults use FloatArrays to allow for the fast primitive.
	This means precision is limeted to 32bit Float.
	For 64 bit Double precision, use #basicRealData:imagData:"
	self basicRealData: real asFloatArray imagData: imag asFloatArray! !

!FFT methodsFor: 'transforming' stamp: 'jmv 8/10/2010 11:03'!
transformForward: forward
	| lev lev1 ip theta realU imagU realT imagT i |

	"Use the primitive if available and if data is in the correct format (FloatArrays and a WordArray for permTable)"
	<primitive: 'primitiveFFTTransformData' module: 'FFTPlugin'>

	self permuteData.
	1 to: nu do:
		[:level |
		lev _ 1 bitShift: level.
		lev1 _ lev // 2.
		1 to: lev1 do:
			[:j |
			theta _ j-1 * (n // lev).   "pi * (j-1) / lev1 mapped onto 0..n/2"
			theta < (n//4)  "Compute U, the complex multiplier for each level"
				ifTrue:
					[realU _ sinTable at: sinTable size - theta.
					imagU _ sinTable at: theta + 1]
				ifFalse:
					[realU _ (sinTable at: theta - (n//4) + 1) negated.
					imagU _ sinTable at: (n//2) - theta + 1].
			forward ifFalse: [imagU _ imagU negated].
"
			Here is the inner loop...
			j to: n by: lev do:
				[:i |   hand-transformed to whileTrue...
"
			i _ j.
			[i <= n] whileTrue:
				[ip _ i + lev1.
				realT _ ((realData at: ip) * realU) - ((imagData at: ip) * imagU).
				imagT _ ((realData at: ip) * imagU) + ((imagData at: ip) * realU).
				realData at: ip put: (realData at: i) - realT.
				imagData at: ip put: (imagData at: i) - imagT.
				realData at: i put: (realData at: i) + realT.
				imagData at: i put: (imagData at: i) + imagT.
				i _ i + lev]]].
	forward ifFalse: [self scaleData]  "Reverse transform must scale to be an inverse"! !

!FFT methodsFor: 'testing' stamp: 'jmv 8/10/2010 11:11'!
testFullPrecision
	 "
	Display restoreAfter: [(FFT new nu: 15) testFullPrecision].
	  --  Test on an array of 32768 samples"
	"Initialize to pure (co)Sine Wave, plot, transform, plot, invert and plot again.
	Allow for full 64bit Double precision, do not use the fast plugin."
	self basicRealData: ((1 to: n) collect: [:i | (Float pi * (i-1) / (n/8)) cos]).
	self plot: realData in: (100@20 extent: 256@60).
	Transcript cr; print: (Time millisecondsToRun: [ self transformForward: true ]); endEntry.
	self plot: realData in: (100@100 extent: 256@60).
	self plot: imagData in: (100@180 extent: 256@60).
	Transcript cr; print: (Time millisecondsToRun:[self transformForward: false]); endEntry.
	self plot: realData in: (100@260 extent: 256@60)! !

!FFT methodsFor: 'plugin-testing' stamp: 'jmv 8/10/2010 11:12'!
pluginTest  
	"
	Display restoreAfter: [(FFT new nu: 15) pluginTest].
	"
	"Test on an array of 32768 samples"
	"Initialize to pure (co)Sine Wave, plot, transform, plot, invert and plot again"
	self realData: ((1 to: n) collect: [:i | (Float pi * (i-1) / (n/8)) cos]).
	self plot: realData in: (100@20 extent: 256@60).
	Transcript cr; print: (Time millisecondsToRun:[self pluginTransformData: true]); endEntry.
	self plot: realData in: (100@100 extent: 256@60).
	self plot: imagData in: (100@180 extent: 256@60).
	Transcript cr; print: (Time millisecondsToRun:[self pluginTransformData: false]); endEntry.
	self plot: realData in: (100@260 extent: 256@60)! !

!FFT methodsFor: 'bulk processing' stamp: 'jmv 8/10/2010 10:42'!
setSize: anIntegerPowerOfTwo
	"Initialize variables and tables for performing an FFT on the given number of samples. The number of samples must be an integral power of two (e.g. 1024). Prepare data for use with the fast primitive."

	self nu: (anIntegerPowerOfTwo log: 2) asInteger.
	n = anIntegerPowerOfTwo ifFalse: [self error: 'size must be a power of two']
! !

!FFT methodsFor: 'bulk processing' stamp: 'jmv 8/10/2010 11:02'!
transformDataFrom: anIndexableCollection startingAt: index
	"Forward transform a block of real data taken from from the given indexable collection starting at the given index. Answer a block of values representing the normalized magnitudes of the frequency components."

	| j real imag out |
	j _ 0.
	index to: index + n - 1 do: [:i |
		realData at: (j _ j + 1) put: (anIndexableCollection at: i)].
	realData *= window.
	imagData _ FloatArray new: n.
	self transformForward: true.

	"compute the magnitudes of the complex results"
	"note: the results are in bottom half; the upper half is just its mirror image"
	real _ realData copyFrom: 1 to: (n / 2).
	imag _ imagData copyFrom: 1 to: (n / 2).
	out _ (real * real) + (imag * imag).
	1 to: out size do: [:i | out at: i put: (out at: i) sqrt].
	^ out
! !


!FloatArray methodsFor: 'testing' stamp: 'jmv 8/10/2010 10:55'!
is: aSymbol
	^aSymbol = #FloatArray or: [ super is: aSymbol ]! !


!FloatArray reorganize!
('accessing' at: at:put: defaultElement length squaredLength)
('arithmetic' * *= + += - -= / /= adaptToNumber:andSend: dot: negated)
('comparing' = hash)
('primitives-plugin' primAddArray: primAddScalar: primDivArray: primDivScalar: primMulArray: primMulScalar: primSubArray: primSubScalar: sum)
('converting' asFloatArray)
('private' replaceFrom:to:with:startingAt:)
('user interface' inspect inspectWithLabel:)
('testing' is:)
!

FFT removeSelector: #pluginPrepareData!
ArrayLiteralTest removeSelector: #array!
