Bit Umkehr

[Mary]

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();

}