#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; }