Proton FFT örneğini kullanan oldumu

[Mucit23]

Arkadaşlar Protonun sample klasörünü karıştırırken Proton da uygulanmış FFT örneğini gördüm.

Aşağıdaki kodlar verilmiş.

FFT.bas

'
' 256 Point (128 point result) FFT Routine using Fixed Point Calculations
'
' Accepts an audio signal on analogue input pin AN0 (PORTA.0)
' Displays the spectrum on a Samsung KS0108 Graphic LCD
'
' For 18F PICmicro devices only using the Crownhill Proton+ BASIC compiler version 3.3.0.7 Onwards
'
	Device = 18F452
	Xtal = 8
'
' Setup the Graphic LCD for use with the Proton Development Board
'        
	LCD_DTPort = PORTD
	LCD_RSPin = PORTC.1
	LCD_ENPin = PORTE.0
	LCD_RWPin = PORTC.0
	LCD_CS1Pin = PORTE.1
	LCD_CS2Pin = PORTE.2
	LCD_Type = GRAPHIC
	Internal_Font = On
	Font_Addr = 0
	
	Warnings = OFF
	Reminders = OFF
	PLL_Req = On										' Operate at 32MHz from an 8MHz crystal
	Optimiser_Level = 3									' Maximum Optimisation
	Dead_Code_Remove = On								' Tighten the ASM code even further
	Unsigned_Dwords = On								' We're not using signed Dwords
 
'--------------------------------------------------------------------------------------------------

	Include "Variable Definitions.Inc"					' Load the Variable and Constant Definitions into the program
	
'--------------------------------------------------------------------------------------------------

	GoTo MAIN_PROGRAM									' Jump Over the FFT Subroutines

'--------------------------------------------------------------------------------------------------

	Include "FFT Subroutines.Inc"						' Load the FFT Subroutines into memory

'--------------------------------------------------------------------------------------------------
' Capture the audio input
' And fill the FFT_DATA Array with the scaled ADC Samples
'
AquireADC:
	NumberOfSamples = 0
	Repeat
		GO_DONE = 1
		While GO_DONE = 1 : Wend						' Wait for a 10-bit Left Justified (Fractional) ADC value
		AD_RESULT = AD_RESULT >> 2
		AD_RESULT.HighByte = AD_RESULT.HighByte - $20   ' Subtract $2000 (fixed point 1.0) So it's now from -1.0 to +1.0
		FFT_DATA[NumberOfSamples] = AD_RESULT			' Place the value into the array
		DelayUS 100										' Alter this value for lesser or greater capture speeds
		Inc NumberOfSamples
	Until OVERFLOW_FLAG = 1								' A loop of 256
	Return
'--------------------------------------------------------------------------------------------------
' The Main Program Starts here
'
' Notes		: In order to display the S:2:13 fixed point value, use a variation of the code:
'			:	Rlcf AccumAH,w : TempWord = AccumA : If STATUS.0 = 1 Then TempWord = -TempWord : Hrsout "-" 	' Check and display the Sign Bit
'			: 	TempFloat = TempWord / 8192 : Hrsout Dec3 TempFloat,13 											' Convert to a Float and Display
'			: Where AccumA holds the Fixed Point Value
'
MAIN_PROGRAM:
	Clear										' Clear all RAM before we start
'
' Setup the ADC for a 18F452 Device
'
	ADCON0 = %01000001							' Fosc/16, Channel 0, ADC Enabled
	ADCON1 = %01000010							' Left Justified, All PORTA Analogue, PORTE Digital
'
'
' Setup the ADC for a 18F4520 Device
'
'	ADCON0 = %00000001							' Channel 0 (AN0), ADC Enabled
'	ADCON1 = %00001000							' All PORTA Analogue, PORTE Digital
'	ADCON2 = %00000101							' Left Justified, Fosc/16
	
	Cls											' Clear the LCD
'
' Perform the FFT and display the spectrum on the Graphic LCD
'
	While 1 = 1	
		GoSub AquireADC												' Fill the FFT array with ADC readings	
		PerformFFT													' Perform the FFT on the data received
		CalculateLog												' Scale the result using a Logarithmic scale (Remark if not used)
		'CalculateLin												' Scale the result using a Linear scale (Remark if not used)
		'
		' Display the spectrum on the graphic LCD
		' Alter this routine as required.
		'
		NumberOfSamples = 127										' Ignore the first sample (DC Content)	
		RealIndex = 1
		Repeat		
			TempByte = RealBuffer[RealIndex]						' Read a value from the array
			If TempByte >= 63 Then TempByte = 63
			PreviousRealValue = ShadowRealBuffer[RealIndex]			' Read a value from the shadow array					
			If PreviousRealValue <> TempByte Then
				Line 0,NumberOfSamples,0,NumberOfSamples, PreviousRealValue
			EndIf
			If TempByte > 1 Then
				Line 1,NumberOfSamples,0,NumberOfSamples, TempByte
			EndIf
			ShadowRealBuffer[RealIndex] = TempByte					' Keep a record of the line for erasure
			Inc RealIndex
			Dec NumberOfSamples
		Until NumberOfSamples = 0									' Perform a loop of 128
	Wend
'--------------------------------------------------------------------------------------------------

Variable Definitions.inc

'
' Variable and Constant Definitions for Fixed Point FFT
'
' Arithmetic working variables
'
	Dim AccumA as Word SYSTEM  			' General Purpose Fixed Point 1:2:13 Format Accumulator 
	Dim AccumB as Word SYSTEM 			' General Purpose Fixed Point 1:2:13 Format Accumulator 
	Dim DivTemp as Byte SYSTEM 			' Divisor for FixedDivide
