Skip to content

GitLab

Commit 51631c47 authored by Blue Bird's avatar Blue Bird
Browse files

Aenderungen wegen FFT Variable

parent 3e818820
Hide whitespace changes
Inline Side-by-side
common_source/mathlib/fouriertransform.cc 0 → 100644
View file @ 51631c47
// solid laser for solid state laser simulation:
// Copyright (C) 2017 ASLD GmbH
// This program is not a free software;
// do not redistribute it and/or modify it
/////////////////////////////////////////////////
/*
* fouriertransform.cc
*
* Created on: Apr 23, 2015
* Author: zhabiz
*/
#include <string>
#include <iostream>
#include <fstream>
#include <math.h>
#include <complex>
#include "fouriertransform.h"
using namespace std;
# define PI 3.1415926535897
complex<double>* FourierTransform::temp = NULL;
int FourierTransform::size_temp = 0;
int FourierTransform::nk = 0;
FourierTransform::FourierTransform(int n_){
nk = n_;
}
/*
* to check whether a number is 2^n
*
*
*/
bool FourierTransform::isPowerOfTwo(unsigned int x)
{
return (!(x == 0) && !(x & (x - 1)));
}
/*
* to ckeck whether a number is even
*
*/
bool FourierTransform::isEven(int x){
if (x % 2)
return false;
else
return true;
}
/*
* to copy vectors
*
*/
void FourierTransform::scalarProduct(complex<double> *x, double scalar){
int i;
for(i=0; i<nk; i++){
x[i].real( x[i].real() * scalar );
x[i].imag( x[i].imag() * scalar );
}
}
/*
* to copy vectors
*
*/
void FourierTransform::copy(complex<double> *x, complex<double> *y){
int i;
for(i=0; i<nk; i++)
y[i] = x[i];
}
/*
* to calculate intensity of the field or spectrum
*
*/
void FourierTransform::calcIntensity(complex<double> *x, complex<double> *intensity){
int i;
for(i=0; i<nk; i++){
intensity[i].real( x[i].real()*x[i].real() + x[i].imag()*x[i].imag() );
intensity[i].imag( 0.);
}
}
void FourierTransform::calcIntensityField(const char* filename, complex<double> *x, complex<double> *intensity, double *xAxis, int n){
////////////////////////////////////////////////////////
// Calculate the intensity of a beam in time domain
//////////////////////////////////////////////////////
iFFT(x, n);
////////////////////////////
//-shift iFFT
////////////////////////////
//fft.copy(stretched_beam, stretched_beam_iFFT);
FFTShift(x, n);
calcIntensity(x, intensity);
print_vector_to_file( filename, intensity, xAxis);
}
/*
*
*
*/
void FourierTransform::print_vector(const char *title, complex<double> *x)
{
int i;
std::cout.precision(6);
std::cout << "\n" ;
std::cout << title << " " << nk <<std::endl;
for(i=0; i<nk; i++ )
std::cout << x[i].real() << " " << x[i].imag() <<std::endl;
std::cout << "\n" ;
return;
}
/*
*
*
*/
void FourierTransform::print_vector_to_file(const char *filename, complex<double> *x, double *xAxis)
{
int i;
ofstream OUT;
OUT.precision(4);
OUT.setf(std::ios::fixed,std::ios::floatfield);
OUT.open(filename, std::ios :: out);
for(i=0; i<nk; i++ ){
OUT << xAxis[i] << " " << x[i].real() << " " << x[i].imag() << endl;
//cout << xAxis[i] << " " << x[i].real() << " " << x[i].imag() <<std::endl;
}
OUT.close();
return;
}
/*
*
*
*
*/
void FourierTransform::resizeTemp(int n) {
if (!isPowerOfTwo(n)) {
std::cout << " Error: number of points must be a power of 2. It is:" << n << std::endl;
return;
}
if(size_temp < n) {
size_temp = n;
delete[] temp;
temp = new complex<double>[n];
}
}
/*
* FFt shift
* equal to FFTshift of Matlab
*
*/
void FourierTransform::FFTShift(complex<double> *v, int n){
int n2;
n2 = n / 2;
complex<double> tmp;
for (int i = 0; i < n2; i++) // swap what can be swapped
{
tmp = v[i];
v[i] = v[i+n2];
v[i+n2] = tmp;
}
if (n & 1) // odd n, shift the rest
{
tmp = v[n-1];
v[n-1] = v[0];
for (int i = 1; i < n2; i++)
{
v[i-1] = v[i];
}
v[n2-1] = tmp;
}//end if
} //end function FFTShift
/*
* fast Fourier transformation of a 1D vector
* *x : contains real part
* *y : contains imaginary part
*
*
*/
void FourierTransform::FFT(complex<double> *v, int n) {
resizeTemp(n);
FFT(v, n, temp);
FFTShift(v, n);
}
//private member function
void FourierTransform::FFT(complex<double> *v, int n, complex<double> *tmp) {
if(n>1) { /* otherwise, do nothing and return */
int k,m;
complex<double> z, w, *vo, *ve;
ve = tmp; vo = tmp+n/2;
for(k=0; k<n/2; k++) {
ve[k] = v[2*k];
vo[k] = v[2*k+1];
}
FFT( ve, n/2, v ); /* FFT on even-indexed elements of v[] */
FFT( vo, n/2, v ); /* FFT on odd-indexed elements of v[] */
for(m=0; m<n/2; m++) {
w.real( cos(2*PI*m/(double)n) );
w.imag( -sin(2*PI*m/(double)n) );
z.real( w.real()*vo[m].real() - w.imag()*vo[m].imag() ); /* Re(w*vo[m]) */
z.imag( w.real()*vo[m].imag() + w.imag()*vo[m].real() ); /* Im(w*vo[m]) */
v[ m ].real( ve[m].real() + z.real() );
v[ m ].imag( ve[m].imag() + z.imag() );
v[m+n/2].real( ve[m].real() - z.real() );
v[m+n/2].imag( ve[m].imag() - z.imag() );
}
}
}//end function FFT
/*
* fast inverse Fourier transformation of a 1D vector
*
*
*
*
*/
void FourierTransform::iFFT(complex<double> *v, int n){
resizeTemp(n);
FFTShift(v, n);
iFFT(v,n,temp);
for(int i = 0; i<n;i++){
v[i].real( v[i].real() / n );
v[i].imag( v[i].imag() / n );
}
}
/*
* fast inverse Fourier transformation of a 1D vector
*
*/
void FourierTransform::iFFT(complex<double> *v, int n, complex<double> *tmp){
if(n>1) {
int k,m; complex<double> z, w, *vo, *ve;
ve = tmp; vo = tmp+n/2;
for(k=0; k<n/2; k++) {
ve[k] = v[2*k];
vo[k] = v[2*k+1];
}
iFFT( ve, n/2, v ); /* FFT on even-indexed elements of v[] */
iFFT( vo, n/2, v ); /* FFT on odd-indexed elements of v[] */
for(m=0; m<n/2; m++) {
w.real( cos(2*PI*m/(double)n) );
w.imag( sin(2*PI*m/(double)n) );
z.real( w.real() * vo[m].real() - w.imag() * vo[m].imag() ); /* Re(w*vo[m]) */
z.imag( w.real() * vo[m].imag() + w.imag() * vo[m].real() ); /* Im(w*vo[m]) */
v[ m ].real( ve[m].real() + z.real() );
v[ m ].imag( ve[m].imag() + z.imag() );
v[m+n/2].real( ve[m].real() - z.real() );
v[m+n/2].imag( ve[m].imag() - z.imag() );
}
} //end if
}
/*
* fast Fourier transformation of a 1D vector
* *x : contains real part
* *y : contains imaginary part
*
*
*/
void FourierTransform::FFT_efficient(int dir,int m,double *x,double *y) {
long nn,i,i1,j,k,i2,l,l1,l2;
double c1,c2,tx,ty,t1,t2,u1,u2,z;
// Calculate the number of points
nn = 1;
for (i=0;i<m;i++)
nn *= 2;
//cout << " n : " << nn <<std::endl;
i2 = nn >> 1;
j = 0;
for (i=0;i<nn-1;i++) {
if (i < j) {
tx = x[i];
ty = y[i];
x[i] = x[j];
y[i] = y[j];
x[j] = tx;
y[j] = ty;
}
k = i2;
while (k <= j) {
j -= k;
k >>= 1;
}
j += k;
} // end for loop
// Compute the FFT
c1 = -1.0;
c2 = 0.0;
l2 = 1;
for (l=0;l<m;l++){
l1 = l2;
l2 <<= 1;
u1 = 1.0;
u2 = 0.0;
for (j=0;j<l1;j++) {
for (i=j;i<nn;i+=l2) {
i1 = i + l1;
t1 = u1 * x[i1] - u2 * y[i1];
t2 = u1 * y[i1] + u2 * x[i1];
x[i1] = x[i] - t1;
y[i1] = y[i] - t2;
x[i] += t1;
y[i] += t2;
}
z = u1 * c1 - u2 * c2;
u2 = u1 * c2 + u2 * c1;
u1 = z;
} // end for loop
c2 = sqrt((1.0 - c1) / 2.0);
if (dir == 1)
c2 = -c2;
c1 = sqrt((1.0 + c1) / 2.0);
} // end outer for loop
// Scaling for forward transform
if (dir == 1) {
for (i=0;i<nn;i++) {
x[i] /= (double)nn;
y[i] /= (double)nn;
}
} // end if
}//end function FFT_efficient
common_source/mathlib/fouriertransform.h 0 → 100644
View file @ 51631c47
// solid laser for solid state laser simulation:
// Copyright (C) 2017 ASLD GmbH
// This program is not a free software;
// do not redistribute it and/or modify it
/////////////////////////////////////////////////
/*
* fouriertransform.h
*
* Author: zhabiz
*/
#ifndef FOURIERTRANSFORM_H_
#define FOURIERTRANSFORM_H_
#include <complex>
using namespace std;
class FourierTransform {
public:
FourierTransform(int n_);
//recursive implementation
static void FFT(complex<double> *v, int n);
static void iFFT(complex<double> *v, int n);
void FFT_efficient(int dir,int m, double *x,double *y);
static void FFTShift(complex<double> *v, int n);
static bool isPowerOfTwo(unsigned int x);
static bool isEven(int x);
static void scalarProduct(complex<double> *x, double scalar);
static void copy(complex<double> *x, complex<double> *y);
static void calcIntensity(complex<double> *x, complex<double> *intensity);
static void calcIntensityField(const char* filename, complex<double> *x, complex<double> *intensity, double *xAxis, int n);
static void print_vector(const char *title, complex<double> *x);
static void print_vector_to_file(const char* filename, complex<double> *x, double *xAxis);
static int nk;
private:
static void FFT(complex<double> *v, int n, complex<double> *tmp);
static void iFFT(complex<double> *v, int n, complex<double> *tmp);
static void resizeTemp(int n);
static complex<double> *temp;
static int size_temp;
};
#endif /* FOURIERTRANSFORM_H_ */
program2D/source/extemp/variableFFT.cc 0 → 100644
View file @ 51631c47
/**********************************************************************************
* Copyright 2010 Christoph Pflaum
* Department Informatik Lehrstuhl 10 - Systemsimulation
* Friedrich-Alexander Universität Erlangen-Nürnberg
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
**********************************************************************************/
// ------------------------------------------------------------
//
// variableFFT.cc
//
// ------------------------------------------------------------
#include <string>
#include <stdexcept>
#include "../../../program/source/mympi.h" // von 3D UGBlocks
#include "../abbrevi2D.h"
//#include "../../../linAlgExptemp/source/linAlg.h"
//#include "../parameter.h"
#include "../math_lib/math_lib.h"
#include "../grid/elements2D.h"
#include "../grid/parti2D.h"
#include "../grid/ug2D.h"
#include "../grid/blockgrid2D.h"
#include "../grid/marker2D.h"
//#include "parallel.h"
#include "extemp2D.h"
#include "variable2D.h"
#include "variableFFT.h"
#include "variable2D_cc.h"
#include "update2D.h"
#include "print_var2D_vtk_cc.h"
#include "../grid/examples_ug2D.h"
#include "../../../common_source/mathlib/fouriertransform.h"
#include "../interpol/interpolTwoD.h"
int VariableFFT::powerTwo(int t) {
int n=1;
for(int l=0;l<t;++l) n=n*2;
return n;
}
void VariableFFT::operator=(VariableFFT& v) {
Variable2D<std::complex<double> >*pv = static_cast<Variable2D<std::complex<double> >* >(&v);
static_cast<Variable2D<std::complex<double> >* >(this)->operator=(*pv);
}
void VariableFFT::interpolate(Variable2D<std::complex<double> >& v,
std::complex<double> defaultInterpolation) {
bool createInterpolator = (assignmentBlockgrid==NULL);
if(createInterpolator==false) {
if(assignmentBlockgrid->getId() != v.Give_blockgrid()->getId()) {
delete assignmentInterpolator;
createInterpolator=true;
}
}
if(createInterpolator) {
assignmentBlockgrid = v.Give_blockgrid();
assignmentInterpolator
= new Interpolate_on_structured_2Dgrid(Nx+1,Ny+1,
D2vector(ug->MinimumX(),ug->MinimumY()),
D2vector(ug->MaximumX(),ug->MaximumY()),
*assignmentBlockgrid);
}
assignmentInterpolator->interpolate(v,data_rectangles[0],defaultInterpolation);
Update_back(0);
}
VariableFFT::VariableFFT(int tx, int ty, Rectangle& ugRectangle) :
Variable2D<std::complex<double> >(
*(blockgrid = new Blockgrid2D(&ugRectangle,powerTwo(tx)-1,powerTwo(ty)-1))) {
ownBlockGrid = true;
init();
}
VariableFFT::VariableFFT(VariableFFT* other) :
Variable2D<std::complex<double> >(*(other->Give_blockgrid())) {
ownBlockGrid = false;
init();
}
VariableFFT::VariableFFT(const VariableFFT& other) :
Variable2D<std::complex<double> >(*(other.Give_blockgrid())) {
ownBlockGrid = false;
init();
}
void VariableFFT::init() {
Nx = blockgrid->Give_Nx_rectangle(0);
Ny = blockgrid->Give_Ny_rectangle(0);
startValueY = new complex<double>[Ny+1];
FourierTransform intF(MAX(Nx+1,Ny+1));
assignmentInterpolator = NULL;
assignmentBlockgrid = NULL;
}
VariableFFT::~VariableFFT() {
if(ownBlockGrid) blockgrid;
delete[] startValueY;
}
void VariableFFT::FFT() {
Update<rectangleEl>(0);
for(int j=0;j<=Ny;++j) {
std::complex<double>* startValue = &(data_rectangles[0][( Nx+1 ) *j]);
FourierTransform::FFT(startValue,Nx+1);
}
for(int i=0;i<=Nx;++i) {
for(int j=0;j<=Ny;++j) {
startValueY[j] = data_rectangles[0][i + ( Nx+1 ) *j];
}
FourierTransform::FFT(startValueY,Ny+1);
for(int j=0;j<=Ny;++j) {
data_rectangles[0][i + ( Nx+1 ) *j] = startValueY[j];
}
}
Update_back(0);
}
void VariableFFT::inversFFT() {
Update<rectangleEl>(0);
for(int i=0;i<=Nx;++i) {
for(int j=0;j<=Ny;++j) {
startValueY[j] = data_rectangles[0][i + ( Nx+1 ) *j];
}
FourierTransform::iFFT(startValueY,Nx+1);
for(int j=0;j<=Ny;++j) {
data_rectangles[0][i + ( Nx+1 ) *j] = startValueY[j];
}
}
for(int j=0;j<=Ny;++j) {
complex<double>* startValue = &(data_rectangles[0][( Nx+1 ) *j]);
FourierTransform::iFFT(startValue,Nx+1);
}
Update_back(0);