fourier reihe II elit sehr! bitte helfen!!!

[Heinrich der IT-genie 78]

Ich habe ebenso versucht ein programm für die berechnung der fourier transformation, bloß kein erfolg. ich muss dass bis montag fertigen. Hilft mir bitte

#include <studio.h>

#include <math.h>

::::

double fct( int chf, double x)

// Funktionen zur Auswahl

{

x=(-1.0)+x/M_PI; // symmetrisiert [-1,1] auf [0, 2 pi]

if (chf==0)

{

if (x>=0.1)

return(sqrt(x));

else

return(sqrt(-x));

}

if (chf==1) // sollte exakt reproduziert werden,

{ // wenn n>=4

return(cos(4.0*M_PI*(x+1.0)));

}

if (chf==2) // Gaussglocke

{

return(exp(-10.0*x*x));

}

if (chf==3) // Runge

{

return(2.0/(1.0+15.0*x*x));

}

if (chf==4) // Betrag

{

return(fabs(x));

}

if (chf==5) // Sprung

{

if (x>=0.0)

return(1.0);

else

return(-1.0);

}

return 0.0;

}

int bitinverse(int k, int n)

// bildet Bitinversion von k < n=2 hoch p

{ // unter Vermeidung von Shiftoperatoren.

int c=0; // Das kann effizienter gemacht werden!

if (n<=1) // Effizientere Versionen sind leider

return(0); // noch schwerer verstaendlich....

// c sammelt das Ergebnis auf,

// n laeuft absteigend durch die Zweierpotenzen

while(n>1) // Logarithmische Ausfuehrungszeit!

{

n/=2; // fuer hoechstes Bit

if (k%2) // hole unteres Bit von k

c+=n ; // baue es vorn in c ein

k/=2; // reduziere k am unteren Ende

}

return©;

}

void split(int n, int inv,

// Nichtrekursive Verteilungsroputine.

Vector real, Vector imag)

// inv = 1/-1: FFT/inverse FFT.

{ // Die Vektoren werden ueberschrieben!

if (RSNMoutput)

{

cout<<"Aufruf split fuer n = "<<n<<"\n";

for (int i=0; i<n; i++)

cout.form("%3d : %8.4f + i * %8.4f \n"

,i,real,imag);

}

int m, mstart, mend, mpj;

double c, s, rj, ij, c1, s1, ch, signum=(double)inv;

// Ab hier: Verteilungsschritte:

for (m=n/2; m>=1; m/=2) // es gibt log(n) Schritte...

{ // Ab hier kommt die nichtrekursive Form

// des naiven Programms.

// Es werden Teilarrays der Laenge 2*m bearbeitet,

// und zwar in zwei Segmente der Laenge m aufgeteilt,

// wobei die ueblichen FFT-Formeln angewendet werden:

// rho = 0 und 1, siehe Buch, Formel (12.2.27).

// Dabei laeuft man ueber die arrays weg (Datenlokalitaet!)

// und arbeitet auf Indizes in [mstart,mstart+2m).

for (mstart=0; mstart<n-m; mstart+=2*m) // Laufschleife...

{

mend=mstart+m;

// die Schleife unten laeuft nur ueber m Elemente,

if (RSNMoutput) // bearbeitet aber 2*m Elemente!

{

cout<<"Aufruf split fuer m = "<<m<<

" zwischen "<<mstart<<

" und "<<mend+m-1<<"\n";

for (int i=mstart; i<mend+m; i++)

cout.form("%3d : %8.4f + i * %8.4f \n",

i,real,imag);

}

// Die Drehfaktoren werden hier durch trigonometrische

// Formeln iterativ berechnet, umd sin/cos-Aufrufe zu sparen.

// Ob das effizienter ist, ist fraglich, weil das Timing

// eines sin/cos-Aufrufes durchaus mit dem einer Multiplikation

// vergleichbar sein kann (bei intelligenten Coprozessoren).

c1=cos(signum*M_PI/((double)m));

// Das ist \zeta_{2m}^{signum},

s1=sin(signum*M_PI/((double)m));

// Real- und Imag.teil

// NUR OBEN: Unterschied FFT-inverse FFT im Verteilungsschritt

c=1.0; // Das ist \zeta_{2m}^{0},

s=0.0; // Real- und Imag.teil

for (int j=mstart; j<mend; j++) // Arbeitsschleife

// Das sollte schon bekannt sein

{ // aus der naiven Routine.

rj=real[j];// Die Schleife laeuft ueber m Werte,

ij=imag[j];// veraendert aber 2*m Werte im array

mpj=m+j;

real[j]+=real[mpj];// (12.2.27), rho = 0

imag[j]+=imag[mpj];

real[mpj]=rj-real[mpj];

// Beginn von (12.2.27), rho = 1

imag[mpj]=ij-imag[mpj];

rj =real[mpj]*c-imag[mpj]*s;

imag[mpj]=real[mpj]*s+imag[mpj]*c;

real[mpj]=rj;

ch=c1*c-s1*s; // Additionstheoreme fuer sin/cos

s= c1*s+s1*c;

c=ch;

}

if (RSNMoutput)

{

cout<<"Nach Aufruf split fuer m = "<<m<<

" zwischen "<<mstart<<" und "<<mend+m-1<<"\n";

for (int i=mstart; i<mend+m; i++)

cout.form("%3d : %8.4f + i * %8.4f \n",

i,real,imag);

}

}

}

if (RSNMoutput)

{

cout<<"Ende split fuer n = "<<n<<"\n";

for (int i=0; i<n; i++)

cout.form("%3d : %8.4f + i * %8.4f \n",

i,real,imag);

}

}

