Hi All,
I am trying to implement a Fast Fourier Transform, followed by an
Inverse FFT on the iPhone in Objective-C. To do this I am using the
functions within the 'rvfft.cpp' file available here:
http://www.musicdsp.org/archive.php?classid=2#79
I am trying to perform the operation on an array of doubles of length
131072 (padded with zeros to ensure it is divisible by 2). My code is
below, and the operation seems to be working, Theoretically performing
the IFFT directly on the output of the FFT (with no additional
processing) should reproduce the original signal (obviously not in
infinite precision detail, but this is irrelevant at the moment).
However, the output (which I have written to a WAV file using a
function i know to work correctly) contains a LOT of noise. You can
hear the content from the original array of doubles (which is a sine
wave sweep), but under the sound of what appears to be noise.
My Code is below, I wonder if any of you kind people have any
suggestions?
double * CFArray3;
CFArray3 = (double *)malloc((CFArray1Size)* sizeof(double));
CFArray3 = CFArray1;
long size = CFArray1Size;
//FFT Algo's
//TODO: Currently introduce noise?!
realfft_radix2(CFArray3,CFArray1Size);
//IFFT
irealfft_split(CFArray3,CFArray1Size);
The realfft_radix2 function:
void realfft_radix2(double *data,long n){
float pi = 3.141592653589793238462643;
double xt,a,e, t1, t2, cc, ss;
long i, j, k, n1, n2, n3, n4, i1, i2, i3, i4;
n4=n-1;
//data shuffling
for (i=0,j=0,n2=n/2; i<n4 ; i++){
if (i<j){
xt=data[j];
data[j]=data[i];
data[i]=xt;
}
k=n2;
while (k<=j){
j-=k;
k>>=1;
}
j+=k;
}
/* -------------------- */
for (i=0; i<n; i += 2)
{
xt = data[i];
data[i] = xt + data[i+1];
data[i+1] = xt - data[i+1];
}
/* ------------------------ */
n2 = 1;
for (k=n;k>2;k>>=1){
n4 = n2;
n2 = n4 << 1;
n1 = n2 << 1;
e = 2*pi/(n1);
for (i=0; i<n; i+=n1){
xt = data[i];
data[i] = xt + data[i+n2];
data[i+n2] = xt-data[i+n2];
data[i+n4+n2] = -data[i+n4+n2];
a = e;
n3=n4-1;
for (j = 1; j <=n3; j++){
i1 = i+j;
i2 = i - j + n2;
i3 = i1 + n2;
i4 = i - j + n1;
cc = cos(a);
ss = sin(a);
a += e;
t1 = data[i3] * cc + data[i4] * ss;
t2 = data[i3] * ss - data[i4] * cc;
data[i4] = data[i2] - t2;
data[i3] = -data[i2] - t2;
data[i2] = data[i1] - t1;
data[i1] += t1;
}
}
}
//division with array length
for(i=0;i<n;i++) data[i]/=n;
}
The irealfft_split function:
void irealfft_split(double *data,long n){
long i,j,k,i5,i6,i7,i8,i0,idd,i1,i2,i3,i4,n2,n4,n8,n1;
double t1,t2,t3,t4,t5,a3,ss1,ss3,cc1,cc3,a,e,sqrt2;
float pi = 3.141592653589793238462643;
sqrt2=sqrt(2.0);
n1=n-1;
n2=n<<1;
for(k=n;k>2;k>>=1){
idd=n2;
n2>>=1;
n4=n2>>2;
n8=n2>>3;
e = 2*pi/(n2);
i1=0;
do{
for (; i1<n; i1+=idd){
i2=i1+n4;
i3=i2+n4;
i4=i3+n4;
t1=data[i1]-data[i3];
data[i1]+=data[i3];
data[i2]*=2;
data[i3]=t1-2*data[i4];
data[i4]=t1+2*data[i4];
if (n4!=1){
i0=i1+n8;
i2+=n8;
i3+=n8;
i4+=n8;
t1=(data[i2]-data[i0])/sqrt2;
t2=(data[i4]+data[i3])/sqrt2;
data[i0]+=data[i2];
data[i2]=data[i4]-data[i3];
data[i3]=2*(-t2-t1);
data[i4]=2*(-t2+t1);
}
}
idd<<=1;
i1=idd-n2;
idd<<=1;
} while ( i1<n1 );
a=e;
for (j=2; j<=n8; j++){
a3=3*a;
cc1=cos(a);
ss1=sin(a);
cc3=cos(a3);
ss3=sin(a3);
a=j*e;
i=0;
idd=n2<<1;
do{
for (; i<n; i+=idd){
i1=i+j-1;
i2=i1+n4;
i3=i2+n4;
i4=i3+n4;
i5=i+n4-j+1;
i6=i5+n4;
i7=i6+n4;
i8=i7+n4;
t1=data[i1]-data[i6];
data[i1]+=data[i6];
t2=data[i5]-data[i2];
data[i5]+=data[i2];
t3=data[i8]+data[i3];
data[i6]=data[i8]-data[i3];
t4=data[i4]+data[i7];
data[i2]=data[i4]-data[i7];
t5=t1-t4;
t1+=t4;
t4=t2-t3;
t2+=t3;
data[i3]=t5*cc1+t4*ss1;
data[i7]=-t4*cc1+t5*ss1;
data[i4]=t1*cc3-t2*ss3;
data[i8]=t2*cc3+t1*ss3;
}
idd<<=1;
i=idd-n2;
idd<<=1;
} while(i<n1);
}
}
/*----------------------*/
i0=0;
idd=4;
do{
for (; i0<n1; i0+=idd){
i1=i0+1;
t1=data[i0];
data[i0]=t1+data[i1];
data[i1]=t1-data[i1];
}
idd<<=1;
i0=idd-2;
idd<<=1;
} while ( i0<n1 );
/*----------------------*/
//data shuffling
for (i=0,j=0,n2=n/2; i<n1 ; i++){
if (i<j){
t1=data[j];
data[j]=data[i];
data[i]=t1;
}
k=n2;
while (k<=j){
j-=k;
k>>=1;
}
j+=k;
}
}