fft.c

FFT的实现在LPC2378上。首先采用了拉个朗日插值算法对采集的电网波形频率跟踪

//-----------------------------------------------------------------------------
// IntFFT_BRIN.c
//-----------------------------------------------------------------------------
// Copyright 2003 Cygnal Integrated Products, Inc.
//
// AUTH: BD
// DATE: 30 JAN 03
//
// This program collects data using ADC0 at <SAMPLE_RATE> Hz and performs
// an FFT on the data. The Real and Imaginary parts of the results are then
// sent to the UART peripheral at <BAUDRATE> bps, where they can be displayed
// or captured using a terminal program.
//
// Note that the FFT performed in this software is optimized for storage space
// (RAM). The resulting Frequency-domain data is not suitable for analyzing
// Signal-to-noise or distortion performance.
//
// This program uses a 22.1184 MHz crystal oscillator multiplied by (9/4)
// for an effective SYSCLK of 49.7664 Mhz. This program also initializes and
// uses UART0 at <BAUDRATE> bits per second.
//
// Target: C8051F12x
// Tool chain: KEIL C51 6.03
//
//-----------------------------------------------------------------------------
// Includes
//-----------------------------------------------------------------------------
#include <stdio.h>
#include "FFT_Code_Tables.h" // Code Tables for FFT routines

//-----------------------------------------------------------------------------
// Global CONSTANTS and Variable Type Definitions
//-----------------------------------------------------------------------------
#define NUM_BITS 16 // Number of Bits in Data

//-----------------------------------------------------------------------------
// Function PROTOTYPES
//-----------------------------------------------------------------------------
void WindowCalc(INT16 Win_Array[], unsigned char SE_data);
void Int_FFT(INT16 ReArray[], INT16 ImArray[]);
void Bit_Reverse(INT16 BR_Array[]);

//-----------------------------------------------------------------------------
// MAIN Routine
//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------
// MAIN Routine
//-----------------------------------------------------------------------------
/*void main()
{
	// disable watchdog timer
	WDTCN = 0xde;
	WDTCN = 0xad;
	SYSCLK_Init(); // initialize external clock and PLL
	PORT_Init (); // set up Port I/O
	UART0_Init (); // initialize UART0
	TIMER3_Init (SYSCLK/SAMPLE_RATE); // initialize Timer3 to overflow at
	// <SAMPLE_RATE>
	ADC0_Init (); // init ADC0
	EA = 1; // globally enable interrupts
	while (1)
	{
		ADC_Index = 0;
		Conversion_Set_Complete = 0;
		EIE2 |= 0x02; // enable ADC interrupts
		SFRPAGE = LEGACY_PAGE;
		while(!Conversion_Set_Complete);
		SFRPAGE = UART0_PAGE;
		printf(“\nCollected Data\nSample\tValue\n”);
		for (BinNum = 0; BinNum < NUM_FFT; BinNum++)
		{
			// Print Data in the format: Sample <tab> Value <tab>
			printf(“%d\t%u\n”, BinNum, Real[BinNum]);
		}
		WindowCalc(Real, 1); // Window Real Data, and convert to
		// differential if it is single-ended
		Bit_Reverse(Real); // Sort Real (Input) Data in bit-reverse
		// order
		Int_FFT(Real, Imag); // Perform FFT on data
		SFRPAGE = UART0_PAGE;
		printf(“\nBin\tReal\tImag\n”);
		// Output the FFT data to the UART
		for (BinNum = 0; BinNum < NUM_FFT; BinNum++)
		{
			// Print Data in the format: Bin <tab> Real <tab> Imaginary
			printf(“%d\t%d\t%d\n”, BinNum, Real[BinNum], Imag[BinNum]);
		}
		if (RUN_ONCE)
			while(1);
	}
} // END MAIN
*/
//-----------------------------------------------------------------------------
// WindowCalc
//-----------------------------------------------------------------------------
//
// Uses the values in WindowFunc[] to window the stored data.
//
// The WindowFunc[] Array contains window coefficients between samples
// 0 and (NUM_FFT/2)-1, and samples from NUM_FFT/2 to NUM_FFT-1 are the mirror
// image of the other side.
// Window values are interpreted as a fraction of 1 (WindowFunc[x]/65536).
// The value at NUM_FFT/2 is assumed to be 1.0 (65536).
//
// If SE_data = 1, the input data is assumed to be single-ended, and is
// converted to a 2’s complement, differential representation, to cancel the DC
// offset.
//
void WindowCalc(INT16 Win_Array[], unsigned char SE_data)
{
    INT16 index;
	#if (WINDOW_TYPE != 0) // Use this section if a window has been specified
	IBALONG NewVal;
	if (SE_data) // If data is single-ended,
		Win_Array[0] ^= 0x8000; // convert it to differential
	NewVal.l = (long)Win_Array[0] * WindowFunc[0];
	if ((NewVal.l < 0)&&(NewVal.i[0]))
		Win_Array[0] = NewVal.i[1] + 1;
	else Win_Array[0] = NewVal.i[1];
	if (SE_data) // If data is single-ended,
		Win_Array[NUM_FFT/2] ^= 0x8000; // convert it to differential
	for (index = 1; index < NUM_FFT/2; index++)
	{
		// Array positions 1 to (NUM_FFT/2 - 1)
		if (SE_data) // If data is single-ended,
			Win_Array[index] ^= 0x8000; // convert it to differential
		NewVal.l = (long)Win_Array[index] * (long)WindowFunc[index];
		if ((NewVal.l < 0)&&(NewVal.i[0]))
			Win_Array[index] = NewVal.i[1] + 1;
		else Win_Array[index] = NewVal.i[1];
		// Array positions (NUM_FFT/2 + 1) to (NUM_FFT - 1)
		if (SE_data) // If data is single-ended,
			Win_Array[NUM_FFT-index] ^= 0x8000; // convert it to differential
		NewVal.l = (long)Win_Array[NUM_FFT-index] * (long)WindowFunc[index];
		if ((NewVal.l < 0)&&(NewVal.i[0]))
			Win_Array[NUM_FFT-index] = NewVal.i[1] + 1;
		else Win_Array[NUM_FFT-index] = NewVal.i[1];
	}
	#endif
	#if (WINDOW_TYPE == 0) // Compile this if no window has been specified
	if (SE_data) // If data is single-ended,
	{ // convert it to differential
		for (index = 0; index < NUM_FFT; index++)
		{
			Win_Array[index] ^= 0x8000; // XOR MSB with ‘1’ to invert
		}
	}
	#endif
} // END WindowCalc
	//-----------------------------------------------------------------------------