void fft(int n, int inv,

// FFT ueber Verteilung und Bitinversion.

Vector real, Vector imag)

// inv =1/-1 : FFT/inverse FFT

{

if (RSNMoutput)

{

cout<<"Aufruf fft fuer n = "<<n<<

" und inv = "<<inv<<"\n";

for (int i=0; i<n; i++)

cout.form("%3d : %8.4f + i * %8.4f \n"

,i,real,imag);

}

split(n, inv, real, imag);

// das macht die Verteilungsschritte

int bitinv; // Der Rest ist Umsortieren der Ergebnisse.

double help; // Dazu wird die Bitinversion benutzt.

for (int i=0; i<n; i++)

// Das ist fuer FFT und inverse FFT gleich.

{

bitinv=bitinverse(i,n);// Kann man auch inline setzen....

if (bitinv>i) // Vermeidet Doppelvertauschung....

{

help=real; // Vertauschen....

real=real[bitinv];

real[bitinv]=help;

help=imag;

imag=imag[bitinv];

imag[bitinv]=help;

}

}

if (inv==(-1)) // Sonderfall inverse FFT

{

for (int i=0; i<n; i++)

{

real/=n;

imag/=n;

}

}

if (RSNMoutput)

{

cout<<"Ende Aufruf fft fuer n = "<<n<<"\n";

for (int i=0; i<n; i++)

cout.form("%3d : %8.4f + i * %8.4f \n",

i,real,imag);

}

}

double fftreal (int n, double t,

// wertet die FFT auf primitive Weise aus

// an der Stelle z=exp(i*t), nur Realteil.

Vector realout, Vector imagout) // FFT output

{ // Setzt voraus, dass eine trigonometrische

// Interpolation der Laenge n vorliegt

double dreal, dimag;

double c, s;

dreal=0.0;

dimag=0.0;

for (int k=0; k<n; k++)

{

c=cos(t*((double) k));

// Auch das kann man durch Drehen machen,

// wenn man sin/cos-Aufrufe sparen will

s=sin(t*((double) k));

dreal+=c*realout[k]-s*imagout[k];

dimag+=s*realout[k]+c*imagout[k];

}

c=cos(t*0.5*((double)n));

// Achtung: es ist an Interpolation gedacht

s=sin(-t*0.5*((double)n));

return(dreal*c-dimag*s);

}

int main()

