Implementierung von Sinustransformation



  • Hallo zusammen,

    ich habe Schwierigkeiten die Sinustransformation als Programm zu implementieren. (Dies *war* eine Programmieraufgabe für die Uni, deren Deadline heute morgen gewesen ist).

    Leider funktioniert meine Implementierung nicht so richtig. Wenn ich zu gegeben y-Werte die alphas berechne und sie danach wieder zurücktransformiere, komme ich nicht auf die Ursprungsfunktion.

    Hat schonmal jemand die Sinustransformation implementiert und kann mir sagen, wo etwas falsch ist? (Ich weiß, dass der ein oder andere free fehlt und man noch ein paar Sachen abfangen müsste (fopen zB)).

    #include <stdio.h>
    #include <stdlib.h>
    #include <time.h>
    #include <math.h>
    
    double* fft(double*x,int k)
    {
    	int m=pow(2,k)/2;
    	int count=0;
    	double* diag;
    	diag=malloc(m*sizeof(double));
    	double* u_1;
    	u_1=malloc(m*sizeof(double));
    	double* u_2;
    	u_2=malloc(m*sizeof(double));
    	double* v;
    	v=malloc(m*sizeof(double));
    	double* w;	
    	w=malloc(m*sizeof(double));
    	double* y;
    	y=malloc(2*m*sizeof(double));
    	if((diag==NULL)||(u_1==NULL)||(u_2==NULL)||(v==NULL)||(w==NULL))
    	{
    		printf("\n\n\n OOM \n\n\n");
    		return NULL;
    	}
    	if(k>1)
    	{
    		double q = 3.141592653589793238462643383/m;
    		for(count=0;count<m;count++)
    		{
    			diag[count]=sin(q*count);
    		}
    		for(count=0;count<m;count++)
    		{
    			u_1[count]=x[count]+x[count+m];
    			u_2[count]=(x[count]-x[count+m])*diag[count];;
    		}
    		w=fft(u_1,k-1);
    		v=fft(u_2,k-1);
    		for(count=0;count<m;count++)
    		{
    			y[2*count]=v[count];
    			y[2*count+1]=w[count];
    		}
    	}
    	else
    	{
    		for(count=0;count<2*m;count++)
    		{
    			y[count]=x[count];
    		}
    	}
    	free(diag);
    	free(u_1);
    	free(u_2);
    	free(v);
    	free(w);
    	return y;
    }
    
    void data(double**x,double**y,double** alpha,int k,double start,double end)
    {
    
    	srand(time(NULL));
    	int count=0;
    	double pi=3.141592653589793238462643383;
    	double dif=end-start;
    	double anz=pow(2,k);
    	double h=dif/anz;	
    	FILE* before;
    	FILE* origin;
    	before=fopen("before","w");
    	origin=fopen("origin","w");
    	(*x)=malloc(anz*sizeof(double));
    	(*y)=malloc(anz*sizeof(double));
    	(*alpha)=malloc(anz*sizeof(double));	
    	for(count=0;count<anz;count++)
    	{
    		double zufall=rand()%100;
    		zufall=zufall/100;
    		(*x)[count]=start+count*h;
    		(*y)[count]=0.7*sin(((*x)[count])*2*pi)+0.94*sin(((*x)[count])*16*pi)+2*(zufall-0.5);
    		fprintf(before,"%f %f\n",(*x)[count],(*y)[count]);
    		fprintf(origin,"%f %f\n",(*x)[count],0.7*sin(((*x)[count])*2*pi)+0.94*sin(((*x)[count])*16*pi));
    	}
    	fclose(before);
    	fclose(origin);
    }
    
    void filter(double*alpha,double gamma, int k)
    {
    	int count=0;
    	double max=0;
    	for(count=0;count<(pow(2,k));count++)
    	{
    		if(max<alpha[count])
    		{
    			max=alpha[count];
    		}
    	}
    	for(count=0;count<(pow(2,k));count++)
    	{
    		if(alpha[count]<(gamma*max))
    		{
    			alpha[count]=0;
    		}
    	}
    }
    
    int main()
    {
    	double* x;
    	double* y;
    	double* alpha;
    	double start=0;
    	double end=3.141592653589793238462643383;
    	double gamma=0.1;
    	int k=11;
    	int count=0;
    	data(&x,&y,&alpha,k,start,end);	
    	alpha=fft(y,k);
    	for(count=0;count<pow(2,k);count++)
    	{
    		alpha[count]=-1*alpha[count]/pow(2,k);
    	}
    	FILE* fourier;
    	fourier=fopen("fourier","w");
    	for(count=0;count<pow(2,k);count++)
    	{
    		fprintf(fourier,"%f %f\n",x[count],alpha[count]);
    	}
    	//filter(alpha,0,k);
    	y=fft(alpha,k);
    	FILE* after;
    	after=fopen("after","w");
    	for(count=0;count<pow(2,k);count++)
    	{
    		fprintf(after,"%f %f\n",x[count],y[count]);
    	}
    	free(x);
    	free(y);
    	free(alpha);
    }
    

    Gruß
    Isolde123



  • isolde123 schrieb:

    u_1[count]=x[count]+x[count+m]; /* hier fehlt doch wohl was */
    
    for(count;count<anz;count++) /* hier auch */
    


  • Hallo Wutz,

    erstmal vielen Danke dass du dir Zeit genommen hast, den Code zu durchsuchen und mir zu antworten!

    Die 2. Schleife habe ich verbessert. Du meintest doch wohl, dass es

    for(count=0;count<anz;count++)
    

    heißen soll oder? Das hat aber zum Glück keinen Schaden angerichtet, da ich die Variable anständig initialisiert habe.

    Zu dem anderem Punkt:
    Es steht so in unserem Skript drin, nur leider scheint das falsch zu sein. Wir haben nämlich:

    F8=[F_4F_4RR]F_8 = \begin{bmatrix} F\_4 & F\_4\\ R & -R \end{bmatrix}

    ,
    was ja gar nicht stimmt. Das müsste nämlich nicht F4F_4 sein, sondern die Matrix die aus dieser entsteht, wenn man alle Einträge quadriert.

    Wir haben den Algorithmus nicht selber hergeleitet, sondern vom Prof. in Pseudocode vorgesetzt bekommen, daher weiß ich die genaue Herleitung nicht.

    Mein Schuß ins Blaue wäre jetzt, dass in der Zeile auch die Quadrate fehlen, kann das sein?



  • Ob da ein Quadrat fehlt oder nicht, kann ich dir leider nicht sagen, sowas wäre etwas für das Mathematik-Subforum hier.
    Aber mind. sollte

    for(count=0;count<m;count++)
            {
                u_1[count]=x[count]+x[count+m];
                u_2[count]=(x[count]-x[count+m])*diag[count];;
            }
    

    doch wohl eher

    for(count=0;count<m;count++)
            {
                u_1[count]=(x[count]+x[count+m])*diag[count];
                u_2[count]=(x[count]-x[count+m])*diag[count];
            }
    

    heißen.
    Außerdem solltest du für Pi eine symbolische Konstante verwenden, z.B.

    #ifndef M_PI
    #define M_PI 3.141592653
    #endif
    

    Weiterhin ist deine Speicherbehandlung, wie du selber schon bemerkt hast, stark verbesserungsbedürftig, du brauchst kein extra Array zurückgeben, nehme besser einfach das übergebene Array zur Abspeicherung der Ergebniswerte, also in etwa (ungetestet):

    void fft(double*x,int k)
    {
      if(k>1)
      {
        int count, m=pow(2,k)/2;
        double* diag=malloc(m*sizeof*diag);
        double* v =malloc(m*sizeof*v);
        double* w =malloc(m*sizeof*w);
        double q = M_PI/m;
        for(count=0;count<m;count++)
        {
          diag[count]=sin(q*count);
        }
        for(count=0;count<m;count++)
        {
          v[count]=(x[count]+x[count+m])*diag[count];
          w[count]=(x[count]-x[count+m])*diag[count];
        }
        fft(v,k-1); /* nach Ausführung der Funktion steht das Ergebnis gleich in v */
        fft(w,k-1); /* nach Ausführung der Funktion steht das Ergebnis gleich in w */
        for(count=0;count<m;count++)
        {
          x[2*count]=v[count];
          x[2*count+1]=w[count];
        }
        free(w);free(v);free(diag);
      }
    }
    

Log in to reply