FFT Radix 3 [DSP]

Prev: Capacity for QAM
Next: Details on ESPS' get_f0 function work

From: Akshaya Mukund on 19 Mar 2010 08:01

I have been trying to code a normal C program for a radix 3 fft problem. My
program works for 3,6,9,27 points but conks off for the 81 point fft
onwards. I would like to know if there is a mistake in my program or is it
because VC can't support such intensive calculations ;P. Well I really
don't know how to go ahead.I thought maybe I had calculated my twiddle
factors wrong and I have even tried writing out the complete equation(81
terms) for one particular output. Nothing wrong with my algo (but I could
definitely be wrong as well). Can anyone please check it out. Here is my
code.

The possible sources of error could be my butterfly function (you can chk
it out below), or maybe the fact that I have messed up some type casting.

If anyone gets lucky with this please let me know. thank you for your
time.

Cheers
Akshaya

/////////////////////////////////////////////////////////////////

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

struct dft
{
double r,i;
};

#define pi (double)2*asin((double)1)

double* data_input(char*, double*); //used to input data
int trittoint(int *trit,int m); //convert trit to int
void tritrev(int *trit, int m); //reverse the trit
void tritconvert(int m, int i,int* trit); //convert to trit
double twiddle_geni(int N,int k,int m); //imag part of twiddle factor
double twiddle_genr(int N,int k,int m); //real part
void butterfly(struct dft *a,int aidx,int stage,int bno,int sets,int N);
//butterfly structure (radix-3)and also twiddle factor mult

int main()
{//to do a 3,6,9,27 just change N to no of pts and m to no. of trits
FILE *fid, *fid1;
int i,j,k,l,N=81,m=4,trit[4],limit,limit1,aidx=0; //3^m = N (no of
trits)
struct dft a[81],b[81],*temp;
char *path = "data2.txt";
double temp1,t2[81],*t1;

//init data for fft
for(i=0;i<N;i++)
{
b[i].r = b[i].i = 0;
a[i].r = 0;
a[i].i = 0;
c[i].r=c[i].i=0;
}

//data entered into input struct var from file
//t1 = data_input(path,t2);
for(i=0;i<N;i++) // i'm entering arbitrary data here
a[i].r = i;

//reordering loop
//reordered into proper order trit conversion then reversal
// this part is fine but you are welcome to chk it out

for(i=1;i<2*N/3;i++)
{
tritconvert(i,m,trit); //convert to trits
tritrev(trit,m); //reverse the trit
j = trittoint(trit,m); //convert to integer

if(j>i) //swap
{
temp1 = a[j].r;
a[j].r = a[i].r;
a[i].r = temp1;
}
}

//fft loop

for(i=0;i<m;i++) //loop takes care of the stage
{
limit = (int)pow((double)3,(double)(m-1-i)); //sets of
//radix 3
butterflies
aidx = 0;
for(j=0;j<limit;j++) //no of sets of butterflies
{
limit1 = (int)pow((double)3,(double)i);
for(k=0;k<limit1;k++) //no of butterflies per set
{
butterfly(a,aidx,i,k,limit,N);
aidx++;
}
aidx += (int)pow((double)3,(double)(i+1))-limit1;
}
}

//writing the fft onto text file
fid = fopen("datafftreal.txt","wb");
fid1 = fopen("datafftimag.txt","wb");
for(i=0;i<27;i++)
{
temp = &(a[i].r);
fwrite(temp,sizeof(double),1,fid);
temp = &(a[i].i);
fwrite(temp,sizeof(double),1,fid1);
}

fclose(fid);
fclose(fid1);

return (0);
}

double* data_input(char* path, double* input)
{
FILE *fid;
int i;

fid = fopen(path,"rb");
fread(input,8,27,fid);

for(i=0;i<27;i++)
{
//svbarve
input[i] = 00.0;
printf("%lf\t",input[i]);
}

fclose(fid);

return (input);
}

void butterfly(struct dft *a,int aidx,int stage,int bno,int sets,int N)
{
int l,k,i,loopinc,n=3;
struct dft temp[3];

loopinc = (int)pow((double)3,(double)stage);

for(l=aidx,i=0;i<3;i++)
{
temp[i].r = a[l].r;
temp[i].i = a[l].i;
a[l].r=0;
a[l].i=0;
l += loopinc;
}

k = aidx;
for(i=0;i<3;i++)
{
for(l=0;l<3;l++)
{
a[k].r += (temp[l].r*twiddle_genr(n,l,i) -
temp[l].i*twiddle_geni(n,l,i))*twiddle_genr(N,l,bno*sets);//doubt
a[k].r -= (temp[l].r*twiddle_geni(n,l,i) +
temp[l].i*twiddle_genr(n,l,i))*twiddle_geni(N,l,bno*sets);
a[k].i += (temp[l].r*twiddle_geni(n,l,i) +
temp[l].i*twiddle_genr(n,l,i))*twiddle_genr(N,l,bno*sets);//doubt
a[k].i += (temp[l].r*twiddle_genr(n,l,i) -
temp[l].i*twiddle_geni(n,l,i))*twiddle_geni(N,l,bno*sets);
}
k += loopinc;
}
}

