finding the frequencies in an analog signal

Last post
23 posts / 0 new
Author
Message
#1
  • 1
  • 2
  • 3
  • 4
  • 5
Total votes: 0

I am using an at90usb1287 and working in asm. Currently, my program collects 256 data points (1 byte each) and stores them to SRAM. I would like to add a block to my program that takes these 256 data points and finds all of the frequencies in the signal. FFT comes to mind, but I am open to other ideas. I have looked at some FFT code in asm and cannot seem to figure it out. I am very new to all of this stuff, so maybe someone with more experience could help me figure out how to implement some FFT code. Also, at this point, I am not particularly concerned with accuracy or optimization...I just want something that works.

Thanks,
Sam

  • 1
  • 2
  • 3
  • 4
  • 5
Total votes: 0

 

  • 1
  • 2
  • 3
  • 4
  • 5
Total votes: 0

The FFT is the usual technique, I don't think you will find anything else unless you know the frequencies in advance. Experiment with it on the PC before you try getting it to run on an AVR.

Leon

Leon Heller
G1HSM

  • 1
  • 2
  • 3
  • 4
  • 5
Total votes: 0

Clawson, I went to the url you provided and looked at the code. Some of it seems to be for using the adc and some seems to be for output to the lcd. I don't really even understand where the data points are loaded. Can someone explain where the raw data goes in, and where the transformed data comes out? Could someone send me a stripped-down version of this code (only containing the actual fft) so that I can start to understand what is going on? Maybe fft is just too complicated for me to understand in asm given my level of experience. Any help would be appreciated.

  • 1
  • 2
  • 3
  • 4
  • 5
Total votes: 0

Quote:

Maybe fft is just too complicated for me to understand in asm

Sounds like it.

 

  • 1
  • 2
  • 3
  • 4
  • 5
Total votes: 0

Here is a very well commented FFT function from the Microchip dsPIC library:

;*********************************************************************
;                                                                    *
;                       Software License Agreement                   *
;                                                                    *
;   The software supplied herewith by Microchip Technology           *
;   Incorporated (the "Company") for its dsPIC controller            *
;   is intended and supplied to you, the Company's customer,         *
;   for use solely and exclusively on Microchip dsPIC                *
;   products. The software is owned by the Company and/or its        *
;   supplier, and is protected under applicable copyright laws. All  *
;   rights are reserved. Any use in violation of the foregoing       *
;   restrictions may subject the user to criminal sanctions under    *
;   applicable laws, as well as to civil liability for the breach of *
;   the terms and conditions of this license.                        *
;                                                                    *
;   THIS SOFTWARE IS PROVIDED IN AN "AS IS" CONDITION.  NO           *
;   WARRANTIES, WHETHER EXPRESS, IMPLIED OR STATUTORY, INCLUDING,    *
;   BUT NOT LIMITED TO, IMPLIED WARRANTIES OF MERCHANTABILITY AND    *
;   FITNESS FOR A PARTICULAR PURPOSE APPLY TO THIS SOFTWARE. THE     *
;   COMPANY SHALL NOT, IN ANY CIRCUMSTANCES, BE LIABLE FOR SPECIAL,  *
;   INCIDENTAL OR CONSEQUENTIAL DAMAGES, FOR ANY REASON WHATSOEVER.  *
;                                                                    *
;   (c) Copyright 2003 Microchip Technology, All rights reserved.    *
;*********************************************************************

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
; This Complex DIF FFT implementation expects that the input data is a
; complex vector such that the magnitude of the real and imaginary parts
; of each of its elements is less than 0.5. If greater or equal to this
; value the results could produce saturation.
;
; Also, the program performs an implicit scaling of 1/2 for every stage
; to the intermediate values, so that the output values are scaled by a
; factor of 1/N, with N the length of the FFT.
;
; NOTE: input is expected in natural ordering, while output is produced
; in bit reverse ordering.
;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


	; Local inclusions.
	.nolist
	.include	"dspcommon.inc"		; fractsetup,CORCON
						; PSVPAG, COEFFS_IN_DATA,
	.list

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

	.section .libdsp, code

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
; _FFTComplexIP: Complex (in-place) DIF FFT.
;
; Operation:
;	F(k) = 1/N*sum_n (f(n)*WN(kn)), WN(kn) = exp[-(j*2*pi*k*n)/N],
;
; n in {0, 1,... , N-1}, and
; k in {0, 1,... , N-1}, with N = 2^m.
;
; Input:
;	w0 = number stages in FFT (log2N)
;	w1 = ptr to complex source vector (srcCV)
;	w2 = ptr to complex twiddle factors (twidFactors)
;	w3 = COEFFS_IN_DATA, or memory program page with twiddle factors.
; Return:
;	w0 = ptr to source vector (srcCV)
;
; System resources usage:
;	{w0..w7}	used, not restored
;	{w8..w13}	saved, used, restored
;	 AccuA		used, not restored
;	 AccuB		used, not restored
;	 CORCON		saved, used, restored
;	 PSVPAG		saved, used, restored (if factors in P memory)
;
; DO and REPEAT instruction usage.
;	2 level DO intruction
;	no REPEAT intructions
;
; Program words (24-bit instructions):
;	60
;
; Cycles (including C-function call and return overheads):
;
;	transform	factors		factors
;	   size		in X-mem	in P-mem
;	-----------------------------------------
;	 32-point	   1664		   1827
;	 64-point	   3802		   4189
;	128-point	   8612		   9511
;	256-point	  19310		  21361
;   512-point     42872       47483
;
;............................................................................

	; Local equates.
	; Commented to avoid IAR assembler problems...
