\r
\r
\r
+ FFT1::FFT1()\r
+ {\r
+ construct_default();\r
+ }\r
+ FFT1::FFT1(int wid)\r
+ {\r
+ construct_default();\r
+ set(wid);\r
+ }\r
+ FFT1::FFT1(const Matrix &source)\r
+ {\r
+ construct_default();\r
+ set(source);\r
+ }\r
+ FFT1::~FFT1()\r
+ {\r
+ release();\r
+ }\r
+ void FFT1::construct_default()\r
+ {\r
+ width_=0;\r
+ img_spc=0; frq_spc=0;\r
+ }\r
+ void FFT1::release()\r
+ {\r
+ if(img_spc!=0) { free(img_spc); img_spc = 0; }\r
+ if(frq_spc!=0) { free(frq_spc); frq_spc = 0; }\r
+ width_=0;\r
+ }\r
+ void FFT1::set(int wid)\r
+ {\r
+ if(img_spc==0 || frq_spc==0 || (width_ != wid) ) {\r
+ release();\r
+ img_spc = (fftw_complex*)malloc(sizeof(fftw_complex) * wid);\r
+ frq_spc = (fftw_complex*)malloc(sizeof(fftw_complex) * wid);\r
+ }\r
+ width_ = wid;\r
+ x_zero = wid + 1;\r
+ left = (wid)/2;\r
+ }\r
+ void FFT1::set(const Matrix &source)\r
+ {\r
+ if(source.getRows()==1) {\r
+ set(source.getCols());\r
+ for(int x=0; x<width_; x++) {\r
+ img_spc[x][REAL] = source(1, x+1);\r
+ img_spc[x][IMAGINARY] = 0.0;\r
+ }\r
+ }\r
+ }\r
+\r
+ double FFT1::getDC()\r
+ {\r
+ return sqrt(frq_spc[0][REAL]*frq_spc[0][REAL]+frq_spc[0][IMAGINARY]*frq_spc[0][IMAGINARY]);\r
+ }\r
+ double FFT1::setDC(double l)\r
+ {\r
+ frq_spc[0][REAL] = l;\r
+ frq_spc[0][IMAGINARY] = 0;\r
+ return l;\r
+ }\r
+\r
+\r
+ void FFT1::fft()\r
+ {\r
+ fftw_plan plan = fftw_plan_dft_1d(width_, img_spc, frq_spc, FFTW_FORWARD, FFTW_ESTIMATE);\r
+ fftw_execute(plan);\r
+ normalizeFFT();\r
+ }\r
+ void FFT1::ifft()\r
+ {\r
+ fftw_plan plan = fftw_plan_dft_1d(width_, frq_spc, img_spc, FFTW_BACKWARD, FFTW_ESTIMATE);\r
+ fftw_execute(plan);\r
+ }\r
+ void FFT1::normalizeFFT()\r
+ {\r
+ double num = width_;\r
+ for (int i=0; i<num; i++) {\r
+ frq_spc[i][REAL] /= num;\r
+ frq_spc[i][IMAGINARY] /= num;\r
+ }\r
+ }\r
+ void FFT1::getImage(Matrix &target)\r
+ {\r
+ double rl, im;\r
+ int h, v;\r
+ if( (target.getCols()!=width_)) {\r
+ target.release();\r
+ target.set(1, width_);\r
+ }\r
+ for (int x=0; x<width_; x++) {\r
+ rl = img_spc[x][REAL];\r
+ im = img_spc[x][IMAGINARY];\r
+ target(1, x+1) = sqrt(rl*rl+im*im);\r
+ }\r
+ }\r
+ void FFT1::getSpectrum(Matrix &target, double gamma)\r
+ {\r
+ const int right = width_%2==0 ? left : left+1;\r
+ double rl, im;\r
+ Color col;\r
+ if( (target.getCols()!=width_) ) {\r
+ target.release();\r
+ target.set(1, width_);\r
+ }\r
+ int h, v = 1;\r
+ for (int x=0; x<right; x++) {\r
+ h = x+left + 1;\r
+ rl = frq_spc[x][REAL];\r
+ im = frq_spc[x][IMAGINARY];\r
+ target(v, h) = pow( sqrt(rl*rl+im*im), gamma );\r
+ }\r
+ for (int x=right; x<width_; x++) {\r
+ h = x-right + 1;\r
+ rl = frq_spc[x][REAL];\r
+ im = frq_spc[x][IMAGINARY];\r
+ target(v, h) = pow( sqrt(rl*rl+im*im), gamma );\r
+ }\r
+ }\r
+\r
+\r
\r
FFT2::FFT2()\r
{\r
double cumulativeLog2NormalDistibution(double log_x, double octave_mu, double octave_sigma);\r
\r
\r
+ class FFT1 {\r
+ protected:\r
+ int width_;\r
+ int left, x_zero;\r
+ fftw_complex *img_spc;\r
+ fftw_complex *frq_spc;\r
+ void construct_default();\r
+ public:\r
+\r
+ // Initialize\r
+ FFT1();\r
+ FFT1(int width);\r
+ FFT1(const Matrix &source);\r
+ ~FFT1();\r
+ void release();\r
+ void set(int wid);\r
+ void set(const Matrix &source);\r
+\r
+ // Accesser to elements\r
+ double getDC();\r
+ double setDC(double l);\r
+\r
+ // Core FFT Execution\r
+ void fft();\r
+ void ifft();\r
+ void normalizeFFT();\r
+\r
+ void getImage(Matrix &absolute);\r
+ void getSpectrum(Matrix &absolute, double gamma = 1.0);\r
+ //void getImage(Matrix &reali, Matrix &imagi);\r
+ //void getSpectrum(Matrix &reali, Matrix &imagi);\r
+\r
+ };\r
+\r
class FFT2 {\r
protected:\r
int width_, height_;\r
#include <psychlops.h>\r
using namespace Psychlops;\r
\r
+\r
+void psychlops_main() {\r
+ Canvas cnvs(1024, 768, Canvas::window);\r
+\r
+ FFT1 fftworkspace;\r
+ Matrix source_mat, filtered_mat, mat;\r
+ source_mat.load("voltage_spec_vaio_ie.csv");\r
+ source_mat.transpose();\r
+\r
+ //fftworkspace.set(64,64);\r
+ //fftworkspace.makeNoise();\r
+ fftworkspace.set(source_mat);\r
+ fftworkspace.fft();\r
+ fftworkspace.getSpectrum(filtered_mat, 1.0);\r
+\r
+ int q = filtered_mat.getCols() / 4;\r
+ Interval col;\r
+ mat = filtered_mat(1, q*2+3<=col<=q*3);\r
+ mat.transpose();\r
+ mat.save("voltage_spec.csv");\r
+\r
+ double mean = 0.0;\r
+ for(int i=0; i<source_mat.getCols(); i++) {\r
+ mean += source_mat(1, i+1);\r
+ }\r
+ mean /= source_mat.getCols();\r
+ int counter = 0;\r
+ for(int i=1; i<source_mat.getCols()-3; i++) {\r
+ if( source_mat(1, i+1)<mean\r
+ && source_mat(1, i+2)<mean\r
+ && source_mat(1, i+3)>mean\r
+ && source_mat(1, i+4)>mean\r
+ ) { counter++; }\r
+ }\r
+ std::cout << counter << std::endl;\r
+}\r
+\r
+\r
+/*\r
const int S = 128;\r
\r
void psychlops_main() {\r
\r
}\r
+*/\r