Bug in FFT-Routine

[GraveCH]

Hallo Leute
Ich habe mich nun seit einiger Zeit mit dem FFT beschäftigt. Um mit meinem rechner den fft benützen zu können habe ich die funktionen von http://teacher.schule.at/krammer/mist.htm heruntergeladen. danach habe ich den code nach c++ portiert, damit ich das selbe auch auf meinem PC machen kann.
Nun ist mein problem, dass die Resultate von denen meines rechners abweichen. es muss sich beim portieren ein fehler eingeschleicht haben, den ich nicht finde.

hier ist der code in c++ (sorry wegen dem etwas chaotischen stil):

void fft(unsigned const int n, double *inreal, double *inimag);
void fftw(unsigned int anz, double *wreal, double *wimag);
void fftrev(unsigned const int n, double *inreal, double *inimag);

void fft(unsigned const int n, double *inreal, double *inimag)
{
	fftrev(n, inreal, inimag);

	int a = 1;
	int b = n;
	int d = 0, c = 0, i = 0, k = 0;
	double t;

	double *wreal = new double[n/2];
	double *wimag = new double[n/2];
	fftw(n, wreal, wimag);

	while (a<n)
	{
		c = a+a;
		b = b/2;
		d = 1;
		for (int m = 1; m<=a; m++)
		{
			i = m;
			while (i<=n)
			{
				k = i+a;
				t = wreal[d-1]*inreal[k-1];
				inreal[k-1] = inreal[i-1]-t;
				inreal[i-1] = inreal[i-1]+t;

				t = wimag[d-1]*inimag[k-1];
				inimag[k-1] = inimag[i-1]-t;
				inimag[i-1] = inimag[i-1]+t;

				i = i+c;
			}
			d = d+b;
		}
		a = c;
	}
	for (i=0; i<n;i++)
	{
		inreal[i] = inreal[i]/(n);
		inimag[i] = inimag[i]/(n);
	}
	delete wreal;
	delete wimag;
}

void fftrev(unsigned const int n, double *inreal, double *inimag)
{
	unsigned int i;
	unsigned int j=1;
	unsigned int m=0;

	double temp;

	for (i=1; i<n; i++)
	{
		if(j>i)
		{
			temp = inreal[j-1];
			inreal[j-1] = inreal[i-1];
			inreal[i-1] = temp;

			temp = inimag[j-1];
			inimag[j-1] = inimag[i-1];
			inimag[i-1] = temp;
		}
		m = n/2;
		while (m>=2 && j>m)
		{
			j = j-m;
			m = m/2;
		}
		j = j+m;
	}
}

void fftw(unsigned int anz, double *wreal, double *wimag)
{
	/*anz ist die Anzahl Elemente. 
	Es gibt doch bloss anz/2 Elemente bei wreal und wimag. 
	Für den IFFT muss wimag mit -1 multipliziert werden.*/

	unsigned int i;
	double c;
	const double PI = 3.1415926535;

	anz = anz/2.0;
	c = PI/anz;

	for (i=0;i<anz; i++)
	{
		wreal[i] = cos(c*i);
		wimag[i] = -sin(c*i);
	}
}

ich habe also 3 funktionen. fftw verlangt die anzahl elemente und 2 zeiger auf ein array mit double werten. diese werte werden dann geschrieben und stellen den realteil und den imaginärteil der drehwerte dar.

fftrev ist zuständig für das bitreversing.

fft ist die hauptfunktion.

ich habe bereits etwas herumexperimentiert und die werte aus fftrev und fftw sind fehlerfrei.
wenn ich in der fft-funktion das:

k = i+a;
				t = wreal[d-1]*inreal[k-1];
				inreal[k-1] = inreal[i-1]-t;
				inreal[i-1] = inreal[i-1]+t;

				t = wimag[d-1]*inimag[k-1];
				inimag[k-1] = inimag[i-1]-t;
				inimag[i-1] = inimag[i-1]+t;

zu:

k = i+a;
				t = 1.0*inreal[k-1];
				inreal[k-1] = inreal[i-1]-t;
				inreal[i-1] = inreal[i-1]+t;

				t = 1.0*inimag[k-1];
				inimag[k-1] = inimag[i-1]-t;
				inimag[i-1] = inimag[i-1]+t;

ändere sind die resultate die selben wie auf meinem rechner(der code auf dem rechner wurde auch dementsprechend verändert).

der fehler muss also in dieser multiplikation sein. inreal[] und inimag[] sind die einganswerte und diese sind logischerweise nicht für den fehler verantwortlich.

ich vermute, dass vielleicht rundungsfehler entstehen und diese ausreichen um meine resultate dementsprechend zu verfälschen.

seht ihr gerade woran es liegt?
danke vielmals für die hilfe!
grs Grave

[GraveCH]

das problem hat sich erledigt.
da es beim fft eine KOMPLEXE multiplikation benötigt kann ich nicht sagen, dass x*y = real(x)*real(y) + imag(x)*imag(y)*i --> x und y sind komplexe zahlen

die komplexe multiplikation geht folgendermassen:
x*y = (real(x)*real(y) - (imag(x)*imag(y))) + (real(x)*imag(y) + imag(x)*real(y))*i

ein wirklich dummer fehler, der mir viele nerven gekostet hat....
danke trotzdem. falls jemand gerne nähere informationen haben möchte, wie nun meine funktionen funktionieren ( ) darf sich gerne bei mir melden.
grs grave