Go Down

Topic: Arduino FFT code not expected results... (Read 7821 times) previous topic - next topic

Judd_Foster

Alright thanks everyone, I'll try to use an already existing FFT library first, and when I understand it maybe I can try re-writing my own code.

jremington

If you don't want to use the assembly-coded library or the fixed point fft library, here is the floating point code I've been using for some time. It came out of "Numerical Recipes" and was implemented using WinAVR. This version is just a test example with a square wave input and the answer is correct to about 5 decimal places. Note the unusual Fortran-like indexing of the data array (begins with element 1), so to be safe allocate N+1 elements for data.
Code: [Select]
/*
REAL FFT tests.
atmega64 @ 11.059 MHz with 32k external sram
910 milliseconds for 2048 point floating point transform (100 ms/256 points)
// relocated heap to external SRAM via linker command
// -Wl,--defsym=__heap_start=0x801100,--defsym=__heap_end=0x807fff
// also must include floating point printf and lib:
// -Wl,-u,vfprintf -lprintf_flt -lm

*/

#include <avr/io.h>
#include <avr/interrupt.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <avr/pgmspace.h>
#define F_CPU 11059200UL
#include <util/delay.h>      


#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
#define NP 256

///////////////////////////////////////////////////////
// Global Variables

float *data; //allocated at runtime in main
int ndata=NP; //size of transform buffer

// initialization sequence in .init3

void init_extSRAM(void) __attribute__((naked,section(".init3"))) ;

void init_extSRAM(void) {
//Enable external SRAM: set MCUCR, XMCRA and XMCRB for
// External SRAM page configuration: 0000h - 7FFFh
// Lower page wait state(s): None
// Upper page wait state(s): None
MCUCR=0x80;
XMCRA=0x00;    
XMCRB &=~7 ;  //Make bits 0-2 of XMCRB low to enable all of external SRAM
}

volatile unsigned short ticks=0; //0.01 second clock ticks updated in ISR

//  delay for time_ms milliseconds by looping
// _delay_ms() has max delay about 13 ms @ 20 MHz

void delay_ms( unsigned short time_ms )
{
unsigned int i;
for ( i = 0; i < time_ms; i++ )  _delay_ms( 1 );
}

// global COM=0 or 1 specifies UART0 or UART1

int COM=1;

#include "uputget.c"

int uart_putchar(char c, FILE *stream) {
/*
* Send character c to the LCD for use with printf
*/
  uputchar(c);
  return 0;
}

//setup outpuf file for printf

FILE uart_str = FDEV_SETUP_STREAM(uart_putchar, NULL, _FDEV_SETUP_WRITE);

/*
* Timer0 overflow interrupt handler.  
* provides global ticks of 10 ms
*/

ISR(TIMER0_OVF_vect)
{
TCNT0=256-108;
ticks++;
}

void timer_init(void) {

// set up Timer0 to provide 100 tick/sec

TCCR0  = (7<<CS00);  //clock/1024
TIMSK |= (1<<TOIE0); // enable int on overflow
TCNT0=256-108; //exact divisor for 11.059 MHz
}

void get_data(void)
{
int j;
for (j=0;j<NP/2;j++) {
data[j]= -1.0;
data[j+NP/2]= 1.0;
}
}
void four1(int nn, int isign)
{
int n,mmax,m,j,istep,i;
/*double*/
float wtemp,wr,wpr,wpi,wi,theta;

float tempr,tempi;

n=nn << 1;
j=1;
for (i=1;i<n;i+=2) {
if (j > i) {
SWAP(data[j],data[i]);
SWAP(data[j+1],data[i+1]);
}
m=n >> 1;
while (m >= 2 && j > m) {
j -= m;
m >>= 1;
}
j += m;
}
mmax=2;
while (n > mmax) {
istep=mmax << 1;
theta=isign*(6.2831853/mmax);
wtemp=sin(0.5*theta);
wpr = -2.0*wtemp*wtemp;
wpi=sin(theta);
wr=1.0;
wi=0.0;
for (m=1;m<mmax;m+=2) {
for (i=m;i<=n;i+=istep) {
j=i+mmax;
tempr=wr*data[j]-wi*data[j+1];
tempi=wr*data[j+1]+wi*data[j];
data[j]=data[i]-tempr;
data[j+1]=data[i+1]-tempi;
data[i] += tempr;
data[i+1] += tempi;
}
wr=(wtemp=wr)*wpr-wi*wpi+wr;
wi=wi*wpr+wtemp*wpi+wi;
}
mmax=istep;
}
}