;	.equ	oTwidd,w3
;	.equ	pTwidd,w8
;	.equ	nGrps,w9
;	.equ	pUpper,w10
;	.equ	pLower,w11
;	.equ	groupCntr,w12
;	.equ	buttCntr,w13
	; The symbols have been replaced in the code by the corresponding
	; working registers; the comments maintain the symbols for clarity.

	.global	_FFTComplexIP	; export
_FFTComplexIP:

;............................................................................

	; Save working registers.
	push.d	w8				; {w8,w9} to TOS
	push.d	w10				; {w10,w11} to TOS
	push.d	w12				; {w12,w13} to TOS

;............................................................................

	; Prepare CORCON for fractional computation.
	push	CORCON
	fractsetup	w7

;............................................................................

	; Prepare CORCON and PSVPAG for possible access of data
	; located in program memory, using the PSV.
	push	PSVPAG

	mov	#COEFFS_IN_DATA,w7		; w7 = COEFFS_IN_DATA
	cp	w7,w3				; w7 - w3
	bra	z,_noPSV			; if w3 = COEFFS_IN_DATA
						; no PSV management
						; else
	psvaccess	w7			; enable PSV bit in CORCON
	mov	w3,PSVPAG			; load PSVPAG with program
						; space page offset
_noPSV:

;............................................................................

	push	w1				; save return value (srcCV)

;............................................................................

	; FFT proper.
	mov	#0x1,w9				; to be shifted...
	sl	w9,w0,w9			; w9 = N (1<<log2N)
	mov	#0x1,w3				; initialize twidds offset,
						; also used as num butterflies
	mov	w2,w8				; w8-> WN(0) (real part)

	; Preform all k stages, k = 1:log2N.
	; NOTE: for first stage, 1 butterfly per twiddle factor (w3 = 1)
	; and N/2 groups  (w9 = N/2) for factors WN(0), WN(1), ..., WN(N/2-1).
	mov	w0,w13				; w13= num stages
						; w0 available for reuse