'
' Temporary variables for FFT and Powers subroutines
'
	Dim wr as Word SYSTEM				' General Purpose Fixed Point 1:2:13 Format
	Dim wpr as Word SYSTEM				' General Purpose Fixed Point 1:2:13 Format
	Dim wi as Word SYSTEM				' General Purpose Fixed Point 1:2:13 Format
	Dim wpi as Word SYSTEM				' General Purpose Fixed Point 1:2:13 Format
	Dim Theta as Word SYSTEM			' Used for FixedSin
	Dim mmax as Byte SYSTEM				' General Counter
	Dim m as Byte SYSTEM				' General Counter
	Dim j as Byte SYSTEM				' General Counter
	Dim istep as Byte SYSTEM			' General Counter
	Dim i as Byte SYSTEM				' Mainly used as an index
	Dim Int1 as Byte SYSTEM				' General Counter
	Dim Int2 as Byte SYSTEM				' General Counter
	Dim Int3 as Byte SYSTEM				' General Counter
	Dim Int4 as Byte SYSTEM				' General Counter
'
' Temporary variables
'
	Dim Temp1r as Word SYSTEM			' General Purpose Fixed Point 1:2:13 Format
	Dim Temp1i as Word SYSTEM			' General Purpose Fixed Point 1:2:13 Format
	Dim Temp2r as Word SYSTEM			' General Purpose Fixed Point 1:2:13 Format
	Dim Temp2i as Word SYSTEM			' General Purpose Fixed Point 1:2:13 Format
	Dim TempWord as Word SYSTEM			' General Purpose Fixed Point 1:2:13 Format	
	Dim MultTemp1 as Byte SYSTEM		' Used in FixedMult subroutine
	Dim MultTemp2 as Byte SYSTEM
	Dim MultTemp3 as Byte SYSTEM
	Dim x0 as Word SYSTEM 				' Used in FixedSin subroutine
	Dim x1 as Word SYSTEM
	Dim ScaleTemp as x0					' Used in ScaleLin and ScaleLog subroutines
	Dim LogResult as x1					' Result of ScaleLog
	Dim NumberOfSamples as Byte SYSTEM	' Holds the amount of ADC readings
	Dim TempByte as Byte SYSTEM
	
	Dim PreviousRealValue as Int1
	Dim PreviousYpos as Int2
	Dim PreviousYposShadow as Int3
	Dim RealIndex as Int4
'
' Variables to make the FFT faster
'	
	Dim wr_X_Temp2r as Word SYSTEM
	Dim wi_X_Temp2i as Word SYSTEM	
	Dim wr_X_Temp2i as Word SYSTEM
	Dim wi_X_Temp2r as Word SYSTEM	
	Dim DataSaveJ as wr_X_Temp2r
	Dim DataSaveJ1 as wi_X_Temp2i
	Dim DataSaveI as wr_X_Temp2i
	Dim DataSaveI1 as wi_X_Temp2r
	
	Dim DataSave1 as wr_X_Temp2r
	Dim DataSave2 as wi_X_Temp2i
	Dim DataSave3 as wr_X_Temp2i
	Dim DataSave4 as wi_X_Temp2r
'
' Complex Value FFT data buffer (256 * 2 bytes per value)
'
	Dim FFT_DATA[256] as Word
'	
' Final Real FFT Values. Uses the same RAM location as the FFT_DATA array
'
	Dim RealBuffer[128] as Byte at FFT_DATA
'
' Holds a record of the previous LCD display, for erasing the LCD line without flicker
'
	Dim ShadowRealBuffer[128] as Byte
'
' Aliases and Constants
'	
	Dim PROD as PRODL.Word
	Dim FSR0 as FSR0L.Word
	Dim TBLPTR as TBLPTRL.Word
	Dim AD_RESULT as ADRESL.Word
'
' For 18F452
'
	Symbol GO_DONE = ADCON0.2			' Change for different PICmicro's
'
' For 18F4520
'
'	Symbol GO_DONE = ADCON0.1			' Change for different PICmicro's
	
	Symbol Fixed_PI = $6488				' Fixed Point Format 1:2:13  PI (3.14159265459)
	Symbol MaxRealSize = 127			' Maximum scaled pixel value
	Symbol n = 256						' Number of FFT samples
	Symbol OVERFLOW_FLAG = STATUS.0		' Carry Flag

Misc Subroutines.inc

'
' Subroutines used in the Fixed Point FFT
'
'--------------------------------------------------------------------------------------------------
' Fixed Point Division of AccumA with DivTemp. (AccumA = AccumA / DivTemp)
' Input		: AccumA\AccumAH = Dividend  
' 			: DivTemp = Divisor
' Output	: AccumA\AccumAH
'
' Notes		: Uses Bit-0 of System Variable "PP0" as the Sign Holder
'			: Macro UDIV1608 uses PRODL\H as temporary variables
'			: Rolled out code for speed
'
	Symbol _SignBit = PP0.0
FixedDivide:
	_SignBit = 0				' Clear the Sign holder
	If AccumA.15 = 1 Then 		' Check the Sign of AccumA
		_SignBit = 1			' Set the Sign holder
		AccumA = -AccumA		' Negate AccumA
	Endif
	UDIV1608					' Do the division macro
	If _SignBit = 1 Then		' If AccumA was < 0
		AccumA = -AccumA		' Negate AccumA
	Endif
	Return
