Hab nun das raus... Hat vlt einer eine Idee wo der Fehler sein kann??
#include<stdio.h>
#include<math.h>
# define pi 3.14
struct komplex
{
double re;
double im;
};
komplex Cadd(komplex x, komplex y)
{
komplex res;
res.re=x.re+y.re;
res.im=x.im+y.im;
return res;
}
komplex Csub(komplex x, komplex y)
{
komplex res;
res.re=x.re-y.re;
res.im=x.im-y.im;
return res;
}
komplex Cmult(komplex x, komplex y)
{
komplex res;
res.re=x.re*y.re-x.im*y.im;
res.im=x.im*y.re+x.re*y.im;
return res;
}
int potenz(int i)
{
int j,k=1;
for(j=0;j<i;j++)
{
k=k*2;
}
return k;
}
int main()
{
int j,n,k,p,r,m,t,M,q=3;
n=potenz(q);
komplex x[n];
//Stützstellen
for(j=0;j<n;j++)
{
x[j].im=0.;
x[j].re=((2*pi*j)/n)-pi;
}
komplex f[n];
/*
//Sägezahn
for (j=0;j<n;j++)
{
if(-pi<=x[j].re<pi)
{
f[j].im=0.;
f[j].re=x[j].re;
}
}
//Rechteck
for (j=0;j<n;j++)
{
f[j].im=0.;
if(x[j].re<0.)
{
f[j].re=pi;
}
else
{
f[j].re=-pi;
}
}
*/
//Dreieck
for (j=0;j<n;j++)
{
f[j].im=0.;
if(x[j].re<-pi/2.)
{
f[j].re=2.*x[j].re+2.*pi;
}
else if(x[j].re>=-pi/2. && x[j].re<pi/2.)
{
f[j].re=-2.*x[j].re;
}
else if(pi/2.<=x[j].re)
{
f[j].re=2.*x[j].re-2.*pi;
}
}
//Umsortierung......DIE IST FALSCH UND ICH WILL WISSEN WIE ES RICHTIG FUNKTIONIERT
for(k=0;k<=n/2-1;k++)
{
t=2*k;
m=n/2;
while(t>=m)
{
t=t-m;
m=m/2;
}
t=t+m;
f[k].re=f[t].re;
printf("%i\n",t);
}
for(k=n/2;k<=n-1;k++)
{
t=2*k+1;
m=n/2;
while(t>=m)
{
t=t-m;
m=m/2;
}
t=t+m;
f[k].re=f[t].re;
printf("%i\n",t);
}
//FFT
for(r=0;r<=q-1;r++)
{
M=potenz(r);
komplex w,x;
w.re=cos(pi/M);
w.im=-sin(pi/M);
for(k=0;k<=M-1;k++)
{
for(j=0;j<=potenz(q-r-1)-1;j++)
{
x=Cmult(w,f[2*j*M+M+k]);
f[2*j*M+M+k]=Csub(f[2*j*M+k],x);
f[2*j*M+k]=Cadd(f[2*j*M+k],x);
}
}
}
//Inverse FFT:Wir machen aus f[j] die inversen: f[j]/n!!!
komplex umkehr;
umkehr.re=1/n;
umkehr.im=0.;
for(j=0;j<=potenz(q)-1;j++)
{
f[j]=Cmult(umkehr,f[j]);
}
//Berechnung der koeffizienten
double A[n/2], B[n/2-1];
for (j=0;j<=n/2;j++)
{
if(j==0) {A[0]=2.*f[0].re;}
else {A[j]=f[j].re+f[n-j].re;
printf("A(%i)=(%3.3f)\n",j,A[j]);
}
}
printf("\n\n");
for (j=1;j<(n/2)-1;j++)
{
B[j]=f[n-1].im-f[j].im;
printf("B(%i)=(%3.3f)\n",j,B[j]);
}
//Auswertung an der Stelle z
double wert,summe,z=pi/2.;
for(j=1;j<=n/2-1;j++)
{
summe=0;
summe=summe+A[j]*cos(j*z)+B[j]*sin(j*z);
}
wert=A[0]/2.+summe+cos(z*n/2.)*A[n/2]/2.;
printf("p(%5.5f)=%5.5f\n",z,wert);
getchar();
}