// Bit_Reverse
//-----------------------------------------------------------------------------
//
// Sorts data in Bit Reversed Address order
//
// The BRTable[] array is used to find which values must be swapped. Only
// half of this array is stored, to save code space. The second half is
// assumed to be a mirror image of the first half.
//
void Bit_Reverse(INT16 BR_Array[])
{
	#if (NUM_FFT >= 512)
		WORD swapA, swapB, sw_cnt; // Swap Indices
	#endif
	#if (NUM_FFT <= 256)
	unsigned char swapA, swapB, sw_cnt; // Swap Indices
	#endif
	INT16 TempStore;
	// Loop through locations to swap
	for (sw_cnt = 1; sw_cnt < NUM_FFT/2; sw_cnt++)
	{
		swapA = sw_cnt; // Store current location
		swapB = BRTable[sw_cnt] * 2; // Retrieve bit-reversed index
		if (swapB > swapA) // If the bit-reversed index is
		{ // larger than the current index,
			TempStore = BR_Array[swapA]; // the two data locations are
			BR_Array[swapA] = BR_Array[swapB]; // swapped. Using this comparison
			BR_Array[swapB] = TempStore; // ensures that locations are only
		} // swapped once, and never with
		// themselves
		swapA += NUM_FFT/2; // Now perform the same operations
		swapB++; // on the second half of the data
		if (swapB > swapA)
		{
			TempStore = BR_Array[swapA];
			BR_Array[swapA] = BR_Array[swapB];
		BR_Array[swapB] = TempStore;
		}
	}
} // END Bit Reverse Order Sort


//-----------------------------------------------------------------------------
// Int_FFT
//-----------------------------------------------------------------------------
//
// Performs a Radix-2 Decimation-In-Time FFT on the input array ReArray[]
//
// During each stage of the FFT, the values are calculated using a set of
// “Butterfly” equations, as listed below:
//
// Re1 = Re1 + (Cos(x)*Re2 + Sin(x)*Im2)
// Re2 = Re1 - (Cos(x)*Re2 + Sin(x)*Im2)
// Im1 = Im1 + (Cos(x)*Im2 - Sin(x)*Re2)
// Im2 = Im1 - (Cos(x)*Im2 - Sin(x)*Re2)
//
// The routine implements this calculation using the following values: