離散フーリエ変換
前回のフーリエ変換のプログラムだと入力されるデータの間隔が1.0でないとうまく動かないという問題がありました。
それと、周波数を求める箇所が間違えてました。
ってことで改良したものを
#include#include using namespace std; class Dft { public: int m_sum; double* m_cos_coefficient; double* m_sin_coefficient; double m_length; public: Dft(); ~Dft(); void computing(double* f, int sum, double dx); void display(void); void reconstruct(void); void reconstruct(double start, double end, int number_divid); }; Dft::Dft() { m_sum = 0; m_cos_coefficient=NULL; m_sin_coefficient=NULL; } Dft::~Dft() { if( m_cos_coefficient!=NULL){ delete m_cos_coefficient; } if( m_sin_coefficient!=NULL){ delete m_sin_coefficient; } } void Dft::computing(double* f, int sum, double dx) { this->m_sum=sum; this->m_cos_coefficient = new double[sum]; this->m_sin_coefficient = new double[sum]; this->m_length = dx*(double)(sum); for( int i = 0; i < sum; i ++){ this->m_cos_coefficient[i]=0.0; this->m_sin_coefficient[i]=0.0; double T = 2*(double)(i)*M_PI/this->m_length; double tmp_cos = 0.0; double tmp_sin = 0.0; // To substitute sum for integral operation for( int j = 0; j < sum; j ++){ double x = (double)(j)*dx; tmp_cos +=f[j]*cos(dx*T)*dx; tmp_sin +=f[j]*sin(dx*T)*dx; } this->m_cos_coefficient[i] = tmp_cos/this->m_length; this->m_sin_coefficient[i] = tmp_sin/this->m_length; } } void Dft::display(void) { for( int i = 0; i < this->m_sum; i ++){ double T = 2*(double)(i)*M_PI/this->m_length; cout << T << " "; cout << this->m_cos_coefficient[i] << " "; cout << this->m_sin_coefficient[i] << endl; } } void Dft::reconstruct(void) { double x = this->m_length/(this->m_sum); cout << x << endl; for( int i = 0; i < this->m_sum; i ++){ double value = 0.0; double w = (double)(i)*x; for( int j = 0; j < this->m_sum; j ++){ double T = 2*(double)(j)*M_PI/this->m_length; value += m_cos_coefficient[j]*cos(w*T)+ m_sin_coefficient[j]*sin(w*T); } cout << (double)i*x << " " << value << endl; } } void Dft::reconstruct(double start, double end, int number_divid) { double dx = (end - start)/(double)number_divid; for( int i = 0; i < number_divid; i ++){ double value = 0.0; for( int j = 0; j < this->m_sum; j ++){ double T = 2*(double)(j)*M_PI/this->m_length; double w = start + (double)(i)*(dx); value += m_cos_coefficient[j]*cos(w*T)+ m_sin_coefficient[j]*sin(w*T); } cout << start + dx*i << " " << value << endl; } } int main(void) { Dft dft; double dx = 3.0; double tmp[] = { 0, 0, 0, 0, 1, 0, 0, 0, 0, }; dft.computing(tmp, sizeof(tmp)/8, dx); //dft.display(); //dft.reconstruct(); dft.reconstruct(0,(sizeof(tmp)/8)*dx,18); }