From a4370edc1368c8544a0912b3a3247216001d64ac Mon Sep 17 00:00:00 2001 From: pauline Date: Mon, 13 Jul 2020 12:36:24 +0200 Subject: [PATCH] Hallo --- abgabe6.cpp | 285 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 285 insertions(+) create mode 100644 abgabe6.cpp diff --git a/abgabe6.cpp b/abgabe6.cpp new file mode 100644 index 0000000..7b16682 --- /dev/null +++ b/abgabe6.cpp @@ -0,0 +1,285 @@ +//g++ abgabe6.cpp -std=c++14 -o abgabe6 +// ./abgabe6 +// Pauline Maas, Ariane Ufer + +#include +#include +#include +#include +#include +#include +#include +using namespace std; +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif +typedef struct { + long n; + double b, a, beta, alpha; + double tmin, tmax,h,tl; + int m; + int sign; +} fft_param; + + +//Pseudo-Zufallszahlen-Generator +//Part of a C-program for MT19937, with initialization improved 2002/1/26. +//Coded by Takuji Nishimura and Makoto Matsumoto. + +/* Period parameters */ +#define N 624 +#define M 397 +#define MATRIX_A 0x9908b0dfUL /* constant vector a */ +#define UPPER_MASK 0x80000000UL /* most significant w-r bits */ +#define LOWER_MASK 0x7fffffffUL /* least significant r bits */ + +static unsigned long mt[N]; /* the array for the state vector */ +static int mti=N+1; /* mti==N+1 means mt[N] is not initialized */ + +/* initializes mt[N] with a seed */ +void init_genrand(unsigned long s) +{ + mt[0]= s & 0xffffffffUL; + for (mti=1; mti> 30)) + mti); + /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */ + /* In the previous versions, MSBs of the seed affect */ + /* only MSBs of the array mt[]. */ + /* 2002/01/09 modified by Makoto Matsumoto */ + mt[mti] &= 0xffffffffUL; + /* for >32 bit machines */ + } +} + +/* initialize by an array with array-length */ +/* init_key is the array for initializing keys */ +/* key_length is its length */ +/* slight change for C++, 2004/2/26 */ +void init_by_array(unsigned long init_key[], int key_length) +{ + int i, j, k; + init_genrand(19650218UL); + i=1; j=0; + k = (N>key_length ? N : key_length); + for (; k; k--) { + mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL)) + + init_key[j] + j; /* non linear */ + mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ + i++; j++; + if (i>=N) { mt[0] = mt[N-1]; i=1; } + if (j>=key_length) j=0; + } + for (k=N-1; k; k--) { + mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL)) + - i; /* non linear */ + mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ + i++; + if (i>=N) { mt[0] = mt[N-1]; i=1; } + } + + mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */ +} + +/* generates a random number on [0,0xffffffff]-interval */ +unsigned long genrand_int32(void) +{ + unsigned long y; + static unsigned long mag01[2]={0x0UL, MATRIX_A}; + /* mag01[x] = x * MATRIX_A for x=0,1 */ + + if (mti >= N) { /* generate N words at one time */ + int kk; + + if (mti == N+1) /* if init_genrand() has not been called, */ + init_genrand(5489UL); /* a default initial seed is used */ + + for (kk=0;kk> 1) ^ mag01[y & 0x1UL]; + } + for (;kk> 1) ^ mag01[y & 0x1UL]; + } + y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK); + mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL]; + + mti = 0; + } + + y = mt[mti++]; + + /* Tempering */ + y ^= (y >> 11); + y ^= (y << 7) & 0x9d2c5680UL; + y ^= (y << 15) & 0xefc60000UL; + y ^= (y >> 18); + + return y; +} + +/* generates a random number on [0,1) with 53-bit resolution*/ +double genrand_res53(void) +{ + unsigned long a=genrand_int32()>>5, b=genrand_int32()>>6; + return(a*67108864.0+b)*(1.0/9007199254740992.0); +} +// End of C-program for MT19937 + + + + +complex expo(fft_param *fp,int k) { + int sign = fp->sign; + int n=fp->n; + return exp(sign*2.0 * M_PI * 1i* ((double)(k/n))); + +} + +vector> fourierTrafo(fft_param *fp, vector> x) { + long n = x.size(); //Länge des x vektors + int sign=fp->sign; + vector> gerPkt(n / 2); //Vektoren um x in gerade und ungerade Teile aufzuteilen + vector> ungerPkt(n / 2); //Haben dann Länge n/2 + + if (n == 1) { + + return x; + } else { + //Aufteilen + for (int k = 0; k < n / 2; k++) { + gerPkt[k] = x[2 * k]; + ungerPkt[k] = x[2 * k + 1]; + } + + vector> ftger = fourierTrafo(fp, gerPkt); //rekursives Aufrufen der Funktion um weiter zu Teilen + vector> ftunger = fourierTrafo(fp, ungerPkt); + vector> out(n); + + for (int k = 0; k < n / 2; k++) { + //In out ergebnis nach FFT routine Speichern + if(sign==-1){ + out[k] = (ftger[k] + ftunger[k] * expo(fp,k)); + out[k + n / 2] = (ftger[k] - ftunger[k] *expo(fp,k)); + } + else{ + out[k] = (ftger[k] + ftunger[k] * expo(fp,k)); + out[k + n / 2] = (ftger[k] - ftunger[k] *expo(fp,k)); + } + } + + return out; + } +} + +vector> makeS(fft_param *fp){ + double h=fp->h; //Variabeln aus struct holen + double tmax=fp->tmax; + double tmin=fp->tmin;; + double t; + double tl=fp->tl; + double alpha=fp->alpha; + int n=fp->n; + vector> s(n); //Vektor zum speichern von s + for(int k=0;k= 0 && t-tl <= tmax){ + s[k]=sin(2*M_PI*alpha*(t-tl)); //s berechnen + + } + else{ + s[k]=0; + } + + } + return s; +} +vector> makeE(fft_param *fp){ + double h=fp->h; //Variabeln aus struct holen + double tmax=fp->tmax; + double tmin=fp->tmin; + double t; + double beta=fp->beta; + double tl=fp->tl; + double a=fp->a; + double b=fp->b; + int n=fp->n; + vector> e(n); //Vektor zum speichern von e + vector> s=makeS(fp); //Vektor s wird benötigt + //FILE *fp3 = fopen("H10_2_30.txt" ,"w"); + for(int k=0;k> aufgabe3 (fft_param *fp){ + double h=fp->h; //Varabeln aus Struct holen + double tmin=fp->tmin; + double tmax=fp->tmax; + double t=tmin; + int n=fp->n; + vector> e=makeE(fp); //e und s berechnen + vector> s=makeS(fp); + vector> Fe=fourierTrafo(fp,e); //e Fouriertransformieren + + vector> Fs=fourierTrafo(fp,s); //s Fouriertransformieren + vector> FFes(n); //Vektor machen zum speichern + + for(int k=0;k> FFFes=fourierTrafo(&fps, FFes); //Rücktransformation des Produkts + return FFFes; + //init_genrand((long)seed); +} +void plotAufgabe3(fft_param *fp){ + int n = fp->n; + double tmin=fp->tmin; + double h=fp->h; + vector> es = aufgabe3(fp); + FILE *fp4 =fopen("es4.csv","w"); + for (int k=0;k