Page MenuHomec4science
Files F67487391

surface_generator_filter_fft.cpp
No OneTemporary
Actions

File Metadata

Created
Sat, Jun 22, 16:56

surface_generator_filter_fft.cpp
View Options

/* -------------------------------------------------------------------------- */
#include "surface_generator_filter_fft.hh"
#include "surface_statistics.hh"
/* -------------------------------------------------------------------------- */
#include <iostream>
#include <cmath>
#include <fstream>
#include <fftw3.h>
/* -------------------------------------------------------------------------- */
SurfaceGeneratorFilterFFT::SurfaceGeneratorFilterFFT(const std::string & inputfile):
SurfaceGenerator(inputfile),SurfaceGeneratorFilter(inputfile){
}
/* -------------------------------------------------------------------------- */
SurfaceGeneratorFilterFFT::SurfaceGeneratorFilterFFT(){
}
/* -------------------------------------------------------------------------- */
complex SurfaceGeneratorFilterFFT::computeFractalPSDProfile(int i, int j){
Real q = sqrt(i*i+j*j);
complex value(Real(0.),Real(0.));
if (q < q0) value = 0.;
else if (q <= q1) value = 1.;
else if (q > q2) value = 0.;
else value = pow(q/q1,-2.*(hurst+1));
return value;
}
/* -------------------------------------------------------------------------- */
void SurfaceGeneratorFilterFFT::buildCorrelationInput(){
int n = surface->size();
// all entries are set to zero automatically
Surface<Real> power_input(n,1.);
std::cout << "rms: " << rms << " , hurst: " << hurst << " , q0: " << q0 << " , q1: " << q1 << " , q2: " << q2 << std::endl;
power_input(0,0) = 0.;
for (int i = 1 ; i < n/2 ; ++i)
for (int j = 1 ; j < n/2 ; ++j){
complex value = computeFractalPSDProfile(i,j);
power_input(i,j) = value.real();
power_input(n-i,n-j) = value.real();
power_input(n-i,j) = value.real();
power_input(i,n-j) = value.real();
}
// external line (i or j == 0)
for (int i = 1 ; i < n/2 ; ++i){
complex value = computeFractalPSDProfile(i,0);
power_input(i,0) = value.real();
power_input(0,i) = value.real();
power_input(0,n-i) = value.real();
power_input(n-i,0) = value.real();
value = computeFractalPSDProfile(i,n/2);
power_input(i,n/2) = value.real();
power_input(n-i,n/2) = value.real();
power_input(n/2,i) = value.real();
power_input(n/2,n-i) = value.real();
}
{
power_input(n/2,n/2) = computeFractalPSDProfile(n/2,n/2).real();
power_input(0,n/2) = computeFractalPSDProfile(0,n/2).real();
power_input(n/2,0) = computeFractalPSDProfile(n/2,0).real();
}
power_input.dumpToParaview("power_spectrum_input");
Real average_spectrum = SurfaceStatistics::computeAverageAmplitude(*white_noise_FFT);
std::cout << "white noise average spectrum is " << average_spectrum << std::endl;
// compute fourier transform of filter coefficients
h_coeff_power = new SurfaceComplex<Real>(n,1.);
Real C = average_spectrum;
for (int i = 0 ; i < n ; ++i)
for (int j = 0 ; j < n ; ++j){
h_coeff_power->at(i,j) = sqrt(power_input(i,j)/C);
}
}
/* -------------------------------------------------------------------------- */
void SurfaceGeneratorFilterFFT::computeFilterCoefficients(){
buildCorrelationInput();
}
/* -------------------------------------------------------------------------- */
/* Parse an input file */
void SurfaceGeneratorFilterFFT::parseInputFile(std::ifstream & fp) {
std::string line;
while (! fp.eof() )
{
getline (fp,line);
size_t pos = line.find("#"); //remove the comment
line = line.substr(0, pos);
trim(line); // remove unnecessary spaces
/* Read all options one by one */
readInput(line, "seed" , random_seed);
readInput(line, "gridsize", gridsize );
readInput(line, "rms" , rms );
readInput(line, "l0" , q0 );
readInput(line, "l1" , q1 );
readInput(line, "l2" , q2 );
readInput(line, "hurst" , hurst );
}
return;
}

Event Timeline

c4science · Help