Fft

Fast Fourier Transform Computation Using a C8051 MCU

Fan Wang Huazhong University of Science & Technology

Photo of the System

Abstract ☼ A method of calculating FFT using a C8051 MCU is

introduced in this paper. The software is aimed at computational efficiency and accuracy. ☼ By three key techniques: (1)avoidance of multiplication whenever possible; (2)16-bit integer storage; (3)on-chip PLL, a real-time audio spectrum is achieved. ☼ The results are displayed on a LCD screen and carefully analyzed. Basic theories of FFT and common problems in practice are discussed in detail.

BACKGROUNDS

A. B. C. D.

Discrete Fourier Transform and Fast Fourier Transform FFT Bit-Reversal Algorithm Anti-aliasing Filtering and Windowing Function Overview of the C8051 MCU

A. Discrete Fourier Transform and Fast Fourier Transform

☼ For a length N complex sequence x(n), n=0,1,2,…N-1, the discrete Fourier transform(DFT) is defined by N −1

X (k ) = ∑ x(n)WNnk ,

k = 0,1,2,...N − 1

W N = e − j 2π / N

where

n= 0

,

☼ According to periodicity and symmetry k. N

X ( k + N / 2) = G ( k ) − W H ( k )

X (k ) = G (k ) + WNk H (k ) where

G (k ) =

k = 0,1,2,...N / 2 − 1 k = 0,1,2,...N / 2 − 1

N −1 2

∑ x(2n)WNkn/ 2 n =0

H (k ) =

N −1 2

∑ x(2n + 1)W n =0

nk N /2

A. Discrete Fourier Transform and Fast Fourier Transform We have the FFT process as the following:

X (0) = G (0) + W40 H (0)

X (1) = G (1) + W41 H (1) X (2) = G (0) − W40 H (0)

X (3) = G (1) − W41 H (1) (N=4)

(N=8)

B.FFT Bit-Reversal Algorithm

Previous Data

Previous Index

New Index

Desired Data

x(0)

000

000

x(0)

x(1)

001

100

x(4)

x(2)

010

010

x(2)

x(3)

011

110

x(6)

x(4)

100

001

x(1)

x(5)

101

101

x(5)

x(6)

110

011

x(3)

x(7)

111

111

x(7)

C. Anti-aliasing Filtering and Windowing ☼ Anti-aliasing Filtering

Solution: A passive low-pass filter.

C. Anti-aliasing Filtering and Windowing ☼ Windowing

C. Anti-aliasing Filtering and Windowing ☼ Windowing

Type

Equation

0: None(Rectangular)

W(n)=1

1: Triangular