'--------------------------------------------------------------------------------------------------
' Fixed Point 1:2:13 Format Sin calculation
' AccumA = Sin(AccumA)
'
' Sin is calculated by interpolation of table values based upon the High Byte of AccumA
' value_msb and value_msb + 1 are the table values corresponding to High Byte of AccumA and the next value in the table
' Will rollover on 0
'
' AccumA = (((v_msb + 1) * lsb AccumA + (v_msb) * (1 - lsb AccumA))
'
SinLUT:
	CDATA as Word  	0,_
				   	255,   511,   766,   1021,  1274,  1527,  1777,  2026,_
					2273,  2518,  2760,  3000,  3237,  3470,  3700,  3927,_
					4150,  4368,  4583,  4793,  4998,  5198,  5393,  5583,_
					5768,  5947,  6120,  6287,  6448,  6603,  6751,  6893,_
					7028,  7156,  7277,  7391,  7498,  7597,  7689,  7774,_
					7850,  7920,  7981,  8035,  8081,  8119,  8149,  8171,_
					8185,  8191,  8189,  8179,  8162,  8136,  8102,  8060,_
					8011,  7953,  7888,  7815,  7735,  7647,  7551,  7448,_
					7338,  7221,  7097,  6965,  6827,  6682,  6531,  6373,_
					6210,  6040,  5864,  5682,  5495,  5303,  5105,  4902,_
					4695,  4483,  4266,  4046,  3821,  3593,  3361,  3126,_
					2888,  2647,  2404,  2158,  1910,  1660,  1408,  1156,_
					902,   647,   391,   135,   -120,  -375,  -631,  -886,_
					-1140, -1393, -1644, -1894, -2142, -2388, -2632, -2873,_
					-3111, -3347, -3579, -3807, -4032, -4253, -4469, -4682,_
					-4889, -5092, -5290, -5483, -5671, -5853, -6029, 6199,_
					6029,  5853,  5671,  5483,  5290,  5092,  4889,  4682,_
					4469,  4253,  4032,  3807,  3579,  3347,  3111,  2873,_
					2632,  2388,  2142,  1894,  1644,  1393,  1140,  886,_
					631,   375,   120,  -135,   -391,  -647,  -902,  -1156,_
					-1408, -1660, -1910, -2158, -2404, -2647, -2888, -3126,_
					-3361, -3593, -3821, -4046, -4266, -4483, -4695, -4902,_
					-5105, -5303, -5495, -5682, -5864, -6040, -6210, -6373,_
					-6531, -6682, -6827, -6965, -7097, -7221, -7338, -7448,_
					-7551, -7647, -7735, -7815, -7888, -7953, -8011, -8060,_
					-8102, -8136, -8162, -8179, -8189, -8191, -8185, -8171,_
					-8149, -8119, -8081, -8035, -7981, -7920, -7850, -7774,_
					-7689, -7597, -7498, -7391, -7277, -7156, -7028, -6893,_
					-6751, -6603, -6448, -6287, -6120, -5947, -5768, -5583,_
					-5393, -5198, -4998, -4793, -4583, -4368, -4150, -3927,_
					-3700, -3470, -3237, -3000, -2760, -2518, -2273, -2026,_
					-1777, -1527, -1274, -1021,  -766,  -511, -255
FixedSin:
'
' The ASM code below is a faster equivalent of:
'   x0 = Lread16 SinLUT[AccumA.HighByte]
'   x1 = Lread16 SinLUT[AccumA.HighByte + 1]
'
	TBLPTR = AccumA.HighByte
	TBLPTR = TBLPTR * 2
	Movlw (SinLUT & 255)
	Addwf TBLPTRL,F,0
	Movlw ((SinLUT >> 8) & 255)
	Addwfc TBLPTRH,F,0
	Tblrd*+
	x0.LowByte = TABLAT
	Tblrd*+
	x0.HighByte = TABLAT
	Tblrd*+
	x1.LowByte = TABLAT
	Tblrd*+
	x1.HighByte = TABLAT
						
	AccumA = AccumA * 256				' \ 
	AccumA = AccumA / 8					' / Convert to a weight factor
	TempWord = FixedMult [x1, AccumA]	' Multiply "next" table value by AccumA and save it					
	AccumA = $2000 - AccumA				' Calculate 1.0 - weight factor
	AccumA = FixedMult [AccumA, x0]		' And multiply it by "previous" value
	AccumA = AccumA + TempWord			' Add the 2 values
	Return
	
'**************************************************************************************************
'
' Power spectrum subroutines for the result of the real FFT.
' Result power is scaled (lin or log) and stored in the Real Buffer Byte Array "RealBuffer"
' One byte per frequency (1 to 126)
'
'**************************************************************************************************
' Calculate power spectrum for the complex value stored in AccumA and AccumB
' Scale it with a log scale (result from 0 to MaxRealSize) and store the result in LogResult.HighByte
'
' AccumA = (AccumA / 2) ^ 2 + (AccumB / 2) ^ 2 (Divide by 2 to avoid overflow)
'
CalculateLog MACRO
	Gosub PowerLog
	ENDM
#ifdef CalculateLog#Req
LogLUT:
	CDATA as Word -4652, -4599, -4548, -4496, -4445, -4394, -4344, -4294,_
				  -4245, -4196, -4147, -4098, -4050, -4002, -3955, -3908,_
				  -3861, -3815, -3769, -3723, -3677, -3632, -3587, -3543,_
				  -3498, -3454, -3410, -3367, -3324, -3281, -3238, -3196,_
				  -3154, -3112, -3071, -3029, -2988, -2947, -2907, -2867,_
				  -2826, -2787, -2747, -2708, -2669, -2630, -2591, -2553,_
				  -2514, -2476, -2438, -2401, -2363, -2326, -2289, -2253,_
				  -2216, -2180, -2143, -2107, -2072, -2036, -2001, -1965,_
				  -1930, -1895, -1861, -1826, -1792, -1758, -1724, -1690,_
				  -1656, -1623, -1590, -1556, -1523, -1491, -1458, -1425,_
				  -1393, -1361, -1329, -1297, -1265, -1234, -1202, -1171,_
				  -1140, -1109, -1078, -1047, -1017, -986,  -956,  -926,_
				  -896,  -866,  -836,  -806,  -777,  -748,  -718,  -689,_
				  -660,  -631,  -603,  -574,  -545,  -517,  -489,  -461,_
				  -433,  -405,  -377,  -349,  -322,  -294,  -267,  -240,_
				  -213,  -186,  -159,  -132,  -105,   -79,  -52,   -26