_doStage:

	; Perform all butterflies in each stage, grouped per twiddle factor.

	; Update counter for groups.
	lsr	w9,w9				; w9 = N/(2^k)
	sl	w9,#2,w12			; w12= lower offset
						; nGrps+sizeof(fractcomplex)
						; *2 bytes per element

	; Set pointers to upper "leg" of butterfly:
	mov	w1,w10				; w10-> srcCV (upper leg)

	; Perform all the butterflies in each stage.
	dec	w3,w4				; w4 = butterflies per group
	do	w4,_endBflies		; {	; do 2^(k-1) butterflies

	; Set pointer to lower "leg" of butterfly.
	add	w12,w10,w11			; w11-> srcCV + lower offset
						; (lower leg)
	; Prepare offset for twiddle factors.
	sl	w3,#2,w0			; oTwidd*sizeof(fractcomplex)
						; *2 bytes per element

	; Perform each group of butterflies, one for each twiddle factor.
	dec	w9,w4				; w4 = nGrps-1
	do	w4,_endGroup		; {	; do (nGrps-1)+1 times
						; no more nested do's...
.ifdef YMEM_ERRATA
    nop
.endif
	; Butterfly proper.
	clr	a,[w8]+=2,w6,[w10]+=2,w4	; a  = 0 (for dual fetch only)
						; w6 = Wr
						; pTwidd-> Wi
						; w4 = Ar
						; pUpper-> Ai (=Ci)
	sub	w4,[w11++],w4			; w4 = Ar - Br
						; pLower-> Bi (=Di)
	mpy	w4*w6,a,[w8]-=2,w7,[w10]-=2,w5	; a  = (Ar - Br)*Wr
						; w7 = Wi
						; pTwidd-> Wr
						; w5 = Ai
						; pUpper-> Ar (=Cr)
	; Update twiddle pointer.
	add	w0,w8,w8			; pTwidd-> for next group

	mpy	w4*w7,b				; b  = (Ar - Br)*Wi
	sub	w5,[w11--],w5			; w5 = Ai - Bi
						; pLower-> Br (=Dr)
	msc	w5*w7,a,[w11]+=2,w7		; a  = (Ar - Br)*Wr
						;    - (Ai - Bi)*Wi = Dr
						; w7 = Br
						; pLower-> Bi (=Di)
	mac	w5*w6,b,[w11]-=2,w5		; b  = (Ar - Br)*Wi
						;    + (Ai - Bi)*Wr = Di
						; w5 = Bi
						; pLower-> Br (=Dr)
	sac.r	a,#1,[w11++]			; save 1/2*Dr
						; pLower-> Bi (=Di)
	sac.r	b,#1,[w11++]			; save 1/2*Di
						; pLower-> Br (next)
	lac	[w10++],a			; a  = Ar
						; pUpper-> Ai (=Ci)
	lac	[w10--],b			; b  = Ai
						; pUpper-> Cr
	add	w7,a				; a  = Ar + Br
	add	w5,b				; b  = Ai + Bi
	sac.r	a,#1,[w10++]			; save 1/2*Cr
						; pUpper-> Ci
_endGroup:
	sac.r	b,#1,[w10++]			; save 1/2*Ci
						; pUpper-> Ar (next)
; }

	add	w12,w10,w10			; w10-> upper leg (next set)
_endBflies:
	mov	w2,w8				; rewind twiddle pointer
; }

	; Update offset to factors.
	sl	w3,w3				; oTwidd *= 2

	; Find out whether to perform another stage...
	dec	w13,w13				; w13= log2N - k -1
	bra	gt,_doStage			; if W13 > 0, do next stage
_endStages:					; else, no more stages...

;............................................................................

	pop	w0				; restore return value

;............................................................................

	; Restore PSVPAG and CORCON.
	pop	PSVPAG
	pop	CORCON

;............................................................................

	; Restore working registers.
	pop.d	w12				; {w12,w13} from TOS
	pop.d	w10				; {w10,w11} from TOS
	pop.d	w8				; {w8,w9} from TOS

;............................................................................

	return	

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

	.end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; OEF

It uses the DSP engine, of course. AVRs don't have anything comparable.

Leon

Leon Heller
G1HSM

  • 1
  • 2
  • 3
  • 4
  • 5
Total votes: 0

But you could treat an FFT routine as a black box, you don't have to understand FFTs in intimate detail to be able to use it.

As long as you understand what goes in and what comes out, then how the routine does its trick is not really relevant.

  • 1
  • 2
  • 3
  • 4
  • 5
Total votes: 0

Quote:
But you could treat an FFT routine as a black box

This is exactly what I am looking for. I am just struggling to figure out where the input and output goes in the code I have been looking at.

Leon, thanks. This looks like it might help.

  • 1
  • 2
  • 3
  • 4
  • 5
Total votes: 0

jayjay1974 wrote:
But you could treat an FFT routine as a black box, you don't have to understand FFTs in intimate detail to be able to use it.

As long as you understand what goes in and what comes out, then how the routine does its trick is not really relevant.

Well, both "what goes in" and "what goes out", work better when you actually know what you are doing. Some knowledge about the DFT in general and and the particular FFT algorithm in particular can avoid a lot of tears and hair-pulling.

Stealing Proteus doesn't make you an engineer.

  • 1
  • 2
  • 3
  • 4
  • 5
Total votes: 0

You also need to do things like "windowing", simply applying an FFT won't work.

Leon

Leon Heller
G1HSM

  • 1
  • 2
  • 3
  • 4
  • 5
Total votes: 0

You can find out a LOT more here.

http://www.dspguide.com/pdfbook.htm

It is well written and not constrained to conventional DSPs. All of the algorithms are actually demonstrated in BASIC! I highly recommend it

Jim

Jim Wagner Oregon Research Electronics, Consulting Div. Tangent, OR, USA There are some answers that have no questions.

  • 1
  • 2
  • 3
  • 4
  • 5
Total votes: 0

Thanks, Jim.

  • 1
  • 2
  • 3
  • 4
  • 5
Total votes: 0

Windowing is a technique used when designing digital filters. It is no necessary to do frequency analysis with an FFT. As an aside, no windowing (truncating) is equivalent to applying a rectangular window.

  • 1
  • 2
  • 3
  • 4
  • 5