W(n)=n/(N/2) (0≤n≤N/2) W(n)=2-n/(N/2) (N/2
2: Hanning

W(n)=0.5-0.5cos(2πn/N)

3: Hamming

W(n)=0.54-0.46cos(2πn/N)

4: Blackman

W(n) = 0.42 - 0.5cos(2πn / N) +0.08cos(4πn / N)

D. Function Overview of the C8051 MCU ☼ A high performance C8051 MCU by Silicon Laboratories™— C8051F120. ☼ Rich analog on-chip sources a 12-bit, 8channel, 100kHz sample rate ADC, two 12-bit DACs, two analog comparators, a voltage reference, etc.

☼ On-chip PLL, the peak speed of the chip can reach as high as 100MIPS. ☼ A 2-cycle 16×16 MAC engine accelerates the computation

☼ Pipelined instruction architecture 8051 core enables it to execute most of instruction set in only 1 or 2 system clocks.

SOFTWARE DESIGN

A. ADC0_ISR( ) B. WindowCalc( ) C. Bit_Reverse( ) D. Int_FFT( ) E. Spectrum_Display( )

A. ADC0_ISR( )

☼ Conversation is initiated when a Timer3 overflow happens. ☼ The data is then stored in the ADC0H :ADC0L registers. ☼ The overflow interval of Timer3 is determined by a configurable macro variable SAMPLE_RATE .

B. WindowCalc( )

☼ Pre-calculated The window coefficients are pre-calculated and stored in the header file FFT_Code_Tables.h. The type define variable WINDOW_TYPE can be selected from 0 to 4.

☼ Single-ended to differential The function also converts the single-ended data collected by ADC0 into differential, in order to remove the DC component of the input signal. The process is carried out by XOR the input data with 0x8000, which has the same effect with subtracting 32768 from the input data, so that the data ranges from 32752 to -32768.

C. Bit_Reverse( )

☼ To save data storage space, only half of the bit-reverse array is stored, the other half is useless for the data will be exchanged by the first half of operation.

☼ For example, an 8-point FFT may use the following bit-reverse array: BRTable[ ]={0,2,1,3}; And the sentence NewIndex=BRTable[PreviousIndex]*2 will help us location the nth reversed index.

Start

D. Int_FFT( )

Im2 = Im1 - (cos(x) × Im2 - sin(x) × Re2)

Butterfly loop

Group loop

Re2 = Re1 - (cos(x) × Re2 + sin(x) × Im2) Im1 = Im1 + (cos(x) × Im2 - sin(x) × Re2)

Stage loop

Re1 = Re1 + (cos(x) × Re2 + sin(x) × Im2)

Stage =0 Group =0 Butterfly =0

Butterfly Computation

Butterfly +=1

Butterfly finished ?

N

N

Y

where

Re1 = ReArray[indexA] Re2 = ReArray[indexB] Im1 = ImArray[indexA] Im2 = ImArray[indexB] x =2×pi×(sin_index/NUM_FFT), in radians.

Group +=1 Butterfly =0

Group Finished ? Y Stage +=1 Group =0

Stage Finished ?

Y End

N

E. Spectrum_Display( )

The FFT result is in complex form, of which the real part is stored in Real[], and the imaginary part is stored in Imag[]. Based on the two arrays, the magnitude spectrum and phase spectrum can be easily obtained. This software only calculates power spectrum of the input signal. The spectrum is displayed on a 128×64 point LCD screen.

EXPERIMENT RESULTS AND DISCUSSION

A. Function Change B. Software Configuration Change C. Property of Symmetry

A. Function Change ☼ Magnitude and Frequency Change of a Sinusoidal Signal

Fig.TR.1. A 2.5kHz 1.3V sine.

Fig.TR.3 A 2.5kHz 2.0V sine.

Fig.TR.2 A 2.5kHz 2.0V sine.

Fig.TR.4 A 5kHz 2.0V sine.

A. Function Change

☼ Wave Form Change

Fig.TR.5. Triangular signal.

Fig.TR.6. Rectangular signal.

B. Software Configuration Change

Name of variable

Function

NUM_FFT

Number of points of FFT algorithm

SAMPLE_RATE

Sampling rate of ADC0 in Hz

WINDOW_TYPE

Type of windows

RUN_ONCE

Program runs once when the value is 0,many times when the value is non-zero.

B. Software Configuration Change

☼ NUM_FFT Change

Fig.TR.7. A 256-point square wave

Fig.TR.8. A 128-point square wave

B. Software Configuration Change

☼ WINDOW_TYPE Change

Fig.TR.9. Rectangular without windowing.

Fig.TR.10.Rectangular with Blackman window.

B. Software Configuration Change

☼ SAMPLE_RATE Change

Fig.TR.11. 5k sine at 40kHz sampling rate.

Fig.TR.12. 5k sine at 20kHz sampling rate.

B. Software Configuration Change

☼ RUN_ONCE Change

To check the speed of the system, RUN_ONCE is set to zero. For all audio signals at a sampling rate of 40kHz, the LCD screen never blinks and a steady frequency spectrum of the input time domain signals can be read without any obstacles. This proves the real-time property of the system.

C. Property of Symmetry The real part of FFT output exhibits even symmetry and imaginary part of FFT output exhibits odd symmetry. This results in the symmetry of frequency spectrum. That is, for N-point FFT products, point N/2+1 is the same as point N/2-1, point N/2+2 is the same as point N/2-2, etc. To illustrate this property, set the full range of the screen to 40kHz. Here are the results.

Fig.TR.13. 2kHz square wave at full range.

Fig.TR.14. 13kHz square wave at full range.

REFERENCES [1] Edward W. Kanmen and Bonnie S. Heck, Fundamentals of Signals and Systems: Using the Web and MATLAB, 2nd ed. Pearson Education, 2002. [2] Julius O. Smith Ⅲ, Mathematics of the Discrete Fourier Transform (DFT) with Audio Applications, 2nd ed. http://ccrma.stanford.edu/~jos/mdft/ [3] Junli Zheng, Signals and systems, 2nd ed. Higher Education Press, 2000. [4] AN142 “FFT Routines for the C8051F12X Family”, Rev.1.1. Inc. ©Silicon Laboratories, 2003. [5] “FFT Library-Module user’s Guide C28x Foundation Software”. Inc. ©Texas Instruments, 2002.

Thank you!