* Floating point exception
@ 2004-09-11 11:39 Ankit Jain
0 siblings, 0 replies; only message in thread
From: Ankit Jain @ 2004-09-11 11:39 UTC (permalink / raw)
To: gcc, linux prg
#include<stdlib.h>
#include<string.h>
#include<complex.h>
#include<errno.h>
#include<time.h>
#include<stdio.h>
#define NX 4
#define NY 4
#define N 2*NX*NY //N is the size of 2 D DFT
#include <math.h>
#define num 2
#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
void fft(float data[], unsigned long nn[], int ndim,
int isign)
{
int idim;
unsigned long
i1,i2,i3,i2rev,i3rev,ip1,ip2,ip3,ifp1,ifp2;
unsigned long ibit,k1,k2,n,nprev,nrem,ntot;
float tempi,tempr;
double theta,wi,wpi,wpr,wr,wtemp; //Double precision
for trigonometric recur-rences.
for (ntot=1,idim=1;idim<=ndim;idim++) //Compute
total number of complex val-ues.
ntot *= nn[idim];
nprev=1;
for (idim=ndim;idim>=1;idim--)
{ // Main loop over the dimensions.
n=nn[idim];
nrem=ntot/(n*nprev);
ip1=nprev << 1;
ip2=ip1*n;
ip3=ip2*nrem;
i2rev=1;
for (i2=1;i2<=ip2;i2+=ip1)
{ //This is the bit-reversal section of the
routine.
if (i2 < i2rev)
{
for (i1=i2;i1<=i2+ip1-2;i1+=2)
{
for (i3=i1;i3<=ip3;i3+=ip2)
{
i3rev=i2rev+i3-i2;
SWAP(data[i3],data[i3rev]);
SWAP(data[i3+1],data[i3rev+1]);
}
}
}
ibit=ip2 >> 1;
while (ibit >= ip1 && i2rev > ibit)
{
i2rev -= ibit;
ibit >>= 1;
}
i2rev += ibit;
}
//Here begins the Danielson-Lanczos sec-tion of the
routine.
ifp1=ip1;
while (ifp1 < ip2)
{
ifp2=ifp1 << 1;
theta=isign*6.28318530717959/(ifp2/ip1);
//Initialize for the trig. recur-rence.
wtemp=sin(0.5*theta);
wpr = -2.0*wtemp*wtemp;
wpi=sin(theta);
wr=1.0;
wi=0.0;
for (i3=1;i3<=ifp1;i3+=ip1) {
for (i1=i3;i1<=i3+ip1-2;i1+=2) {
for (i2=i1;i2<=ip3;i2+=ifp2) {
k1=i2; //Danielson-Lanczos formula:
k2=k1+ifp1;
tempr=(float)wr*data[k2]-(float)wi*data[k2+1];
tempi=(float)wr*data[k2+1]+(float)wi*data[k2];
data[k2]=data[k1]-tempr;
data[k2+1]=data[k1+1]-tempi;
data[k1] += tempr;
data[k1+1] += tempi;
}
}
wr=(wtemp=wr)*wpr-wi*wpi+wr; //Trigonometric
recurrence.
wi=wi*wpr+wtemp*wpi+wi;
}
ifp1=ifp2;
}
nprev *= n;
}
}
int main()
{
float rin[NX][NY],iin[NX][NY];
float in[2*NX*NY];
unsigned long nn[]={NX,NY};
int i=0,j=0,k=0;
k = 0;
for (i = 0; i< NX; i++)
{
for (j = 0; j< NY; j++)
{
rin[i][j] = (float) (2*i*i);
iin[i][j] = 0;
in[k] = rin[i][j];
k=k+1;
in[k] = iin[i][j];
k= k+1;
}
}
//for (i = 0; i< 32; i++)
printf("%d\t%f\n",i,in[i]);
fft(in,nn,2,1); //Executes the plan
k = 0;
for (i = 0; i< NX; i++)
{
for (j = 0; j< NY; j++)
{
rin[i][j] = in[k];
printf("real part is%f\n",in[k]);
k=k+1;
iin[i][j] = in[k];
printf("imaginary part is%f\n",in[k]);
k=k+1;
}
}
printf("finished fft\n");
return 0;
}
________________________________________________________________________
Yahoo! Messenger - Communicate instantly..."Ping"
your friends today! Download Messenger Now
http://uk.messenger.yahoo.com/download/index.html
^ permalink raw reply [flat|nested] only message in thread
only message in thread, other threads:[~2004-09-11 11:39 UTC | newest]
Thread overview: (only message) (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2004-09-11 11:39 Floating point exception Ankit Jain
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).