InternalPowerLog:
	Temp1r = Div2 [AccumA]					' Save AccumA
	AccumA = Div2 [AccumB]
	Temp2r = FixedMult [AccumA, AccumA]
	AccumA = FixedMult [Temp1r, Temp1r]
	AccumA = AccumA + Temp2r
	If AccumA < 32 Then
		LogResult.HighByte = 0
		Return
	Endif
	If AccumA > (32 * MaxRealSize) Then
		LogResult.HighByte = MaxRealSize
		Return
	Endif
	LogResult = 0							' Clear result register
											' AccumA is between 32 and (32 * MaxRealSize)		
	ScaleTemp = AccumA * 8					' Multiply ScaleTemp by 8 (ScaleTemp is between 256 and (256 * MaxRealSize))
'	
' Normalisation, until ScaleTemp < 256
'
	While ScaleTemp.HighByte <> 0
		LogResult = LogResult + 4652		' Add 256 * ln(2) * MaxRealSize / ln(MaxRealSize) ( = 4652) to current evaluation
		ScaleTemp = Div2[ScaleTemp]			' Signed Divide ScaleTemp by 2
	Wend									' And loop
'
' ScaleTemp is now between 128 and 255 so substract 128 and use it as index in LogLUT
'
' The ASM Code below is a faster version of: "ScaleTemp = Lread16 LogLUT[ScaleTemp - 128]"
'	
	TBLPTR = ScaleTemp - 128
	TBLPTR = TBLPTR * 2
	Movlw (LogLUT & 255)
	Addwf TBLPTRL,F,0
	Movlw ((LogLUT >> 8) & 255)
	Addwfc TBLPTRH,F,0
	Tblrd*+
	ScaleTemp.LowByte = TABLAT
	Tblrd*
	ScaleTemp.HighByte = TABLAT
	LogResult = LogResult + ScaleTemp 		' Add this value to result
	Return
'--------------------------------------------------------------------------------------------------
' Calculate power spectrum for the complex spectrum stored in Word Array "FFT_DATA"
' Scale it with a log scale (result from 1 to MaxRealSize) and store the result in Byte Array "RealBuffer"
'
' Notes 	: Element 1 (DC) is discarded
'
PowerLog:
	i = 1
	Repeat
		Int2 = i * 2
		AccumA = FFT_DATA[Int2]
		AccumB = FFT_DATA[Int2 + 1]
		Gosub InternalPowerLog
		RealBuffer[i] = LogResult.HighByte
		Inc i
	Until i.7 = 1							' A Loop of 128
	Return	
#endif
'--------------------------------------------------------------------------------------------------
' Calculate power spectrum for the complex value stored in AccumA and AccumB
' Scale it with a linear scale (result from 0 to MaxRealSize) and store the result in ScaleTemp.LowByte
'
' AccumA = (AccumA / 2) ^ 2 + (AccumB / 2) ^ 2 (divide by 2 to exclude overflow)
'
CalculateLin MACRO
	Gosub PowerLin
	ENDM
#ifdef CalculateLin#Req
InternalPowerLinear:
	Temp1r = AccumA							' Save AccumA
	AccumA = Div2 [Temp1r]
	Temp2r = FixedMult [AccumA, AccumA]
	AccumA = Div2 [Temp1r]
	AccumA = FixedMult [AccumA, AccumA]
	AccumA = AccumA + Temp2r
	ScaleTemp = AccumA / 8
	If ScaleTemp > 63 Then ScaleTemp = 63
	Return
'--------------------------------------------------------------------------------------------------
' Calculate power spectrum for the complex spectrum stored in Word Array "FFT_DATA"
' Scale it with a linear scale (result from 1 to MaxRealSize) and store the result in Byte Array "RealBuffer"
'
' Notes 	: Element 1 (DC) is discarded
'
PowerLin:
	i = 1
	Repeat
		Int2 = i * 2
		AccumA = FFT_DATA[Int2]
		AccumB = FFT_DATA[Int2 + 1]
		Gosub InternalPowerLinear
		RealBuffer[i] = ScaleTemp.LowByte
		Inc i
	Until i.7 = 1							' A Loop of 128
	Return
#endif
'--------------------------------------------------------------------------------------------------

FFT Subroutines.inc

'
' 256 point FFT routines, using S:2:13 Fixed Point Format for numeric computations
' All data is stored in the the Word Array "FFT_DATA".
'
	Include "FFT Macros.Inc"					' Load the macros into the program
	Include "Misc Subroutines.Inc"				' Load the General Purpose Subroutines into memory
	
'--------------------------------------------------------------------------------------------------
' Do a Complex FFT, radix-1, decimation in place, 128 points
' Input 	: 128 complex values (4 bytes), stored in FFT_DATA[0..255]
'          		For each data element: 
'					Byte 0 = LSB of Real part
'                   Byte 1 = MSB of Real part
'                  	Byte 2 = LSB of Complex part
'                   Byte 3 = MSB of Complex part
'
' Output	: Replace FFT_DATA[0..255] by its FFT (128 frequencies, Complex Values)
' Notes		: Algorithm adapted from the book "Numerical recipes in C"
'
PerformFFT Macro
	GoSub FFTSub
	Endm