Total votes: 0

I have do not have much experience with FFTs, but the few times I read things about, windowing the data before applying the FFT is very common.

For example, in LTSpice you can select a quite few windowing options for the FFT graph.

  • 1
  • 2
  • 3
  • 4
  • 5
Total votes: 0

tyler.w.gilbert wrote:
Windowing is a technique used when designing digital filters. It is no necessary to do frequency analysis with an FFT. As an aside, no windowing (truncating) is equivalent to applying a rectangular window.

Windowing is needed to avoid errors due to leakage when performing an FFT!

Leon

Leon Heller
G1HSM

  • 1
  • 2
  • 3
  • 4
  • 5
Total votes: 0

leon_heller wrote:
You also need to do things like "windowing", simply applying an FFT won't work.

Leon

Applying an FFT to a block of data (without windowing) DOES work. The resulting data set IS the FFT of the original data.

Yes, I understand the windowing is necessary and convinient from an analysis point of view (frequency leakage as you say on your other post), but not from a mathematical point of view.

  • 1
  • 2
  • 3
  • 4
  • 5
Total votes: 0

I concede. Windowing will reduce leakage and is necessary. However, collecting 256 data points and performing an FFT is equivalent to a de facto rectangular window. Rectangular windows have the worst frequency leakage characteristics, but for a beginning it is probably sufficient.

  • 1
  • 2
  • 3
  • 4
  • 5
Total votes: 0

Be aware that the frequency span and resolution are intimately tied to the sample rate and the total number of samples.

Jim

Jim Wagner Oregon Research Electronics, Consulting Div. Tangent, OR, USA There are some answers that have no questions.

  • 1
  • 2
  • 3
  • 4
  • 5
Total votes: 0

ganzziani wrote:
Applying an FFT to a block of data (without windowing) DOES work. The resulting data set IS the FFT of the original data.
No, an FFT (there are several) is just an algorithm. The result of an FFT is a set of Fourier coefficients of a discrete Fourier series (DFS), representing a periodic discrete time signal.

The crux lies in the periodic property. The way you obtained (by sampling) the original discrete time signal makes it unlikely that you sampled exactly a multiple of periods of the signal. Instead, you have cut the signal somewhere midway. So your samples don't represent your original signal, but a different one.

You now feed the wrong data into an FFT algorithm, and you obtain the coefficients of a DFS for that wrong signal. Yes, you got the right coefficients for your data, only that your data doesn't represent your signal. So the coefficients aren't the right ones for a DFS of your original signal. Garbage in, garbage out.

Simple windowing, or more advanced windowing techniques like Bartlett or Welsh try to mitigate the problem of feeding the wrong data into an FFT algorithm. The result is a set of DFS coefficients that are in some (but not all) aspects better. "better" as in closer to the real DFS coefficients of your original signal.

Stealing Proteus doesn't make you an engineer.

  • 1
  • 2
  • 3
  • 4
  • 5
Total votes: 0

ArnoldB. I understand and agree to what you are saying.
I was trying to differentiate the analitical vs the mathematical aspects of the DFT.

From a mathematical point of view, if you apply a DFT to a data set, you get a set of DFT coefficients that represent the original data set in the frequency domain.

  • 1
  • 2
  • 3
  • 4
  • 5
Total votes: 0

ganzziani is, of course, correct.

Without windowing, the FFT gives you the spectrum of the data as presented with the assumption that the data presented represents a periodic series with a period (or integer multiple of the period) equal to the sample length.

As previously pointed out, windowing helps to reduce the apparent errors that arise when the assumption in not satisfied. These errors are called "frequency leakage" because they result in frequency components in the output spectrum that are not really present in the input and which would NOT be present if the sampling period is the "correct" duration.

Jim

Jim Wagner Oregon Research Electronics, Consulting Div. Tangent, OR, USA There are some answers that have no questions.

  • 1
  • 2
  • 3
  • 4
  • 5
Total votes: 0

ka7ehk wrote:
You can find out a LOT more here.

http://www.dspguide.com/pdfbook.htm

It is well written and not constrained to conventional DSPs. All of the algorithms are actually demonstrated in BASIC! I highly recommend it

Jim

Amazing site! Thank you for sharing Jim.

  • 1
  • 2
  • 3
  • 4
  • 5
Total votes: 0

I quit trying to implement an FFT algorithm. Instead, I wrote code that computes the regular DFT. It's not perfect yet, but it was a lot less complex than the FFT. Of course, in my case, it is also ~32 times slower. I haven't messed around with windowing yet. I suspect that I will have to do this in the future. Thanks for all of the help.