FFT-Problem

Hi!
Mein FFT-Code gibt mir zwar dieselben Fourierkoeffizienten wie Maple, aber bei der Funktionsauswertung (eval) kommt auch bei höherer Stützstellenanzahl nicht die Gerade heraus

was mache ich falsch?

#include <iostream>
#include <vector>
#include <complex>
using namespace std;

typedef vector< complex<double> > fourier_coeffs;

inline unsigned log2(unsigned value)
{
	unsigned result = 0;
	do {
		++result;
		value >>= 1;
	} while (value != 0);
	return result-1;
}

const double Pi = 3.14159265358979323846264338327950288;

unsigned int bit_reverse(unsigned int x, unsigned N) {
	int n = 0;
	int mask = 0x1;
	for (int i=0; i < N; i++) {
		n <<= 1;
		n |= (x & 1);
		x >>= 1;
	}
	return n;
}

fourier_coeffs FFT(const vector<double>& f)
{
	size_t N = f.size();
	size_t n = log2(N);
	fourier_coeffs b;
	b.resize(N);
	for (size_t i = 0; i < N; ++i)
		b[bit_reverse(i,n)] = complex<double>(f[i]);

	for (size_t m = 1; m <= n; ++m)
		for (size_t j = 0; j < (1 << (m-1)); ++j) {
			complex<double> e = exp(complex<double>(0, -Pi*j/(1<<(m-1))));
			//complex<double> e = complex<double>(cos(Pi*j/(1<<(m-1))), -sin(Pi*j/(1<<(m-1))));
			for (size_t r = 0; r < N; r += (1 << m)) {
				complex<double> u = b[r + j],
					v = b[r + j + (1 << (m-1))] * e;
				b[r + j] = u + v;
				b[r + j + (1 << (m-1))] = u - v;
			}
		}

	for (size_t i = 0; i < N; ++i)
		b[i] /= N;

	return b;
}

complex<double> eval(const fourier_coeffs& b, double x)
{
	complex<double> f = 0;
	for (size_t i = 0; i < b.size(); ++i)
		f += b[i]*exp(complex<double>(0, i*x));
	return f;
}

int main()
{
	// f periodisch [-pi, pi)
	const unsigned n = 2;
	const size_t N = 1 << n;
	vector<double> f,x;
	f.resize(N);
	x.resize(N);
	for (size_t i = 0; i < N; ++i)
		f[i] = x[i] = -Pi + 2*Pi*i/N;
	fourier_coeffs b = FFT(f);
	for (size_t i = 0; i < N; ++i)
		cout << b[i] << "\n";

	cout << "\n";
	for (size_t i = 0; i < N; ++i)
		cout << eval(b, -Pi + 2*Pi*i/N) << "\n";

	cout << "\n";
	for (size_t i = 0; i < N; ++i)
		cout << "(" << b[i].real() << "+I*(" << b[i].imag() << "))*exp(I*" << i << "*x) + ";
	return 0;
}

Jester

fft schrieb:

Mein FFT-Code gibt mir zwar dieselben Fourierkoeffizienten wie Maple, aber bei der Funktionsauswertung (eval) kommt auch bei höherer Stützstellenanzahl nicht die Gerade heraus

Die Gerade kommt nicht heraus? Äh, moment. Welche Gerade denn? Kannst Du vielleicht kurz erklären, was Du überhaupt genau machen möchtest?

Die Gerade:

for (size_t i = 0; i < N; ++i)
        f[i] = x[i] = -Pi + 2*Pi*i/N;

FFT berechnet die Fourier-Koeffizienten für die Stützstellen -Pi + 2*Pi*i/N mit Werten f[i]. Eval berechnet dann das Fourierpolynom.

Turing

Hast Du Dir schonmal ansatzweise überlegt, wie das Spektrum Deiner Funktion wohl aussehen mag?

Ich versuchs mal:

Deine Gerade ist ein Dreieck, das mit 'nem Rechteck auf halber Distanz gefaltet wurde. Im Frequenzbereich wäre das eine Multiplikation von drei Si-Fkt. Die Si-Fkt. ist unendlich ausgedehnt und das Abtasttheorem fordert ja bekanntlich Bandbegrenzung. Ergo, kannst Du die fehlerfreie Rücktrafo vergessen.

Die Ergebnisse sind aber total falsch, an den originalen Stützstellen müssen ja die richtigen Ergebnisse rauskommen!

Turing

wieso?

Weil für die Koeffizienten gilt b_k = <f, exp(nix)> (Hermitesches Skalarprodukt, x = x_i, i = 0...n-1), P(x_k) = sum(b_k*exp(nix_k)) = <b, exp(-nix)> = f_k ist.

(also weil exp(nix) Orthogonalbasis)

Turing

Tut mir Leid, was das Hermitesche Skalarprodukt angeht, kann ich Dir leider nicht weiter helfen.

Dafür kann ich Dir aber aus systemtheoretischer Sicht begründen, warum die Fehler auftreten:

Du hast eine diskrete Zeitfunktion, diese hat immer ein periodisches Spektrum. Deine Zeitfunktion ist nicht periodisch. Deswegen ist das Spektrum kontinuierlich.
(Kannst Du z.B. in Digital Signal Processing von Proakis nachlesen oder mir einfach glauben)

Folglich musst Du das gesamte kontinuierliche Spektrum rücktransformieren, um die ursprünglichen diskreten Werte zu erhalten, alles andere führt zu Fehlern.

Außerdem ist Deine Zeitfunktion nicht bandbegrenzt und es kommt zu den sogenannten spektralen Überlappungen, welche wohl eigentlich am Fehler Schuld sind.

Turing schrieb:

Deine Zeitfunktion ist nicht periodisch.

Achso, sorry für die Ungenauigkeit. Meine lineare Funktion wird auf [-pi,pi] so definiert und dann periodisch fortgesetzt.