FFTSub:
'
' Make sure we're pointing to code memory for table reads
'
	EECON1 = 0
	bsf EECON1,PP_EEPGD
'
' Bit Reversal
'
	j = 0
	i = 0											
	Repeat
		If j > i Then		
			AccumA = FFT_DATA[i]			' Save the contents of FFT_DATA[i]
			FFT_DATA[i] = FFT_DATA[j]
			FFT_DATA[j] = AccumA		
			AccumA = FFT_DATA[i + 1]		' Save the contents of FFT_DATA[i + 1]
			FFT_DATA[i + 1] = FFT_DATA[j + 1]
			FFT_DATA[j + 1] = AccumA
		EndIf
		m = n / 2
		While 1 = 1
			If m < 2 Then Break
			If j < m Then Break
			j = j - m
			m = m / 2
		Wend
		j = j + m	
		i = i + 2
	Until OVERFLOW_FLAG = 1 				' A Loop of 0 to 255 i.e. n - 1
'
' The Actual FFT routine
'
  	mmax = 2
	While 1 = 1
		If mmax = 0 Then 
			wpr = $C008
			wpi = $FFFF
			Break
		EndIf
		istep = mmax * 2
		
		' More Efficient equivalent of: Theta = 2 * pi / mmax
		AccumA = Fixed_PI					' 3.14159265459
		DivTemp = mmax / 2
		GoSub FixedDivide
		Theta = AccumA				

		' More Efficient equivalent of: wpr = -2 * Sin(Theta / 2) ^ 2
		AccumA = Div2 [AccumA]
		GoSub FixedSin				
		AccumB = FixedMult [AccumA, AccumA]						
		wpr = FixedMult [AccumB, $C000]		' AccumB * -2.0

		' More Efficient equivalent of: wpi = Sin(Theta)
		AccumA = Theta				
		GoSub FixedSin				
		wpi = AccumA		
		wr = $2000 ' wr = 1.0
		wi = $0000 						' wi = 0.0
		m = 0
		Repeat
			i = m
Indexloop:
			j = i + mmax
			'
			' Pre-Capture values from the FFT_DATA array to speed up the routine
			'
			DataSaveJ = FFT_DATA[j]
			DataSaveJ1 = FFT_DATA[j + 1]
			DataSaveI = FFT_DATA[i]
			DataSaveI1 = FFT_DATA[i + 1]
				
			' More Efficient equivalent of: Temp1r = wr * FFT_DATA[j] - wi * FFT_DATA[j + 1]
			If wi <> 0 Then					' Save an unnecessary multiply
				AccumA = FixedMult [wr, DataSaveJ]
				Temp1r = FixedMult [wi, DataSaveJ1]
				Temp1r = AccumA - Temp1r
			Else
				Temp1r = FixedMult [wr, DataSaveJ]
			EndIf
				
			' More Efficient equivalent of: Temp1i = wr * FFT_DATA[j + 1] + wi * FFT_DATA[j]
			If wi <> 0 Then					' Save an unnecessary multiply
				Temp1i = FixedMult [wi, DataSaveJ]
				AccumA = FixedMult [wr, DataSaveJ1]
				Temp1i = Temp1i + AccumA
			Else
				Temp1i = FixedMult [wr, DataSaveJ1]
			EndIf
			'
			' Butterfly update with Scaling
			'
			' More Efficient equivalent of: FFT_DATA[j] = (FFT_DATA[i] - Temp1r) / 2
			AccumA = DataSaveI - Temp1r
			AccumA = Div2 [AccumA]			
			FFT_DATA[j] = AccumA			

			' More Efficient equivalent of: FFT_DATA[j + 1] = (FFT_DATA[i + 1] - Temp1i) / 2
			AccumA = DataSaveI1 - Temp1i
			AccumA = Div2 [AccumA]
			FFT_DATA[j + 1] = AccumA

			' More Efficient equivalent of: FFT_DATA[i] = (FFT_DATA[i] + Temp1r) / 2
			AccumA = DataSaveI + Temp1r
			AccumA = Div2 [AccumA]
			FFT_DATA[i] = AccumA

			' More Efficient equivalent of: FFT_DATA[i + 1] = (FFT_DATA[i + 1] + Temp1i) / 2
			AccumA = DataSaveI1 + Temp1i
			AccumA = Div2 [AccumA]		
			FFT_DATA[i + 1] = AccumA
			
			If istep <> 0 Then
				i = i + istep
				If OVERFLOW_FLAG = 0 Then GoTo Indexloop
			EndIf			
			AccumB = wr						' Save the contents of wr
			
			' More Efficient equivalent of: wr = wr * wpr - wi * wpi + wr
			Temp1r = FixedMult [wi, wpi]
			AccumA = FixedMult [wr, wpr]
			AccumA = AccumA - Temp1r
			wr = wr + AccumA

 			' More Efficient equivalent of: wi = wi * wpr + AccumB * wpi + wi
			Temp1i = FixedMult [AccumB, wpi]
			AccumA = FixedMult [wi, wpr]
			AccumA = AccumA + Temp1i
			wi = wi + AccumA
			m = m + 2
		Until m >= mmax - 1
		mmax = istep
	Wend
