離散フーリエ変換

前回のフーリエ変換のプログラムだと入力されるデータの間隔が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);
}