void realft(int n, int isign)
{

int i,i1,i2,i3,i4,np3;
float c1=0.5,c2,h1r,h1i,h2r,h2i;
/*double*/
float wr,wi,wpr,wpi,wtemp,theta;

theta=3.141592654/(float) (n>>1);
if (isign == 1) {
c2 = -0.5;
four1(n>>1,1);
} else {
c2=0.5;
theta = -theta;
}
wtemp=sin(0.5*theta);
wpr = -2.0*wtemp*wtemp;
wpi=sin(theta);
wr=1.0+wpr;
wi=wpi;
np3=n+3;
for (i=2;i<=(n>>2);i++) {
i4=1+(i3=np3-(i2=1+(i1=i+i-1)));
h1r=c1*(data[i1]+data[i3]);
h1i=c1*(data[i2]-data[i4]);
h2r = -c2*(data[i2]+data[i4]);
h2i=c2*(data[i1]-data[i3]);
data[i1]=h1r+wr*h2r-wi*h2i;
data[i2]=h1i+wr*h2i+wi*h2r;
data[i3]=h1r-wr*h2r+wi*h2i;
data[i4] = -h1i+wr*h2i+wi*h2r;
wr=(wtemp=wr)*wpr-wi*wpi+wr;
wi=wi*wpr+wtemp*wpi+wi;
}
if (isign == 1) {
data[1] = (h1r=data[1])+data[2];
data[2] = h1r-data[2];
} else {
data[1]=c1*((h1r=data[1])+data[2]);
data[2]=c1*(h1r-data[2]);
four1(n>>1,-1);
}
}

////////////////////////////////////////////////////////////////////////
int main ()
{
int i;
DDRD |= 1;   //set PORTD.0 = output
PORTD &= ~1;   //led on
delay_ms(200);
PORTD |= 1; //led off

uart_init();            
stdout = &uart_str;         // associate stream with stdout
timer_init();
sei();

printf("Real FFT 0.1\n\r");   // announce string

if ( (data = (float *)malloc( (ndata+1)*sizeof(float)) ) == NULL )
printf("Unable to allocate data array\r\n");
else
printf("Successfully allocated data array\r\n");

get_data();

for(i=0; i<32; i+=2)
printf( "%7.2f %7.2f \r\n",(double) data[i],(double) data[i+1]);

ticks=0;
realft(NP,1); //transform
i=ticks;
printf(" realfft in %dx10 ms \r\n",i);
printf(" %d data points \r\n",NP);
for(i=0; i<32; i+=2)
printf( "%7.2f %7.2f \r\n",(double) data[i],(double) data[i+1]);
while(1);
}


Judd_Foster

I was just looking at my code and I realized that when I multiply by sina[0] and sina[4] in the 2 butterflies, those values are equal to zero, so the sine half of the 2 butterflies are all 0 values. Could this be why the cosine values are coming out all weird? Because it's the sine that works but cosine doesn't.
Thanks,
Judd

tmd3


It is very unlikely that anyone on the forum will have the time and/or patience to go through your code line by line ...

And I didn't.  I notice, though, that there's nothing in the code that performs complex arithmetic, and the results of the butterflies are, in general, complex numbers.  It looks like the twiddle factors are complex - the sine and cosine terms are kept separate, at least - but there's no complex multiplication implemented, and all the results are added together as real numbers. 

I don't see anything that suggests that you understand the notion of a complex number.  Try these articles:

  • http://en.wikipedia.org/wiki/Imaginary_number

  • http://en.wikipedia.org/wiki/Complex_number


If the mathematical concepts in those articles aren't immediately familiar, I'd recommend that you step back and work on understanding the math behind the FFT.  It starts with those concepts.  If that's not somewhere you want to go, then I'd recommend sticking with FFT programs and libraries written by others.

el_supremo

Don't send me technical questions via Private Message.

Judd_Foster

I think I might just use someone else's library, I'm not up to par with this math about imaginary and complex numbers...
Thanks for your help anyways,
Judd

Go Up