'
' Find the Values from the above FFT result
'
	Theta = $00C9 					    ' Theta = 2 * pi / n
	wpr = $FFFE			 				' wpr = -2 * Sin(Theta / 2) ^ 2
	wpi = $00C8			 				' wpi = Sin(2 * pi / n)			
	wr = $1FFE			  				' wr = wpr + 1.0
	wi = wpi
	i = 2
	Repeat
		Int1 = i + i
		Int1 = Int1 - 2
		Int2 = Int1 + 1
		Int3 = -Int2
		Int3 = Int3 + 1
		Int4 = Int3 + 1
		'
		' Pre-Capture values from the FFT_DATA array to speed up the routine
		'
		DataSave1 = FFT_DATA[Int1]
		DataSave2 = FFT_DATA[Int2]
		DataSave3 = FFT_DATA[Int3]
		DataSave4 = FFT_DATA[Int4]
		
		' More Efficient equivalent of: Temp1r = (FFT_DATA[Int1] + FFT_DATA[Int3]) / 2
		Temp1r = DataSave1 + DataSave3
		Temp1r = Div2 [Temp1r]					
		
		' More Efficient equivalent of: Temp1i = (FFT_DATA[Int2] - FFT_DATA[Int4]) / 2
		Temp1i = DataSave2 - DataSave4
		Temp1i = Div2 [Temp1i]						

		' More Efficient equivalent of: Temp2r = (FFT_DATA[Int2] + FFT_DATA[Int4]) / 2		
		Temp2r = DataSave2 + DataSave4
		Temp2r = Div2 [Temp2r]			
		
		' More Efficient equivalent of: Temp2i = -(FFT_DATA[Int1] - FFT_DATA[Int3]) / 2			
		Temp2i = DataSave3 - DataSave1
		Temp2i = Div2 [Temp2i]					
		'
		' Pre-Multiply key variables to streamline the FFT
		'
		wi_X_Temp2i = FixedMult [wi, Temp2i]
		wr_X_Temp2r = FixedMult [wr, Temp2r]	
		wr_X_Temp2i = FixedMult [wr, Temp2i]
		wi_X_Temp2r = FixedMult [wi, Temp2r]
	
		' More Efficient equivalent of: FFT_DATA[Int1] = Temp1r + wr * Temp2r - wi * Temp2i
		AccumA = wr_X_Temp2r - wi_X_Temp2i
		AccumA = AccumA + Temp1r
		FFT_DATA[Int1] = AccumA

		' More Efficient equivalent of: FFT_DATA[Int2] = Temp1i + wr * Temp2i + wi * Temp2r
		AccumA = wr_X_Temp2i + wi_X_Temp2r
		AccumA = AccumA + Temp1i
		FFT_DATA[Int2] = AccumA

		' More Efficient equivalent of: FFT_DATA[Int3] = Temp1r - wr * Temp2r + wi * Temp2i	
		AccumA = wi_X_Temp2i - wr_X_Temp2r
		AccumA = AccumA + Temp1r
		FFT_DATA[Int3] = AccumA

		' More Efficient equivalent of: FFT_DATA[Int4] = -Temp1i + wr * Temp2i + wi * Temp2r
		AccumA = wi_X_Temp2r + wr_X_Temp2i
		AccumA = AccumA - Temp1i
		FFT_DATA[Int4] = AccumA		
		AccumB = wr						' Save the value of wr

		' More Efficient equivalent of: wr = wr * wpr - wi * wpi + wr
		Temp1r = FixedMult [wi, wpi]
		AccumA = FixedMult [wr, wpr]
		AccumA = AccumA - Temp1r
		wr = wr + AccumA

		' More Efficient equivalent of: wi = wi * wpr + AccumB * wpi + wi
		Temp1r = FixedMult [wi, wpr]
		AccumA = FixedMult [AccumB, wpi]
		AccumA = AccumA + Temp1r
		wi = wi + AccumA
		Inc i
	Until i > (n / 4)
'
' Special case of the first two data elements
'
	AccumA = FFT_DATA#0						' Save the value of FFT_DATA[0]
	FFT_DATA#0 = FFT_DATA#0 + FFT_DATA#1
	FFT_DATA#1 = AccumA - FFT_DATA#1
	Return

FFT Macros.inc

'
' Macros used in the Fixed Point FFT
'
' Note that none of the macros perform any RAM bank switching so ensure that all variables used are in Access RAM
'
' Bring some of the compiler's System Variables into the BASIC program for use with the Fixed Point Multiplier
'
    Dim PP0 as Byte SYSTEM
    Dim PP0H as Byte SYSTEM
    Dim PP1 as Byte SYSTEM
    Dim PP1H as Byte SYSTEM
'--------------------------------------------------------------------------------------------------
' Signed Divide 16-bit Parameter By 2 (Var = Parameter1 / 2)
'
Div2 MACRO Parameter1
#if(DIV2_RETURN == 0)                ; If no return variable given, do the following code
ASM-
    Bcf STATUS,C
    Rrcf Parameter1 + 1,f            ; Divide high byte
    Rrcf Parameter1,f                ; and low byte
    Btfsc Parameter1 + 1,6            ; Test old sign bit
    Bsf Parameter1 + 1,7            ; If set, alter the msb
ENDASM-
#else
    #if(RETURN_VAR == Parameter1)     ; Is the return variable the same as the parameter ?
ASM-
    Bcf STATUS,C
    Rrcf Parameter1 + 1,f            ; Divide high byte
    Rrcf Parameter1,f                ; and low byte
    Btfsc Parameter1 + 1,6            ; Test old sign bit
    Bsf Parameter1 + 1,7            ; If set, alter the msb    
ENDASM
    #else
ASM-
    Movff Parameter1,RETURN_VAR
    Movff Parameter1 + 1,RETURN_VAR + 1
    Bcf STATUS,C
    Rrcf RETURN_VAR + 1,f            ; Divide high byte
    Rrcf RETURN_VAR,f                ; and low byte
    Btfsc RETURN_VAR + 1,6            ; Test old sign bit
    Bsf RETURN_VAR + 1,7            ; If set, alter the msb        
