如何在C++中正确使用带有向量的FFTW

我是C++的初学者，在实习期间，我必须使用fftw库来获得信号的傅立叶变换。

事实上，一开始，我的程序使用python，然后我试图将其翻译成C++。然而，它不起作用，我真的不明白为什么。。。我真的不习惯指针之类的东西，所以也许这就是问题所在。首先，我想知道为什么程序(python和C++中的程序(不能给出相同的结果？此外，如果有人能向我解释如何将FFTW与矢量一起使用，我将不胜感激？我想更好地理解事物。非常感谢：(

python代码

import numpy as np
import matplotlib.pyplot as plt
from math import cos, sin, sqrt, pi, exp as cos, sin, sqrt, pi, exp
#parameters (suprathresold case)
eps = 0.05
f = 0.215
mu = 1.1
D = 0.001
DeltaT = 0.001
period = int(1.0 / f / DeltaT)
num_periods_4_transient = 10
num_periods_4_spectrum = 100

#timewindow = 2 * period * DeltaT #quite long to be independant of the transient 
l_transient = num_periods_4_transient * period
l_measure = num_periods_4_spectrum * period
#num_points = int(timewindow/DeltaT)
N = 100 # number of realization (number of neurons) should be at the end 10^8
DeltaT_ratio = 100
#signal that modulates the firing rate 
#s = [sin(2*3.14*f*t) for t in T]
ts_transient = np.arange(l_transient)*DeltaT
ts_measure = np.arange(l_measure)*DeltaT
s = np.sin(2.0 * np.pi * f *  ts_transient)
#Euler's method with FFT 
v0 = np.zeros(N)
v = np.zeros((l_transient,N))
samples = np.random.normal(0, 1, (l_transient,N))
DeltaF = 1.0 / (l_measure * DeltaT)
#print "running transient..."
for i in range(1,l_transient):
v[i,:] = v[i-1,:] + DeltaT *(-v[i-1,:]+ mu + eps*s[i-1])+ sqrt(2*D*DeltaT)*samples[i,:]
mask_sp = v[i,:] > 1
v[i, mask_sp ] = 0
v0 = v[-1,:]
v = np.zeros((l_measure,N))
v[0,:] = v0

s = np.sin(2.0 * np.pi * f *  ts_measure)
samples = np.random.normal(0, 1, (l_measure,N))
l_compteur = l_measure // DeltaT_ratio + 1
compteur=np.zeros((l_compteur, N))
DeltaT_prime = DeltaT * DeltaT_ratio
#print l_measure,l_measure // DeltaT_ratio
#print "running for measure..."
for i in range(1,l_measure):
v[i,:] = v[i-1,:] + DeltaT *(-v[i-1,:]+ mu + eps*s[i-1])+ sqrt(2*D*DeltaT)*samples[i,:]
mask_sp = v[i,:] > 1
compteur[ int(i / DeltaT_ratio), mask_sp ] = 1.0   / DeltaT_prime        
v[i, mask_sp ] = 0                
l_freq = int(l_compteur / 2) + 1
spectrum_matrix = np.abs(np.fft.rfft(compteur,axis=0)*DeltaT_prime)**2 
powerspectra = np.mean(spectrum_matrix,axis=1)
powerspectra /= (l_measure*DeltaT)
powerspectra[0] = 0.0
plt.plot(np.arange(l_freq)*DeltaF, powerspectra )
plt.show()

C++代码

#include <iostream>
#include <cmath>
#include <complex>
#include <fstream>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <fstream>
#include <fftw3.h>
#define PI 3.1415926535897932384626433832795
#include <vector>
using namespace std;
double s(double x, double );
double AWGN_generator(void);
std::string double2str(const double x)
{
stringstream sss;
sss<<x;
return sss.str();
}
int main()
{   
//    srand (time(NULL));
//parameters
double eps = 0.05;
double f_1 = 0.215 ;
//double f_2 = 0.33;
double mu = 1.1 ;
double D = 0.001 ;
double deltaT = 0.001;
long period = 1.0 / f_1 / deltaT;
int N = 1000; //# number of realization (number of neurons) should be at the end 10^8
long DeltaT_ratio = 100;
double DeltaT_prime  = DeltaT_ratio*deltaT;
int num_periods_4_transient = 10;
int num_periods_4_spectrum = 20;
long l_transient = num_periods_4_transient * period;
long l_measure = num_periods_4_spectrum * period;
long l_compteur = l_measure / DeltaT_ratio + 1 ;
long l_freq = int(l_compteur / 2) + 1;
double deltaf = 1.0 / (l_compteur * DeltaT_prime);
vector<double> powerspectrum(l_freq,0);
cout<<l_measure<<l_measure / DeltaT_ratio<<endl;
double a = 1.0;
ofstream m_out ;
string outfile = "/Users/Pierrot/powerspectrum_test.txt";
m_out.open(outfile.c_str());

vector< double > compteur(l_compteur,0);
std::vector< std::complex<double> > out(l_freq,0);
fftw_plan p;
p = fftw_plan_dft_r2c_1d(l_compteur, &compteur[0], reinterpret_cast<fftw_complex*>(&out[0]), FFTW_MEASURE);

//cout<<"test1"<<endl;
// Transient    
vector< double > v(l_transient,0.0 );
vector<double> Time(l_transient, 0) ;
for (int i = 0 ;  i < l_transient ; i++ )
{   Time[i] = deltaT*i ;
}
int pos_1 ;
for (pos_1 = 1 ; pos_1<l_transient ; pos_1 ++)
{
double t = Time[pos_1-1] ;
v[pos_1] = v[pos_1-1] + deltaT*(-v[pos_1-1]+ mu + eps*s(t, f_1)) + sqrt(2*D*deltaT)*AWGN_generator();
if (v[pos_1]>1)
{
v[pos_1]=0 ;
}
}
vector< double > v_bis(l_measure, v[l_transient-1]);
//    v_bis[0]=v[l_transient-1];    
for( int k=0; k<N ; k++)
{
//copy last value of transient in the new v
v_bis[0] = v_bis[l_measure-1];
double DeltaT_prime = deltaT * DeltaT_ratio ;
vector<double> Time_bis(l_measure) ;

for (int i = 0 ;  i < l_measure ; i++ )
{   Time_bis[i] = deltaT*i ;
}
//cout<<"test2"<<endl;    
for (int pos_1_bis = 1 ; pos_1_bis<l_measure ; pos_1_bis ++)
{
double t = Time_bis[pos_1_bis-1] ;
v_bis[pos_1_bis] = v_bis[pos_1_bis-1] + deltaT*(-v_bis[pos_1_bis-1]+ mu + eps*s(deltaT*(pos_1_bis-1), f_1)) + sqrt(2*D*deltaT)*AWGN_generator();
//cout<<v_bis[pos_1_bis]<<endl;
if (v_bis[pos_1_bis]>1)
{
v_bis[pos_1_bis]=0 ;
compteur[pos_1_bis/DeltaT_ratio] = 1.0/DeltaT_prime ;
//                       cout<<pos_1_bis/DeltaT_ratio<<endl;
}
}


fftw_execute(p);

for (long m(0);m<l_freq;m++)
{
powerspectrum[m] +=  (real( out[m] * conj(out[m]) ) * DeltaT_prime*DeltaT_prime) / (l_measure*deltaT*N);
} 

}

//powerspectrum = powerspectrum *DeltaT_prime*DeltaT_prime/(l_measure*deltaT*N)

fftw_destroy_plan(p);
fftw_cleanup();
for (long m(1);m<l_freq;m++)
{
m_out<<deltaf*m <<"t"<<powerspectrum[m]<<"n";
}
m_out.close();
printf("ca a marche test.txt");
return 0;
}



double s(double x, double f)
{
return sin(2*PI*f*x);
}

double AWGN_generator()
{/* Generates additive white Gaussian Noise samples with zero mean and a standard
deviation of 1. */
double temp1;
double temp2;
double result;
int p;
p = 1;
while( p > 0 )
{
temp2 = ( rand() / ( (double)RAND_MAX ) ); /*  rand() function generates an
integer between 0 and  RAND_MAX,
which is defined in stdlib.h.
*/
if ( temp2 == 0 )
{// temp2 is >= (RAND_MAX / 2)
p = 1;
}// end if
else
{// temp2 is < (RAND_MAX / 2)
p = -1;
}// end else
}// end while()
temp1 = cos( ( 2.0 * (double)PI ) * rand() / ( (double)RAND_MAX ) );
result = sqrt( -2.0 * log( temp2 ) ) * temp1;
return result;        // return the generated random sample to the caller
}