void tritconvert(int i, int m, int* trit)
//used inreordering loop (conv to trit)
{
int j=0;
for(j=0;j<m;j++)
trit[j]=0;

j=0;
while(i>0)
{
trit[j]=i%3;
i = i/3;
j++;
}
}

void tritrev(int *trit,int m) //used in reordering loop (rev the trits)
{
int temp,i,j;
for(i=0,j=m-1;i<m/2;i++,j--)
{
temp = trit[i];
trit[i] = trit[j];
trit[j] = temp;
}
}

int trittoint(int *trit,int m) //used in reordering loop(convert to int)
{
int i,j=0,z=1;
for(i=0;i<m;i++)
{
j = j + trit[i]*z;
z*= 3;
}
return (j);
}

double twiddle_genr(int N,int k,int m)
{
double te;
te = cos((-2*k*m*pi)/N);
return(te);
}
double twiddle_geni(int N,int k,int m)
{
double te;
te = sin((-2*k*m*pi)/N);
return(te);
}

From: Clay on 19 Mar 2010 16:44

On Mar 19, 8:01 am, "Akshaya Mukund"
<akshaya.mukund(a)n_o_s_p_a_m.gmail.com> wrote:
> I have been trying to code a normal C program for a radix 3 fft problem. My
> program works for 3,6,9,27 points but conks off for the 81 point fft
> onwards. I would like to know if there is a mistake in my program or is it
> because VC can't support such intensive calculations ;P. Well I really
> don't know how to go ahead.I thought maybe I had calculated my twiddle
> factors wrong and I have even tried writing out the complete equation(81
> terms) for one particular output. Nothing wrong with my algo (but I could
> definitely be wrong as well). Can anyone please check it out. Here is my
> code.
>
> The possible sources of error could be my butterfly function (you can chk
> it out below), or maybe the fact that I have messed up some type casting.
>
> If anyone gets lucky with this please let me know. thank you for your
> time.
>
> Cheers
> Akshaya
>
> /////////////////////////////////////////////////////////////////
>
> #include <stdio.h>
> #include <stdlib.h>
> #include <math.h>
>
> struct dft
> {
> double r,i;
>
> };
>
> #define pi (double)2*asin((double)1)
>
> double* data_input(char*, double*); //used to input data
> int trittoint(int *trit,int m); //convert trit to int
> void tritrev(int *trit, int m); //reverse the trit
> void tritconvert(int m, int i,int* trit); //convert to trit
> double twiddle_geni(int N,int k,int m); //imag part of twiddle factor
> double twiddle_genr(int N,int k,int m); //real part
> void butterfly(struct dft *a,int aidx,int stage,int bno,int sets,int N);
> //butterfly structure (radix-3)and also twiddle factor mult
>
> int main()
> {//to do a 3,6,9,27 just change N to no of pts and m to no. of trits
> FILE *fid, *fid1;
> int i,j,k,l,N=81,m=4,trit[4],limit,limit1,aidx=0; //3^m = N (no of
> trits)
> struct dft a[81],b[81],*temp;
> char *path = "data2.txt";
> double temp1,t2[81],*t1;
>
> //init data for fft
> for(i=0;i<N;i++)
> {
> b[i].r = b[i].i = 0;
> a[i].r = 0;
> a[i].i = 0;
> c[i].r=c[i].i=0;
> }
>
> //data entered into input struct var from file
> //t1 = data_input(path,t2);
> for(i=0;i<N;i++) // i'm entering arbitrary data here
> a[i].r = i;
>
> //reordering loop
> //reordered into proper order trit conversion then reversal
> // this part is fine but you are welcome to chk it out
>
> for(i=1;i<2*N/3;i++)
> {
> tritconvert(i,m,trit); //convert to trits
> tritrev(trit,m); //reverse the trit
> j = trittoint(trit,m); //convert to integer
>
> if(j>i) //swap
> {
> temp1 = a[j].r;
> a[j].r = a[i].r;
> a[i].r = temp1;
> }
> }
>
> //fft loop
>
> for(i=0;i<m;i++) //loop takes care of the stage
> {
> limit = (int)pow((double)3,(double)(m-1-i)); //sets of
> //radix 3
> butterflies
> aidx = 0;
> for(j=0;j<limit;j++) //no of sets of butterflies
> {
> limit1 = (int)pow((double)3,(double)i);
> for(k=0;k<limit1;k++) //no of butterflies per set
> {
> butterfly(a,aidx,i,k,limit,N);
> aidx++;
> }
> aidx += (int)pow((double)3,(double)(i+1))-limit1;
> }
> }
>
> //writing the fft onto text file
> fid = fopen("datafftreal.txt","wb");
> fid1 = fopen("datafftimag.txt","wb");
> for(i=0;i<27;i++)
> {
> temp = &(a[i].r);
> fwrite(temp,sizeof(double),1,fid);
> temp = &(a[i].i);
> fwrite(temp,sizeof(double),1,fid1);
> }
>
> fclose(fid);
> fclose(fid1);
>
> return (0);
>
> }
>
> double* data_input(char* path, double* input)
> {
> FILE *fid;
> int i;
>
> fid = fopen(path,"rb");
> fread(input,8,27,fid);
>
> for(i=0;i<27;i++)
> {
> //svbarve
> input[i] = 00.0;
> printf("%lf\t",input[i]);
> }
>
> fclose(fid);
>
> return (input);
>
> }
>
> void butterfly(struct dft *a,int aidx,int stage,int bno,int sets,int N)
> {
> int l,k,i,loopinc,n=3;
> struct dft temp[3];
>
> loopinc = (int)pow((double)3,(double)stage);
>
> for(l=aidx,i=0;i<3;i++)
> {
> temp[i].r = a[l].r;
> temp[i].i = a[l].i;
> a[l].r=0;
> a[l].i=0;
> l += loopinc;
> }
>
> k = aidx;
> for(i=0;i<3;i++)
> {
> for(l=0;l<3;l++)
> {
> a[k].r += (temp[l].r*twiddle_genr(n,l,i) -
> temp[l].i*twiddle_geni(n,l,i))*twiddle_genr(N,l,bno*sets);//doubt
> a[k].r -= (temp[l].r*twiddle_geni(n,l,i) +
> temp[l].i*twiddle_genr(n,l,i))*twiddle_geni(N,l,bno*sets);
> a[k].i += (temp[l].r*twiddle_geni(n,l,i) +
> temp[l].i*twiddle_genr(n,l,i))*twiddle_genr(N,l,bno*sets);//doubt
> a[k].i += (temp[l].r*twiddle_genr(n,l,i) -
> temp[l].i*twiddle_geni(n,l,i))*twiddle_geni(N,l,bno*sets);
> }
> k += loopinc;
> }
>
> }
>
> void tritconvert(int i, int m, int* trit)
> //used inreordering loop (conv to trit)
> {
> int j=0;
> for(j=0;j<m;j++)
> trit[j]=0;
>
> j=0;
> while(i>0)
> {
> trit[j]=i%3;
> i = i/3;
> j++;
> }
>
> }
>
> void tritrev(int *trit,int m) //used in reordering loop (rev the trits)
> {
> int temp,i,j;
> for(i=0,j=m-1;i<m/2;i++,j--)
> {
> temp = trit[i];
> trit[i] = trit[j];
> trit[j] = temp;
> }
>
> }
>
> int trittoint(int *trit,int m) //used in reordering loop(convert to int)
> {
> int i,j=0,z=1;
> for(i=0;i<m;i++)
> {
> j = j + trit[i]*z;
> z*= 3;
> }
> return (j);
>
> }
>
> double twiddle_genr(int N,int k,int m)
> {
> double te;
> te = cos((-2*k*m*pi)/N);
> return(te);}
>
> double twiddle_geni(int N,int k,int m)
> {
> double te;
> te = sin((-2*k*m*pi)/N);
> return(te);
>
>
>
> } - Hide quoted text -
>
> - Show quoted text -

I'm sure the problem is in your code and not a limitation of VC.

Clay

From: Jerry Avins on 19 Mar 2010 17:02

Akshaya Mukund wrote:
> I have been trying to code a normal C program for a radix 3 fft problem. My
> program works for 3,6,9,27 points but conks off for the 81 point fft
> onwards. I would like to know if there is a mistake in my program or is it
> because VC can't support such intensive calculations ;P. Well I really
> don't know how to go ahead.I thought maybe I had calculated my twiddle
> factors wrong and I have even tried writing out the complete equation(81
> terms) for one particular output. Nothing wrong with my algo (but I could
> definitely be wrong as well). Can anyone please check it out. Here is my
> code.
>
> The possible sources of error could be my butterfly function (you can chk
> it out below), or maybe the fact that I have messed up some type casting.
>
> If anyone gets lucky with this please let me know. thank you for your
> time.

Do you use signed bytes as counters? 3*81 wraps.

Jerry
--
Discovery consists of seeing what everybody has seen, and thinking what
nobody has thought. .. Albert Szent-Gyorgi
��

|
Pages: 1
Prev: Capacity for QAM
Next: Details on ESPS' get_f0 function work