Probleme mit diskreter Fouriertransformation



  • Hallo,

    ich versuche gerade eine diskrete FourierTransformation (nicht FFT) zu implementieren und ....scheitere. Ich denke es hat was mit den komplexen Zahlen zu tun, aber tatsächlich sind dies auch meine ersten oder zweiten Schritte in C.
    Könnte vielleicht jemand darauf schauen. Das wäre äußerst nett. DankeDanke Max

    /*
     * main.c
     *
     *  Created on: 05.10.2012
     *      Author: herrk
     */
    
    #include<stdio.h>
    #include<math.h>
    #include<complex.h>
    
    /*****************Create Structur for Complex Numbers**************/
    typedef struct {double re; double im;} Complex;
    
    /***************Data in complex Double form gets stored in structure*************/
    void ComplexSplit (complex double* In, Complex* Out, int N){
    	int i;
    
    	for(i=0;i<N;i=i++){
    		(Out+i)->re=creal(*(In+i));
    		(Out+i)->im=cimag(*(In+i));
    	};
    };
    
    /******************ComplexMultiplication**************************/
    Complex Cexpmult(Complex a, Complex b){
    	Complex z;
    	z.re=a.re*b.re-a.im*b.im;
    	z.im=a.re*b.im+a.im*b.re;
    	return z;
    
    }
    
    /****************Discrete FourierTransformation**************************/
    void FTrans (Complex* In, Complex* Out, int N, int dir){
    	int j=0;
    	int k=0;
    	Complex Phase[N];
    	for(j=0;j<N;j=j++){
    		for(k=0;k<N;k=k++){
    
    			Phase[k].re=cos(atan(dir * 2 * M_PI * k * j/N)); /* Phasefactor Realpart*/
    			Phase[k].im=sin(atan(dir * 2 * M_PI * k * j/N));/* Phasefactor Imaginarypart*/
    
    			printf("Phasefactor:Re:%f Im: %f\n", Phase[k].re,Phase[k].im);
    
    			Complex f_k= Cexpmult(Phase[k] ,*(In+k));
    
    			(Out+j)->re=(Out+j)->re +1/sqrt(N) * f_k.re;
    			(Out+j)->im=(Out+j)->im +1/sqrt(N) * f_k.im;
    
    		};
    
    	};
    };
    
    /*********Function to Print out the results**********/
    void PrintOut(Complex* In, int N){
    	int i;
    	for(i=0;i<N;i=i++){
    		printf("Point %i: %f + I * %f\n", i+1, (In+i)->re,(In+i)->im);
    	};
    
    };
    
    int main (){
    	int 	N=10; /* Number of data point*/
    	complex double* Data = calloc(N , sizeof(complex double));
    	Complex* Split = calloc(N , sizeof(Complex));
    	Complex* FDT = calloc(N , sizeof(Complex));
    
    	int i;
    
    	/********Generation of testing data*************/
    	for(i=0;i<=N;i=i++){
    		*(Data+i)=cos(i * M_PI/2);
    		printf("Point %i: Re=%f\n",i,creal(*(Data+i)));
    	};
    /**********************Splitting Data*********************************/
    
    		ComplexSplit(Data , Split, N);
    		for(i=0;i<=N;i=i++){
    		printf("Im=%f,Re=%f\n",Split[i].im,Split[i].re);
    
    	};
    
    /***************Transformation and Print***************/		
    	FTrans(Split ,FDT,N,1);
    	PrintOut(FDT, N);
    	return 0;
    };
    


  • Kannst Du netterweise erklären, was genau Dein Problem ist? Du schreibst, daß Du Probleme mit den komplexen Zahlen hast, was meinst Du damit? Was geschieht bei der Transformation eines Testvektors (z.B. alle Werte 0, oder alle Werte 1)?



  • Was genau ist dein Problem?
    "Ich scheitere" ist keine Fehlerbeschreibung und genau so bescheuert wie "es klappt nicht"

    Welche Meldungen (Warnungen, Fehler) gibt es vom Compiler?
    Welche vom Linker?

    Warum hast du ein #include <complex.h> und definierst danach einen eigenen Complex-Datentyp?



  • Ich definiere einen eigenen Complex typ, weil es mir schien, dass man damit einfacher Zeigerarithmetik machen kann, dies gelang mir mit complex double nicht.

    Ich dachte aber, dass ich complex double am ende besser von einem File einlesen kann in dem die Daten in z.B. polarform stehen.

    Wenn ich einen Testvektor (0,0,0....) nehme, ist die FDT auch 0.

    Aber wenn ich z.b. 10 Punkte von cos(nπ2)cos(\frac{n\pi}{2})
    und hin und zurück transformiere bekomme ich das: (erst eingangs daten, dann hin, dann zurück):

    Im=0.000000,Re=1.000000
    Im=0.000000,Re=0.000000
    Im=0.000000,Re=-1.000000
    Im=0.000000,Re=-0.000000
    Im=0.000000,Re=1.000000
    Im=0.000000,Re=0.000000
    Im=0.000000,Re=-1.000000
    Im=0.000000,Re=-0.000000
    Im=0.000000,Re=1.000000
    Im=0.000000,Re=0.000000
    Im=0.000000,Re=0.000000
    Point 1: 0.316228 + I * 0.000000
    Point 2: 0.216853 + I * 0.050875
    Point 3: 0.250746 + I * 0.017518
    Point 4: 0.269799 + I * 0.008361
    Point 5: 0.280611 + I * 0.004828
    Point 6: 0.287429 + I * 0.003128
    Point 7: 0.292088 + I * 0.002187
    Point 8: 0.295463 + I * 0.001613
    Point 9: 0.298017 + I * 0.001239
    Point 10: 0.300015 + I * 0.000980
    Point 1: 0.887730 + I * 0.028691
    Point 2: 0.406550 + I * -0.689602
    Point 3: 0.290529 + I * -0.745866
    Point 4: 0.242162 + I * -0.763635
    Point 5: 0.215729 + I * -0.771609
    Point 6: 0.199157 + I * -0.775915
    Point 7: 0.187833 + I * -0.778534
    Point 8: 0.179619 + I * -0.780266
    Point 9: 0.173396 + I * -0.781482
    Point 10: 0.168521 + I * -0.782376
    

    beim compilieren bekomme ich das:

    **** Build of configuration Debug for project FDT ****
    
    make all 
    Building file: ../main.c
    Invoking: Cross GCC Compiler
    gcc -O0 -g3 -Wall -c -fmessage-length=0 -MMD -MP -MF"main.d" -MT"main.d" -o "main.o" "../main.c"
    ../main.c: In Funktion »ComplexSplit«:
    ../main.c:21:15: Warnung: Operation auf »i« könnte undefiniert sein [-Wsequence-point]
    ../main.c: In Funktion »FTrans«:
    ../main.c:42:15: Warnung: Operation auf »j« könnte undefiniert sein [-Wsequence-point]
    ../main.c:43:16: Warnung: Operation auf »k« könnte undefiniert sein [-Wsequence-point]
    ../main.c: In Funktion »PrintOut«:
    ../main.c:66:15: Warnung: Operation auf »i« könnte undefiniert sein [-Wsequence-point]
    ../main.c: In Funktion »main«:
    ../main.c:75:2: Warnung: Implizite Deklaration der Funktion »calloc« [-Wimplicit-function-declaration]
    Finished building: ../main.c
    
    Building target: FDT
    Invoking: Cross GCC Linker
    gcc  -o "FDT"  ./main.o   -lm
    ../main.c:75:25: Warnung: Unverträgliche implizite Deklaration der eingebauten Funktion »calloc« [standardmäßig aktiviert]
    ../main.c:84:16: Warnung: Operation auf »i« könnte undefiniert sein [-Wsequence-point]
    ../main.c:91:17: Warnung: Operation auf »i« könnte undefiniert sein [-Wsequence-point]
    Finished building target: FDT
    
    **** Build Finished ****
    


  • ich hatte vergessen bei der rücktransformation komplex zu konjugieren.
    daran liegt es jedoch nicht.



  • Erst mal die Warnungen:

    Entweder 
    i = i+1
    oder
    i+=1
    oder
    i++ oder ++i
    

    Aber nicht i=i++ auf gar keinen Fall.

    Das mit calloc bekommst du mit einem

    #include <stdlib.h>
    

    weg.



  • Danke, jetzt sind alles Fehlermeldungen weg. Die FDT funktioniert noch nicht. Aber immerhin sind die groben Schnitzer raus.



  • Ist das Out-Feld bei der Rücktrannsformation auch leer?

    Mach mal zwischen Zeile 42 und 43 (vor dem for(k=...) )

    (Out+j)->re= 0;
    (Out+j)->im= 0;
    


  • Hier ist nochmal der Code, es änderte nichts.

    /*
     * main.c
     *
     *  Created on: 05.10.2012
     *      Author: herrk
     */
    
    #include<stdio.h>
    #include<math.h>
    #include<complex.h>
    #include<stdlib.h>
    
    /*****************Create Structur for Complex Numbers**************/
    typedef struct {double re; double im;} Complex;
    
    /***************Data in complex Double form gets stored in structure*************/
    void ComplexSplit (complex double* In, Complex* Out, int N){
    	int i;
    
    	for(i=0;i<N;i++){
    		(Out+i)->re=creal(*(In+i));
    		(Out+i)->im=cimag(*(In+i));
    	};
    };
    
    /******************ComplexMultiplication**************************/
    Complex Cexpmult(Complex a, Complex b){
    	Complex z;
    	z.re=a.re*b.re-a.im*b.im;
    	z.im=a.re*b.im+a.im*b.re;
    	return z;
    
    }
    
    /****************Discrete FourierTransformation**************************/
    void FTrans (Complex* In, Complex* Out, int N, int dir){
    	int j=0;
    	int k=0;
    	Complex Phase[N];
    	for(j=0;j<N;j++){
    		(Out+j)->re= 0;
    		(Out+j)->im= 0;
    		for(k=0;k<N;k++){
    
    			Phase[k].re=cos(atan(dir * 2 * M_PI * k * j/N)); /* Phasefactor Realpart*/
    			Phase[k].im=sin(atan(dir * 2 * M_PI * k * j/N));/* Phasefactor Imaginarypart*/
    
    			/*printf("Phasefactor:Re:%f Im: %f\n", Phase[k].re,Phase[k].im);*/
    
    			(In+k)->im=dir*(In+k)->im;
    
    			Complex f_k= Cexpmult(Phase[k] ,*(In+k));
    
    			(Out+j)->re=(Out+j)->re +1/sqrt(N) * f_k.re;
    			(Out+j)->im=(Out+j)->im +1/sqrt(N) * f_k.im;
    
    		};
    
    	};
    };
    
    /*********Function to Print out the results**********/
    void PrintOut(Complex* In, int N){
    	int i;
    	for(i=0;i<N;i++){
    		printf("Point %i: %f + I * %f\n", i+1, (In+i)->re,(In+i)->im);
    	};
    
    };
    
    int main (){
    	int 	N=10; /* Number of data point*/
    	complex double* Data = calloc(N , sizeof(complex double));
    	Complex* Split = calloc(N , sizeof(Complex));
    	Complex* FDT = calloc(N , sizeof(Complex));
    	Complex* Back = calloc(N , sizeof(Complex));
    
    	int i;
    
    	/********Generation of testing data*************/
    	for(i=0;i<N;i++){
    		*(Data+i)= cos(M_PI/2 * i);
    	;
    	};
    /**********************Splitting Data*********************************/
    
    		ComplexSplit(Data , Split, N);
    		for(i=0;i<N;i++){
    		printf("Im=%f,Re=%f\n",Split[i].im,Split[i].re);
    
    	};
    
    /***************Transformation and Print***************/
    	FTrans(Split ,FDT,N,1);
    	PrintOut(FDT, N);
    	FTrans(FDT, Back, N, -1);
    	PrintOut(Back, N);
    	return 0;
    };
    

    Calloc schreibt doch überall 0 , oder?

    Ausgabe war:

    Im=0.000000,Re=1.000000
    Im=0.000000,Re=0.000000
    Im=0.000000,Re=-1.000000
    Im=0.000000,Re=-0.000000
    Im=0.000000,Re=1.000000
    Im=0.000000,Re=0.000000
    Im=0.000000,Re=-1.000000
    Im=0.000000,Re=-0.000000
    Im=0.000000,Re=1.000000
    Im=0.000000,Re=0.000000
    Point 1: 0.316228 + I * 0.000000
    Point 2: 0.216853 + I * 0.050875
    Point 3: 0.250746 + I * 0.017518
    Point 4: 0.269799 + I * 0.008361
    Point 5: 0.280611 + I * 0.004828
    Point 6: 0.287429 + I * 0.003128
    Point 7: 0.292088 + I * 0.002187
    Point 8: 0.295463 + I * 0.001613
    Point 9: 0.298017 + I * 0.001239
    Point 10: 0.300015 + I * 0.000980
    Point 1: 0.887730 + I * -0.028691
    Point 2: 0.406550 + I * -0.689602
    Point 3: 0.241215 + I * -0.772689
    Point 4: 0.242162 + I * -0.763635
    Point 5: 0.160912 + I * -0.787036
    Point 6: 0.199157 + I * -0.775915
    Point 7: 0.131655 + I * -0.789157
    Point 8: 0.179619 + I * -0.780266
    Point 9: 0.116705 + I * -0.789545
    Point 10: 0.168521 + I * -0.782376
    


  • die Formeln waren falsch ich habe den Kopf verloren. Ich muss natürlich nicht den Arkustan von den Phasen ziehen....

    Danke , es war trotzdem eine Hilfe


  • Mod

    Gibt es eigentlich einen tieferen Grund, eine Fouriertrafo selber zu implementieren? Es gibt so gute Bibliotheken dafür…



  • SeppJ schrieb:

    Gibt es eigentlich einen tieferen Grund, eine Fouriertrafo selber zu implementieren? Es gibt so gute Bibliotheken dafür…

    Ja, gibt es.
    Hast du denn etwa selbst nie mehr oder weniger Ähnliches mit mehr oder weniger ähnlich viel Aufwand probiert/gemacht/verwirklicht? 😕


  • Mod

    xyz offgeloggt schrieb:

    SeppJ schrieb:

    Gibt es eigentlich einen tieferen Grund, eine Fouriertrafo selber zu implementieren? Es gibt so gute Bibliotheken dafür…

    Ja, gibt es.
    Hast du denn etwa selbst nie mehr oder weniger Ähnliches mit mehr oder weniger ähnlich viel Aufwand probiert/gemacht/verwirklicht? 😕

    Ich verstehe nicht, was du mir sagen möchtest.



  • Man muss nicht alles wirklich immer alles mit externen Libs lösen, das ist schlichtweg Unsinn. Wenn es sich vermeiden lässt empfiehlt, es sich immer seine Software von externen Bibliotheken zu entkoppeln. Abhängigkeiten nur aufbauen wenn es wirklich unbedingt nötig ist. So sieht es jedenfalls bei kommerziellen Produktionen aus, im privaten Lernumfeld implementiert man eh wirklich fast alles erst einmal selbst, so bekommt man ein tieferes Verständnis für ganze Technologie.


  • Mod

    Klar braucht man nicht für alles eine externe Bibliothek, aber für Fouriertrafos? Die werden so oft gebraucht, sind aber so schwer (gut) zu implementieren. Die Bibliotheken sind zudem in der Regel extrem klein. Ich wüsste wirklich nicht, was dafür sprechen würde, das selber zu implementieren.



  • SeppJ schrieb:

    Klar braucht man nicht für alles eine externe Bibliothek, aber für Fouriertrafos? Die werden so oft gebraucht, sind aber so schwer (gut) zu implementieren. Die Bibliotheken sind zudem in der Regel extrem klein. Ich wüsste wirklich nicht, was dafür sprechen würde, das selber zu implementieren.

    1. Zugriff auf den Quellcode um Fehler zu beheben
    2. Zugriff auf den Quellcode um Geschwindigkeitsoptimierungen für bestimmte Prozessoren durchzuführen
    3. Zugriff auf den Quellcode zur Anpassung des Interfaces (Zur Vermeidung von nötigen Adapter Pattern)
    4. Lizenz-rechtliche Sicherheit
    5. Tieferes intuitives Verständnis
    6. Selbstbestätigung das man "es kann"
    7. Hobby


  • Mod

    MisterX schrieb:

    1. Zugriff auf den Quellcode um Fehler zu beheben

    Open Source? Außerdem: Was hat wohl eher einen Fehler? Die vielbenutzte Bibliothek oder die Eigenimplementierung?

    2. Zugriff auf den Quellcode um Geschwindigkeitsoptimierungen für bestimmte Prozessoren durchzuführen

    Genau wie Punkt 1. Wenn du es schaffst, eine gängige FT-Bibliothek in Eigenarbeit zu verbessern, dann bewirb dich bei denen als Chefentwickler.

    3. Zugriff auf den Quellcode zur Anpassung des Interfaces (Zur Vermeidung von nötigen Adapter Pattern)

    Siehe Punkt 1.

    4. Lizenz-rechtliche Sicherheit

    Dann nimm eben eine mit passender Lizenz.Die gängigen Bibliotheken haben recht laxe Lizenzen. Die Hersteller wissen, dass schnelle FT sehr schwer zu implementieren ist, aber sehr oft gebraucht wird. Wenn sie es frei machen, pushen sie ihre Plattform.

    5. Tieferes intuitives Verständnis

    Ok.

    6. Selbstbestätigung das man "es kann"

    Ok. Guck dir den Code des TE an. Er ist weit von einer guten Implementierung entfernt. Daher sollte man so etwas schon mal vorschlagen.

    7. Hobby

    Ok.


Log in to reply