--- libscidavis/src/FFT.cpp 2017-07-14 04:54:53.000000000 -0300
+++ /home/fellype/compilando/scidavis-github/libscidavis/src/FFT.cpp 2017-07-19 09:36:40.000000000 -0300
@@ -50,6 +50,12 @@
{
init();
setDataFromCurve(curveTitle);
+ // intersperse 0 imaginary components
+ double* tmp=new double[2*d_n];
+ memset(tmp,0,2*d_n*sizeof(double));
+ for (size_t i=0; i<d_n; ++i) tmp[2*i]=d_y[i];
+ delete [] d_y;
+ d_y=tmp;
}
void FFT::init ()
@@ -63,133 +69,13 @@
d_sampling = 1.0;
}
-QList<Column *> FFT::fftCurve()
+QList<Column *> FFT::fftTable()
{
- int i2;
- int n2 = d_n/2;
double *amp = new double[d_n];
- double *result = new double[2*d_n];
-
- if(!amp || !result)
- {
- QMessageBox::critical((ApplicationWindow *)parent(), tr("SciDAVis") + " - " + tr("Error"),
- tr("Could not allocate memory, operation aborted!"));
- d_init_err = true;
- return QList<Column *>();
- }
-
- double df = 1.0/(double)(d_n*d_sampling);//frequency sampling
- double aMax = 0.0;//max amplitude
- QList<Column *> columns;
- if(!d_inverse)
- {
- d_explanation = tr("Forward") + " " + tr("FFT") + " " + tr("of") + " " + d_curve->title().text();
- columns << new Column(tr("Frequency"), SciDAVis::Numeric);
-
- gsl_fft_real_workspace *work=gsl_fft_real_workspace_alloc(d_n);
- gsl_fft_real_wavetable *real=gsl_fft_real_wavetable_alloc(d_n);
-
- if(!work || !real)
- {
- QMessageBox::critical((ApplicationWindow *)parent(), tr("SciDAVis") + " - " + tr("Error"),
- tr("Could not allocate memory, operation aborted!"));
- d_init_err = true;
- return QList<Column *>();
- }
-
- gsl_fft_real_transform(d_y, 1, d_n, real,work);
- gsl_fft_halfcomplex_unpack (d_y, result, 1, d_n);
- gsl_fft_real_wavetable_free(real);
- gsl_fft_real_workspace_free(work);
- }
- else
- {
- d_explanation = tr("Inverse") + " " + tr("FFT") + " " + tr("of") + " " + d_curve->title().text();
- columns << new Column(tr("Time"), SciDAVis::Numeric);
-
- gsl_fft_real_unpack (d_y, result, 1, d_n);
gsl_fft_complex_wavetable *wavetable = gsl_fft_complex_wavetable_alloc (d_n);
gsl_fft_complex_workspace *workspace = gsl_fft_complex_workspace_alloc (d_n);
- if(!workspace || !wavetable)
- {
- QMessageBox::critical((ApplicationWindow *)parent(), tr("SciDAVis") + " - " + tr("Error"),
- tr("Could not allocate memory, operation aborted!"));
- d_init_err = true;
- return QList<Column *>();
- }
-
- gsl_fft_complex_inverse (result, 1, d_n, wavetable, workspace);
- gsl_fft_complex_wavetable_free (wavetable);
- gsl_fft_complex_workspace_free (workspace);
- }
-
- if (d_shift_order)
- {
- for(unsigned i=0; i<d_n; i++)
- {
- d_x[i] = (i-n2)*df;
- int j = i + d_n;
- double aux = result[i];
- result[i] = result[j];
- result[j] = aux;
- }
- }
- else
- {
- for(unsigned i=0; i<d_n; i++)
- d_x[i] = i*df;
- }
-
- for(unsigned i=0;i<d_n;i++)
- {
- i2 = 2*i;
- double real_part = result[i2];
- double im_part = result[i2+1];
- double a = sqrt(real_part*real_part + im_part*im_part);
- amp[i]= a;
- if (a > aMax)
- aMax = a;
- }
-
- // ApplicationWindow *app = (ApplicationWindow *)parent();
-
- columns << new Column(tr("Real"), SciDAVis::Numeric);
- columns << new Column(tr("Imaginary"), SciDAVis::Numeric);
- columns << new Column(tr("Amplitude"), SciDAVis::Numeric);
- columns << new Column(tr("Angle"), SciDAVis::Numeric);
- for (unsigned i=0;i<d_n;i++)
- {
- i2 = 2*i;
- columns.at(0)->setValueAt(i, d_x[i]);
- columns.at(1)->setValueAt(i, result[i2]);
- columns.at(2)->setValueAt(i, result[i2+1]);
- if (d_normalize)
- columns.at(3)->setValueAt(i, amp[i]/aMax);
- else
- columns.at(3)->setValueAt(i, amp[i]);
- columns.at(4)->setValueAt(i, atan(result[i2+1]/result[i2]));
- }
- delete[] amp;
- delete[] result;
- columns.at(0)->setPlotDesignation(SciDAVis::X);
- columns.at(1)->setPlotDesignation(SciDAVis::Y);
- columns.at(2)->setPlotDesignation(SciDAVis::Y);
- columns.at(3)->setPlotDesignation(SciDAVis::Y);
- columns.at(4)->setPlotDesignation(SciDAVis::Y);
- return columns;
-}
-
-QList<Column *> FFT::fftTable()
-{
- int i;
- int rows = d_table->numRows();
- double *amp = new double[rows];
-
- gsl_fft_complex_wavetable *wavetable = gsl_fft_complex_wavetable_alloc (rows);
- gsl_fft_complex_workspace *workspace = gsl_fft_complex_workspace_alloc (rows);
-
if(!amp || !wavetable || !workspace)
{
QMessageBox::critical((ApplicationWindow *)parent(), tr("SciDAVis") + " - " + tr("Error"),
@@ -198,18 +84,18 @@
return QList<Column *>();
}
- double df = 1.0/(double)(rows*d_sampling);//frequency sampling
+ double df = 1.0/(double)(d_n*d_sampling);//frequency sampling
double aMax = 0.0;//max amplitude
QList<Column *> columns;
if(!d_inverse)
{
columns << new Column(tr("Frequency"), SciDAVis::Numeric);
- gsl_fft_complex_forward (d_y, 1, rows, wavetable, workspace);
+ gsl_fft_complex_forward (d_y, 1, d_n, wavetable, workspace);
}
else
{
columns << new Column(tr("Time"), SciDAVis::Numeric);
- gsl_fft_complex_inverse (d_y, 1, rows, wavetable, workspace);
+ gsl_fft_complex_inverse (d_y, 1, d_n, wavetable, workspace);
}
gsl_fft_complex_wavetable_free (wavetable);
@@ -217,11 +103,11 @@
if (d_shift_order)
{
- int n2 = rows/2;
- for(i=0; i<rows; i++)
+ int n2 = d_n/2;
+ for(int i=0; i<int(d_n); i++)
{
d_x[i] = (i-n2)*df;
- int j = i + rows;
+ int j = i + d_n;
double aux = d_y[i];
d_y[i] = d_y[j];
d_y[j] = aux;
@@ -229,11 +115,11 @@
}
else
{
- for(i=0; i<rows; i++)
+ for(size_t i=0; i<d_n; i++)
d_x[i] = i*df;
}
- for(i=0; i<rows; i++)
+ for(size_t i=0; i<d_n; i++)
{
int i2 = 2*i;
double a = sqrt(d_y[i2]*d_y[i2] + d_y[i2+1]*d_y[i2+1]);
@@ -246,7 +132,7 @@
columns << new Column(tr("Imaginary"), SciDAVis::Numeric);
columns << new Column(tr("Amplitude"), SciDAVis::Numeric);
columns << new Column(tr("Angle"), SciDAVis::Numeric);
- for (i=0; i<rows; i++)
+ for (size_t i=0; i<d_n; i++)
{
int i2 = 2*i;
columns.at(0)->setValueAt(i, d_x[i]);
@@ -270,9 +156,7 @@
void FFT::output()
{
QList<Column *> columns;
- if (d_graph && d_curve)
- columns = fftCurve();
- else if (d_table)
+ if (d_y)
columns = fftTable();
if (!columns.isEmpty())
|