ENDASM
    #endif
#endif
    ENDM
'--------------------------------------------------------------------------------------------------
' Signed 16-bit Multiply
' Notes        : Overflow not detected if Fixed Point result > 4.0
'
FixedMult MACRO Multiplier, Multiplicand
#if(Multiplier == Multiplicand)                    ; Multiplier and Multiplicand are the same variable?
    #if((PRM_1 == NUM8) || (PRM_1 == NUM16) || (PRM_1 == NUM32))
        Num_Word Multiplier,PP0
        Num_Word Multiplier,PP1
    #endif
    #if(PRM_1 == BYTE)
ASM-
    Movf Multiplier,W
    Movwf PP0
    Movwf PP1
ENDASM
    #endif
    #if((PRM_1 == WORD) || (PRM_1 == DWORD))
ASM-
    Movf Multiplier,W
    Movwf PP0
    Movwf PP1
    Movf Multiplier + 1,W
    Movwf PP0H
    Movwf PP1H
ENDASM
    #endif
#else
    #if((PRM_1 == NUM8) || (PRM_1 == NUM16) || (PRM_1 == NUM32))
        Num_Word Multiplier,PP0
    #endif
    #if(PRM_1 == BYTE)
        Byte_Word Multiplier,PP0
    #endif
    #if((PRM_1 == WORD) || (PRM_1 == DWORD))
        Word_Word Multiplier,PP0
    #endif
    #if((PRM_2 == NUM8) || (PRM_2 == NUM16) || (PRM_2 == NUM32))
        Num_Word Multiplicand,PP1
    #endif
    #if(PRM_2 == BYTE)
        Byte_Word Multiplicand,PP1
    #endif
    #if((PRM_2 == WORD) || (PRM_2 == DWORD))
        Word_Word Multiplicand,PP1
    #endif
#endif
    Gosub FixedMultSub
ASM-
    Movff MultTemp3, RETURN_VAR + 1    ; Store the result..
    Movff MultTemp2, RETURN_VAR        ; in return variable
ENDASM
    ENDM
'---------------------------------------------------------------------------
' Fixed Point Multiply Subroutine used by the FixedMult Macro
'
' Input        : PP0 = Multiplier
'              PP1 = Multiplicand
' Output    : MultTemp2 and MultTemp3 hold the 16-bit result
'
#ifdef FixedMult#req
FixedMultSub:
    Movf PP0,w
    Mulwf PP1
    Movff PRODH,MultTemp1
    Movf PP0H,w
    Mulwf PP1H
    Movff PRODL,MultTemp2
    Movff PRODH,MultTemp3
    Movf PP0,w
    Mulwf PP1H
    Movf PRODL,w
    Addwf MultTemp1,f
    Movf PRODH,w
    Addwfc MultTemp2,f
    Movlw 0
    Addwfc MultTemp3,f
    Movf PP0H,w
    Mulwf PP1
    Movf PRODL,w
    Addwf MultTemp1,f
    Movf PRODH,w
    Addwfc MultTemp2,f
    Movlw 0
    Addwfc MultTemp3,f
    Btfss PP1H,7      		; Is PP1 < 0?
    Bra $ + 10
    Movf PP0,w            	; Yes. So, adjust value
    Subwf MultTemp2,f
    Movf PP0H,w
    Subwfb MultTemp3,f
    Btfss PP0H,7            ; Is PP0 < 0?
    Bra $ + 10
    Movf PP1,w           	; Yes. So, adjust value
    Subwf MultTemp2,f
    Movf PP1H,w
    Subwfb MultTemp3,f
    Rlcf MultTemp1,f                
    Rlcf MultTemp2,f
    Rlcf MultTemp3,f
    Rlcf MultTemp1,f                
    Rlcf MultTemp2,f
    Rlcf MultTemp3,f
    Rlcf MultTemp1,f                
    Rlcf MultTemp2,f
    Rlcf MultTemp3,f
    Return
#endif    
'---------------------------------------------------------------------------
' Signed Fixed Point divide routine
'
' Input        : Arguments in AccumA and DivTemp
' Output    : Quotient AccumA
'
' Notes        : Rolled out for speed at the cost of code size
'            : Uses PRODL\H as temporary variables
'
UDIV1608  MACRO
ASM-
    Clrf PRODL
NOLIST
 Variable DivEntities = 0
 while (DivEntities < 8)
LIST
    Rlcf AccumAH,w
    Rlcf PRODL,f
    Movf DivTemp,w
    Subwf PRODL,f
    Bc $ + 6
    Addwf PRODL,f
    Bcf STATUS,C
    Rlcf AccumAH,f
NOLIST
 Variable DivEntities = DivEntities + 1
 endw
LIST
    Clrf PRODH
NOLIST
 Variable DivEntities = 8
 while DivEntities < 16
LIST
    Rlcf AccumA,w
    Rlcf PRODL,f
    Rlcf PRODH,f
    Movf DivTemp,w
    Subwf PRODL,f
    Movlw 0
    Subwfb PRODH,f
    Bc $ + 12
    Movf DivTemp,w
    Addwf PRODL,f
    Movlw 0
    Addwfc PRODH,f
    Bcf STATUS,C
    Rlcf AccumA,f
NOLIST
 Variable DivEntities = DivEntities + 1
 endw
LIST
ENDASM
    ENDM
'--------------------------------------------------------------------------------------------------
   

Bi dünya kod var burada. Anlaması zor. Daha önce deneyen oldumu bunları. Öğrenmek istiyorum bu işin mantığını

