OpenGL de プログラミング - OpenCV::離散フーリエ変換::線形性
現在地 >> メニュー >> サンプルコード::OpenCV >> OpenCV::パワースペクトルの描画 >> OpenCV::離散フーリエ変換::線形性

問題


以下の操作を行え。
  1. 入力データを2つ作る
  2. 入力データそれぞれを離散フーリエし、その結果を合成する
  3. 入力データを合成してから離散フーリエする
  4. [2],[3]で得られた結果のパワースペクトルを出力し、離散フーリエの線形性を確認する。

※線形性
F[ f1 + f2 ] = F[f1] + F[f2] = F1 + F2
つまり、
 合成してから離散フーリエしても、離散フーリエしてから合成しても
結果は同じであるという事。

答え


#include <iostream>
#include <cmath>
#include <vector>
#include <algorithm>
#include <cv.h>
#include <highgui.h>

/********* 各種外部変数 ***********/
const CvSize WindowSize = {300,300};
const int SLICE = 8;
const int ELEMENT = SLICE+1;
const double PI_OVER_180 = 0.0174532925; // π/180
const double unit_steps = 360/SLICE;

/*********** プロトタイプ宣言 ***********/
CvMat *GetMat(double *data);

CvMat *GetMatFromVec(std::vector<double>::iterator start,int nums);

void Show_Result(const CvMat *dataM);
IplImage *Get_Power_Spectrum(const CvMat *dataM);
double GET_MAX(const double *power_spectrum);



/************** ここからメイン関数 **********************/
int main()
{

	using namespace std;
	
	//先に領域をある程度確保
	vector<double> data_1(ELEMENT),data_2(ELEMENT),mix_data(ELEMENT);

	double angle;
	for(int loop = 0; loop < ELEMENT; ++loop)
	{
		angle = unit_steps *loop;
		data_1[loop]=sin(angle *  PI_OVER_180);  //sin X
		data_2[loop]=sin( (angle + 10)*PI_OVER_180 ); // sin(X+10)
		mix_data[loop]= ( data_1.at(loop) + data_2.at(loop) );
	}


	CvMat *data_1_M = GetMatFromVec(data_1.begin(),(ELEMENT));
	CvMat *data_2_M = GetMatFromVec(data_2.begin(),(ELEMENT));
	CvMat *mix_data_M = GetMatFromVec(mix_data.begin(),(ELEMENT));


	/* 離散フーリエ変換 */
	cvDFT(data_1_M,data_1_M,CV_DXT_FORWARD,data_1_M->rows);
	cvDFT(data_2_M,data_2_M,CV_DXT_FORWARD,data_2_M->rows);
	cvDFT(mix_data_M,mix_data_M,CV_DXT_FORWARD,mix_data_M->rows);

	CvMat *Mat_plus_Mat = cvCreateMat(1,ELEMENT,CV_64FC2);
	cvAdd(data_1_M,data_2_M,Mat_plus_Mat);

	IplImage *imgA = Get_Power_Spectrum(mix_data_M);
	IplImage *imgB = Get_Power_Spectrum(Mat_plus_Mat);


	cvNamedWindow("data mix before DFT",CV_WINDOW_AUTOSIZE);
	cvShowImage("data mix before DFT",imgA);
	

	cvNamedWindow("data mix after DFT",CV_WINDOW_AUTOSIZE);
	cvShowImage("data mix after DFT",imgB);


	/* メモリ解放 */
	cvReleaseImage(& imgA);
	cvReleaseImage(& imgB);
	
	cvReleaseMat(&data_1_M);
	cvReleaseMat(&data_2_M);
	cvReleaseMat(&mix_data_M);
	cvReleaseMat(&Mat_plus_Mat);

	cvWaitKey(0);
	cvDestroyAllWindows();

	return EXIT_SUCCESS;
}


/********** データから計算用に行列を作成 ***********/
CvMat *GetMatFromVec(std::vector<double>::iterator start,int nums)
{
	using namespace std;
	CvMat *complex = cvCreateMat(1,nums, CV_64FC2);
	CvMat Re = cvMat(1,ELEMENT, CV_64FC1,&(*start));
	CvMat *Im = cvCreateMat(1,ELEMENT, CV_64FC1);
	cvSetZero(Im);

	cvMerge(&Re,Im,NULL,NULL,complex);
	cvReleaseMat(&Im);

	return complex;
}

/**********パワースペクトルの画像を取得*************/
IplImage *Get_Power_Spectrum(const CvMat *dataM)
{
	using std::cout;

	IplImage *imgA = cvCreateImage(WindowSize,IPL_DEPTH_8U,3);//ウィンドウの初期化
	cvSet (imgA, cvScalarAll (255), 0);

	double Power_Spectrum[ELEMENT];//パワースペクトル用配列
	for(int loop = 0; loop < ELEMENT; ++loop) //パワースペクトルの計算
	{
		Power_Spectrum[loop] = pow(dataM->data.db[loop*2],2) + pow(dataM->data.db[loop*2+1],2); 
	}
	
	
	double Max = GET_MAX(&Power_Spectrum[0]);  //パワースペクトルの最大値を取得

	/*** 出力用にスケール変換 ***/
	for(int loop = 0; loop < ELEMENT; ++loop)
	{
		Power_Spectrum[loop] = (Power_Spectrum[loop]/Max) * imgA->height;
	}


	//横幅決定
	int unit_width = cvRound( imgA->width / static_cast<double>(ELEMENT) );


	/*右端の余白からシフト分を計算*/
	int space_shift = cvRound*1*0.5);


	/**** 描画 ***/
	CvScalar color={0,255,0};
	bool flag = false;
	for (int j = 0; j < ELEMENT; j++)
	{
		if (j == 0 || j == (ELEMENT - 1) )
		{
			flag = true; //両端を赤く塗るフラグ
		}

		cvLine(imgA,cvPoint(j * unit_width+space_shift , imgA->height),cvPoint (j * unit_width +space_shift ,
			imgA->height - cvRound (Power_Spectrum[j])), flag ? CV_RGB(255,0,0) : color ,3);

		flag = false;
	}

	return imgA;

}


/*** 配列の最大値を取得 ***/
double GET_MAX(const double *power_spectrum)
{
	double temp = power_spectrum[0];
	for(int loop = 1; loop < ELEMENT; ++loop)
	{
		if(temp < power_spectrum[loop])
		{
			temp = power_spectrum[loop];
		}

	}

	return temp;

}