You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
149 lines
3.5 KiB
149 lines
3.5 KiB
#include <stdlib.h>
|
|
#include <stdio.h>
|
|
#include <string.h>
|
|
#include <ctype.h>
|
|
#include <math.h>
|
|
#include <complex>
|
|
#include<vector>
|
|
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;
|
|
|
|
complex<double> 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<complex<double>> fourierTrafo(fft_param *fp, vector<complex<double>> x) {
|
|
long N = x.size();
|
|
int sign=fp->sign;
|
|
vector<complex<double>> gerPkt(N / 2);
|
|
vector<complex<double>> ungerPkt(N / 2);
|
|
|
|
if (N == 1) {
|
|
|
|
return x;
|
|
} else {
|
|
for (int k = 0; k < N / 2; k++) {
|
|
gerPkt[k] = x[2 * k];
|
|
ungerPkt[k] = x[2 * k + 1];
|
|
}
|
|
|
|
vector<complex<double>> ftger = fourierTrafo(fp, gerPkt);
|
|
vector<complex<double>> ftunger = fourierTrafo(fp, ungerPkt);
|
|
vector<complex<double>> out(N);
|
|
|
|
for (int k = 0; k < N / 2; k++) {
|
|
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<complex<double>> test(fft_param *fp){
|
|
double h=fp->h;
|
|
double tmin=fp->tmin;
|
|
double tmax=fp->tmax;
|
|
double t=tmin;
|
|
int N=fp->N;
|
|
|
|
vector<complex<double>> x(N);
|
|
|
|
for(int k=0;k<N;k++){
|
|
x[k]=sin(tmin+k*h);
|
|
}
|
|
vector<complex<double>> out = fourierTrafo(fp, x);
|
|
fft_param n= *fp;
|
|
n.sign=-1;
|
|
vector<complex<double>> doppelout = fourierTrafo(&n, out);
|
|
FILE *fp2 = fopen("testFFT.csv" ,"w");
|
|
for(int k=0;k<N;k++){
|
|
fprintf(fp2,"%f\t%f\t%f\n",(tmin+k*h),x[k].real(), doppelout[k].real()/(double)N);
|
|
}
|
|
|
|
return x;
|
|
}
|
|
vector<complex<double>> makeS(fft_param *fp){
|
|
double h=fp->h;
|
|
double tmax=fp->tmax;
|
|
double tmin=0;
|
|
double t;
|
|
double tl=fp->tl;
|
|
double alpha=fp->alpha;
|
|
int N=fp->N;
|
|
vector<complex<double>> s(N);
|
|
for(int k=0;k<N;k++){
|
|
t=tmin+k*h;
|
|
s[k]=sin(2*M_PI*alpha*(t-tl));
|
|
}
|
|
return s;
|
|
}
|
|
vector<complex<double>> makeE(fft_param *fp){
|
|
double h=fp->h;
|
|
double tmax=fp->tmax;
|
|
double tmin=0;
|
|
double t;
|
|
double beta=fp->beta;
|
|
double tl=fp->tl;
|
|
double a=fp->a;
|
|
double b=fp->b;
|
|
int N=fp->N;
|
|
vector<complex<double>> e(N);
|
|
vector<complex<double>> s=makeS(fp);
|
|
for(int k=0;k<N;k++){
|
|
t=tmin+k*h;
|
|
e[k]=b*sin(2*M_PI*beta*t)+ 2*a*(rand()-0.5)+s[k];
|
|
//e[k]=b*sin(2*M_PI*beta*t)+s[k];
|
|
}
|
|
return e;
|
|
}
|
|
vector<complex<double>> aufgabe3 (fft_param *fp){
|
|
double h=fp->h;
|
|
double tmin=fp->tmin;
|
|
double tmax=fp->tmax;
|
|
double t=tmin;
|
|
int N=fp->N;
|
|
vector<complex<double>> e=makeE(fp);
|
|
vector<complex<double>> s=makeS(fp);
|
|
vector<complex<double>> Fe=fourierTrafo(fp,e);
|
|
fft_param fps=*fp;
|
|
fps.sign=-1;
|
|
vector<complex<double>> Fs=fourierTrafo(&fps,s);
|
|
vector<complex<double>> FFes(N);
|
|
|
|
for(int k=0;k<N;k++){
|
|
FFes[k]=Fe[k]*conj(Fs[k]);
|
|
|
|
}
|
|
vector<complex<double>> FFFes=fourierTrafo(fp, FFes);
|
|
return FFFes;
|
|
|
|
}
|
|
int main(){
|
|
fft_param fp;
|
|
fp.N= pow(2,13);
|
|
fp.tmin=0;
|
|
fp.tmax = 100;
|
|
fp.h=(fp.tmax-fp.tmin)/fp.N;
|
|
fp.sign=1;
|
|
test(&fp);
|
|
|
|
|
|
}
|