[t2]

Mantığını  ne yapacaksın,  çalışmıyor mu bunlar?
veya beğenmediğin bir yeri var da onu mu değiştirmek istiyorsun?

[Mucit23]

Çalışıp çalışmadığını bilmiyorum. Ama Proton derliyor. İsistede GLCD ile çalıştırdım fakat ekrana ne olduğu anlamsız şeyler çıkıyor.
Vumetre yapmak istiyorum.

[Mucit23]

Aslında daha önce bu konuyu hiç incelemedim. Biraz araştırma yaptım. Goertzel Algoritması daha sade uygulaması kolay görünüyor. FFT ile Goertzel algooritması arasında ne gibi farklar var?

[yamak]

Hocam vumetre için fft gerekmez ki.Auido spektrum analyzer ın için fft ye gerek var.

[yamak]

Fft ile goertzel arasındaki fark fft belli bir aralıktaki sinyal bileşenlerini verir, goertzel ise belirlenen sinyal bileşenlerini gösterir.

[Mucit23]

Vumetre şakaydı Ama basitten bir ( En azından 8 Band) spektrum analyser gerçekten yapmak istiyorum.
Sizin dediğinize göre Gorertzel agooritmasında hangi frekansları istiyorsak bize  Sinyal içerisindeki istediğimiz frekansın ne kadar  baskın olduğunu gösteriyor.
İnternette birkaç yazı gördüm goertzel ile ilgili. Az çok mantığını anladım fakat yapılan bu işlemler 16 Bit sınırlar içerisinde tamamlanabiliyormu, Yoksa 32 bit değişkenlermi gerekiyor.

[Hattusa]

s.a.
mucid usta verdiğin kodun şeması ve derlemesini paylaşabilirmisiniz? bende merak ediyorum açıkcası bu FFT yi ve goertzel algoritmasını

[GreeN]

fft ve Goertzel arasındaki farklardan biride gerekli olan örnek sayısı , fft 2^n kadar örnek işlemesi gerekirken Goertzel'de örnek sayısı önemli olmuyor. Belli bir frekans takip edilecekse Goertzel fft'ye göre hızlı.

[Mucit23]

Şeması yok malesef.
Dün öyle ayaküstü kodu denemek için koda bakarak şemayı çizdim. Analoğ Giriş An0 olacak. GLCD pinleride yazılımda belirtilmiş.

[Mucit23]

fft ve Goertzel arasındaki farklardan biride gerekli olan örnek sayısı , fft 2^n kadar örnek işlemesi gerekirken Goertzel'de örnek sayısı önemli olmuyor. Belli bir frekans takip edilecekse Goertzel fft'ye göre hızlı.

Hocam Mesela 20 Band spektrum analyser için Bellirli bir frekanslar vardır. Biz goertzel ile sadece bu belirlediğimiz frekansların Baskınlığını ölçüyoruz.
Diğer mesajımdada sormuştum.Yapılan işlemler 16 bit sınırlarını geçiyormu?

[GreeN]

uygulamayı pic vs mcu ile hiç denemedim , ama c builder ile yaptım  , ayrıca fxdev uygulamayı pic ile yapmıştı 18f xx kullanmıştı sanırım ,

    k=(freq/Ornek_Freq);
    coeff = 2*cos(2*M_PI*k);
    for (i=0; i<N; i++) {
        s = samples[i] + coeff * Ornek_Onceki - Ornek_Onceki2;
        Ornek_Onceki2 = Ornek_Onceki;
        Ornek_Onceki = s;

    }
      reel=  Ornek_Onceki-Ornek_Onceki2*cos(2*M_PI*k);
     imag = Ornek_Onceki2*sin(2*M_PI*k);
    power = sqrt(reel*reel+imag*imag)*2/N;

programın kodu bundan ibaret. float işlemler için dspic kullanmak (yada dsp türü mcu lar) sanırım daha uygun olur.

[Mucit23]

Bu işlemleri yapmadan önce adc den N kadar örnek almak gerekiyor değilmi. Yani Tüm örneklerimizi alıp N kadar değişkene yerleştirdikten sonra Goertzel işlemlerine başlıyoruz. İşlem arasında Yeniden ADC okuması yapılmıyor.

[GreeN]

Bu işlemleri yapmadan önce adc den N kadar örnek almak gerekiyor değilmi. Yani Tüm örneklerimizi alıp N kadar değişkene yerleştirdikten sonra Goertzel işlemlerine başlıyoruz. İşlem arasında Yeniden ADC okuması yapılmıyor.

Tüm örnekleri değişkende tutulması lazım , algoritmada tüm örnekler işleniyor yani misal adc 10 örnek alacak arkasında algoritma çalışacak ,  hızlı bir algoritma hiç ölçmedim ama hızı konusunda meth ediliyor.

dspic30f4012 ile bir çalışmam olmuştu , 16 adc okuyup kesme üretiyordu , adc çevrimi sırasında işlemleri gayet rahat bitirebilisiniz. adc çevrimi bitince diğer örnekleri alıp işleme devam edebilirsiniz. Hele DMA'lı bir PIC kullanırsanız  işler daha kolay.

[Mucit23]

Anladım şimdi demek istediğinizi.
Yukarıda En sonda Power isimli bir değişken var. Bu değişken Goertzel in sonucu olsa gerek. Biz hangi frekansı istiyorsak o frekansa göre verdiğiniz işlemleri tekrarlamamız gerekiyor. Ve bu işlemlerin sonucunda Power değişkeni sinyalimizin baskınlık değerini veriyor. Peki Bu değer Hangi aralıkta değişiyor. Daha doğrusu 0 ile XX arasında değişen bir değermi?