{

int nmax, chf;

double t, tstep, tnext, s, c;

ofstream pointfile, intfile;

RSNMoutput=1;

pointfile.open("points");

intfile.open("interpol");

do

{

cout<<"Bitte 0/1/2/3/4/5 eingeben"

<<" fuer diverse Funktionen\n";

cin>>chf;

}

while (chf*(6-chf)<0);

cout<<"Bitte Punktzahl (Zweierpotenz) eingeben\n";

cin>>nmax;

Vector realin =NewVector(nmax);

Vector imagin =NewVector(nmax);

for (int i=0; i<nmax; i++)

{

t=((double)i)*2.0*M_PI/((double) nmax);

c=cos(t*0.5*((double)nmax));

// Drehfaktoren fuer Interpolation...

s=sin(t*0.5*((double)nmax)); // siehe (12.2.9)

realin=fct(chf,t);

pointfile<<t<<' '<<realin<<'\n';

realin*=c;

// Drehung vereinfacht wg. Imaginaerteil=0

imagin=s*realin;

}

fft (nmax, -1, realin, imagin);

/* fft (nmax, 1, realin, imagin); */

// zum Testen

// Plotterei

tstep=(0.02*M_PI/((double) nmax));

t=0.0;

for (int i=0; i<nmax; i++)

{

tnext=((double)i+1)*2.0*M_PI/((double) nmax);

while (t <tnext)

{

intfile<<t<<' '<<

fftreal(nmax, t, realin, imagin)<<'\n';

t+=tstep;

}

t=tnext;

}

DeleteVector(imagin);

DeleteVector(realin);

pointfile.close();

intfile.close();

}

[Johannes Buchner]

Lieber "Heinrich der IT-genie 78", auch wenn dein Problem für dich "elit", ist das noch kein Grund, hier unformatierten Code reinzuwerfen. Ich kann so deinen Code nicht lesen und mag nicht indenten! Es wäre besser, wenn du ihn zum Download bereitstellst und konkrete Fragen stellst. Du solltest dein Problem selbst lösen wollen!

siehe http://www.lugbz.org/documents/smart-questions_de.html

[sissy66]

@johannes

DANKESCHÖN.

Bei dieser Bemerkung fällt mir ein Stein vom Herzen.

Liebe Grüße,

sissy66.

[sissy66]

Hi Du krasser Mathematiker,

Du bist übrigens ein Jahr älter als ich, wenn ich die 78 richtig interpretiere.

Du hast Dich also daran getraut, die Fouriertransformation mit imaginären Zahlen zu implementieren. ... ich würde sagen 'Kniefall+Fusskuss'?

Feststellungen als Nicht-Mathematiker-, aber aus Informatikersicht:

- Dieses Programm hat niemals jemals schon ein Compiler gesehen.

- Und Teile Deines Programmes hast Du natürlich nicht mit printf gecheckt oder mal die Klammern von irgendwelchen tollen Tools überprüfen lassen.

... aber das ist nicht Deine Leidenschaft, Deine Leidenschaft ist die Mathe richtig?

-----------------------

Erst mal noch Organisatorisches für Dich:

mathemäßig:

- Diplomarbeiten? kann man aber auch länger planen!!!

- Welches Ausgabe vom Bronstein? benutzt Du eigentlich?

hardware-/softwareseitig:

- Bitte, bitte besorg Dir von irgendwem nen größeren Bildschirm und somit auch nen größeres Fenster.

- Lege Dir bitte einen Editor mit syntax-highlighting+kontetsensitive Hilfe an oder mach in Deinem Emacs? den cpp-mode an (M-x --mode ...).

!!! Programme auscodieren != Mathematische Formeln schreiben !!!

--------------------

Wenn das Programm durchläuft, dann kenne ich ne Menge Leute, die das haben wollen. Das kannste dann patentieren, wenn das sonst so alles von den Beweisen her stimmt und Du das wirklich konsequent mathematisch durchhälts und Dir irgendwen zum Implementieren schnappst.

Ja, es ist durchaus üblich so etwas zu tun! Deswegen diesen Rat hier.

Liebe Grüße,

sissy66.