3
// Copyright (c) 1998-2010 by The VoxBo Development Team
5
// This file is part of VoxBo
7
// VoxBo is free software: you can redistribute it and/or modify it
8
// under the terms of the GNU General Public License as published by
9
// the Free Software Foundation, either version 3 of the License, or
10
// (at your option) any later version.
12
// VoxBo is distributed in the hope that it will be useful, but
13
// WITHOUT ANY WARRANTY; without even the implied warranty of
14
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15
// General Public License for more details.
17
// You should have received a copy of the GNU General Public License
18
// along with VoxBo. If not, see <http://www.gnu.org/licenses/>.
20
// For general information on VoxBo, including the latest complete
21
// source code and binary distributions, manual, and associated files,
22
// see the VoxBo home page at: http://www.voxbo.org/
24
// original version written by Dongbo Hu
28
#include "fitOneOverF.h"
30
double a2fixed = 0; // global variable defined for two-parameter fitting
32
/*************************************************************************
33
* expb_f is the function to be fitted. fittingVar is a gsl vector
34
* which includes three variables to be fitted.
35
* In GSL documentation, this variable is named "x",
36
* which caused some confusion with the input X values.
38
* params is a variable of data type which includes the
39
* number of elements in the input X and Y, input values of X and Y,
40
* and standard deviation at each point.
41
*************************************************************************/
42
int expb_f (const gsl_vector *fittingVar , void *params,
45
size_t n = ((struct data *)params)->n;
46
double *inputX = ((struct data *)params)->inputX;
47
double *inputY = ((struct data *)params)->inputY;
48
double *sigma = ((struct data *) params)->sigma;
49
double a0 = gsl_vector_get (fittingVar, 0);
50
double a1 = gsl_vector_get (fittingVar, 1);
51
double a2 = gsl_vector_get (fittingVar, 2);
53
for (size_t i = 0; i < n; i++) {
54
/* Model Yi = 1.0 / (a0 * (Xi + a2)) + a1 */
56
double Yi = 1.0 / (a0 * (t + a2)) + a1;
57
gsl_vector_set (f, i, (Yi - inputY[i])/sigma[i]);
63
/*************************************************************************
64
* expb_df includes the first order partialderivatives of each parameter *
65
*************************************************************************/
66
int expb_df (const gsl_vector * fittingVar, void *params,
69
size_t n = ((struct data *)params)->n;
70
double *inputX = ((struct data *) params)->inputX;
71
double *sigma = ((struct data *) params)->sigma;
72
double a0 = gsl_vector_get (fittingVar, 0);
73
double a2 = gsl_vector_get (fittingVar, 2);
75
for (size_t i = 0; i < n; i++) {
76
/* Jacobian matrix J(i,j) = dfi / dxj, */
77
/* where fi = (Yi - yi)/sigma[i], */
78
/* Yi = 1.0 / (a0 * (Xi + a2)) + a1 */
79
/* and the xj are the parameters (a0, a1, a2) */
81
double s = 1.0 / sigma[i];
82
double tmp = -1.0 /(a0 * (t + a2));
84
gsl_matrix_set (J, i, 0, tmp / a0 * s);
85
gsl_matrix_set (J, i, 1, s);
86
gsl_matrix_set (J, i, 2, tmp * s / (t + a2));
92
/************************************************************************
93
* expb_fdf combines expb_f anf expb_df together. *
94
************************************************************************/
95
int expb_fdf (const gsl_vector * fittingVar, void *params,
96
gsl_vector * f, gsl_matrix * J)
98
expb_f (fittingVar, params, f);
99
expb_df (fittingVar, params, J);
103
/* Print out three parameters */
104
void print_state (size_t iter, gsl_multifit_fdfsolver * s)
106
printf ("iter: %3d parameters are: %15.8f %15.8f %15.8f |f(x)| = %g\n",
107
(int)iter, gsl_vector_get (s->x, 0), gsl_vector_get (s->x, 1),
108
gsl_vector_get (s->x, 2), gsl_blas_dnrm2 (s->f));
111
/*************************************************************************************
112
* curvefit() is the main function to calculate the three fitting parameters.
113
* This function requires 8 arguments:
114
* #1: vector of x values for fitting
115
* #2: vector of y values for fitting
116
* #3: vector of sigma values at each point (not available in IDL's curvefit function
117
* #4: initial guess for first parameter (steepness)
118
* #5: initial guess for second parameter (y-offset)
119
* #6: initial guess for third parameter (x-shift)
120
* #7: output file's name
121
* #8: print mode (true: print out a summary on stdout; false: no any screen printout
122
**************************************************************************************/
123
VB_Vector *curvefit(VB_Vector *xVec, VB_Vector *yVec, VB_Vector *sigmaVec, double var3min,
124
double var1_init, double var2_init, double var3_init,
125
const char *outputFile, bool printFlag)
127
/* fittingVec will include 8 elements: fittingStatus, # of iteration,
128
* three fitting parameters' value, three deviation values
129
* If status is 1.0, it means fitting is successful. 0 means fitting isn't successful.*/
130
VB_Vector *fittingVec = new VB_Vector(8);
131
size_t N = xVec->getLength();
132
const gsl_multifit_fdfsolver_type *T;
133
gsl_multifit_fdfsolver *s;
139
gsl_matrix *covar = gsl_matrix_alloc (p, p);
141
double inputX[N], inputY[N], sigma[N];
142
struct data d = {n, inputX, inputY, sigma};
143
gsl_multifit_function_fdf f;
145
double fitting_init[3] = { var1_init, var2_init, var3_init };
146
gsl_vector_view fittingVar = gsl_vector_view_array (fitting_init, p);
147
const gsl_rng_type * type;
151
type = gsl_rng_default;
152
r = gsl_rng_alloc (type);
161
for (size_t i = 0; i < n; i++) {
162
inputX[i] = xVec->getElement(i);
163
inputY[i] = yVec->getElement(i);
165
/* Set sigma to be 0.1 at each point. This value is completely arbitrary.
166
* It won't change the fitting variables' values, but the deviation of each
167
* fitting variable (the value after "+/-") will be affected. */
168
sigma[i] = sigmaVec->getElement(i);
171
T = gsl_multifit_fdfsolver_lmsder;
172
s = gsl_multifit_fdfsolver_alloc (T, n, p);
173
gsl_multifit_fdfsolver_set (s, &f, &fittingVar.vector);
176
print_state (iter, s);
180
status = gsl_multifit_fdfsolver_iterate (s);
182
printf ("status = %s\n", gsl_strerror (status));
183
print_state (iter, s);
186
// Make sure var3 is always larger than minimum. If not, call curverfit12()
187
if (gsl_vector_get (s->x, 2) <= var3min) {
188
double var3 = var3min * 0.99;
189
printf("******************************************************************\n");
190
printf("x-shift minimum value reached.\n");
191
printf("Arbitrarily set x-shift to %.5f and fit the other two ...\n", var3);
192
printf("******************************************************************\n");
193
gsl_multifit_fdfsolver_free (s);
194
return curvefit12(xVec, yVec, sigmaVec, var3,gsl_vector_get (s->x, 0),
195
gsl_vector_get (s->x, 1), outputFile, printFlag);
201
status = gsl_multifit_test_delta (s->dx, s->x, 1e-4, 1e-4);
203
while (status == GSL_CONTINUE && iter < 500);
205
gsl_multifit_covar (s->J, 0.0, covar);
207
#define FIT(i) gsl_vector_get(s->x, i)
208
#define ERR(i) sqrt(gsl_matrix_get(covar, i, i))
211
printf("\nElements in covariance matrix are:\n");
212
gsl_matrix_fprintf (stdout, covar, "%g");
214
printf("\nparameter #1: steepness = %.5f +/- %.5f\n", FIT(0), ERR(0));
215
printf("parameter #2: y-offset = %.5f +/- %.5f\n", FIT(1), ERR(1));
216
printf("parameter #3: x-shift = %.5f +/- %.5f\n", FIT(2), ERR(2));
217
printf ("status = %s\n\n", gsl_strerror (status));
220
// Write the fitting result to a file (if output filename is available)
221
if (outputFile != '\0') {
222
FILE * f1 = fopen(outputFile, "w");
223
fprintf (f1, ";VB98\n");
224
fprintf (f1, ";REF1\n");
225
fprintf (f1, "; GSL nonlinear least-squares fitting\n");
226
string myTime, myDate;
227
maketimedate(myTime, myDate); // maketimedate() is defined in libvoxbo/vbutil.cpp
228
fprintf (f1, "; Date created: %s_%s\n", myTime.data(), myDate.data());
229
fprintf (f1, "; Initial guess: %f, %f, %f\n", var1_init, var2_init, var3_init);
230
fprintf (f1, "; Number of iteration: %d\n", (int)iter);
231
fprintf (f1, "; Status = %s\n", gsl_strerror (status));
232
fprintf (f1, "; Deviation: %f, %f, %f\n\n", ERR(0), ERR(1), ERR(2));
233
fprintf (f1, "%f\n%f\n%f\n", FIT(0), FIT(1), FIT(2));
237
// Set elements in fittingVec. If status is 0 (success), set elements one by one.
238
// Otherwise set all elements to be zero.
240
fittingVec->setElement(0, 1.0);
241
fittingVec->setElement(1, iter);
242
fittingVec->setElement(2, FIT(0));
243
fittingVec->setElement(3, FIT(1));
244
fittingVec->setElement(4, FIT(2));
246
fittingVec->setElement(5, ERR(0));
247
fittingVec->setElement(6, ERR(1));
248
fittingVec->setElement(7, ERR(2));
251
fittingVec->setAll(0);
253
gsl_multifit_fdfsolver_free (s);
257
/* Two parameter version of expb_f */
258
int expb_f12 (const gsl_vector *fittingVar, void *params,
261
size_t n = ((struct data *)params)->n;
262
double *inputX = ((struct data *)params)->inputX;
263
double *inputY = ((struct data *)params)->inputY;
264
double *sigma = ((struct data *) params)->sigma;
266
double a0 = gsl_vector_get (fittingVar, 0);
267
double a1 = gsl_vector_get (fittingVar, 1);
269
for (size_t i = 0; i < n; i++) {
270
/* Model Yi = 1.0 / (a0 * (Xi + var3)) + a1 */
271
double t = inputX[i];
272
double Yi = 1.0 / (a0 * (t + a2fixed)) + a1;
273
gsl_vector_set (f, i, (Yi - inputY[i])/sigma[i]);
279
/***************************************************************************
280
* expb_df12 includes the first order partialderivatives of each parameter *
281
***************************************************************************/
282
int expb_df12 (const gsl_vector * fittingVar, void *params, gsl_matrix *J)
284
size_t n = ((struct data *)params)->n;
285
double *inputX = ((struct data *) params)->inputX;
286
double *sigma = ((struct data *) params)->sigma;
287
double a0 = gsl_vector_get (fittingVar, 0);
289
for (size_t i = 0; i < n; i++) {
290
/* Jacobian matrix J(i,j) = dfi / dxj, */
291
/* where fi = (Yi - yi)/sigma[i], */
292
/* Yi = 1.0 / (a0 * (Xi + a2fixed)) + a1 */
293
/* and the xj are the parameters (a0, a1, a2fixed) */
294
double t = inputX[i];
295
double s = 1.0 / sigma[i];
296
double tmp = -1.0 /(a0 * (t + a2fixed));
297
gsl_matrix_set (J, i, 0, s * tmp / a0 );
298
gsl_matrix_set (J, i, 1, s);
304
/************************************************************************
305
* expb_fdf combines expb_f anf expb_df together. *
306
************************************************************************/
307
int expb_fdf12 (const gsl_vector * fittingVar, void *params,
308
gsl_vector * f, gsl_matrix * J)
310
expb_f12 (fittingVar, params, f);
311
expb_df12 (fittingVar, params, J);
316
/* Print out two parameters */
317
void print_state12 (size_t iter, gsl_multifit_fdfsolver * s)
319
printf ("iter: %3d parameters are: %15.8f %15.8f |f(x)| = %g\n", (int)iter,
320
gsl_vector_get (s->x, 0), gsl_vector_get (s->x, 1), gsl_blas_dnrm2 (s->f));
323
/* curvefit12() fits the input X and Y vectors based on this equation:
324
* Y = 1.0 / [var1 * (X + var3min * 0.99)] + var2
325
* Note that var3min * 0.99 is a fixed arbitrary value to make sure
326
* that when 1/f function is expanded to a longer range (totalReps * TR / 2000),
327
* the initial value wouldn't be negative and enormously large. */
328
VB_Vector *curvefit12 (VB_Vector *xVec, VB_Vector *yVec, VB_Vector *sigmaVec,
329
double var3, double var1_init, double var2_init,
330
const char *outputFile, bool printFlag)
333
/* fittingVec will include 8 elements: fittingStatus, # of iteration,
334
* three fitting parameters' value, three deviation values
335
* If status is 1.0, it means fitting is successful. 0 means fitting isn't successful.*/
336
VB_Vector *fittingVec = new VB_Vector(8);
337
size_t N = xVec->getLength();
338
const gsl_multifit_fdfsolver_type *T;
339
gsl_multifit_fdfsolver *s;
345
gsl_matrix *covar = gsl_matrix_alloc (p, p);
347
double inputX[N], inputY[N], sigma[N];
348
struct data d = {n, inputX, inputY, sigma};
349
gsl_multifit_function_fdf f;
351
double fitting_init[2] = {var1_init, var2_init};
352
gsl_vector_view fittingVar = gsl_vector_view_array (fitting_init, p);
353
const gsl_rng_type * type;
357
type = gsl_rng_default;
358
r = gsl_rng_alloc (type);
367
for (size_t i = 0; i < n; i++) {
368
inputX[i] = xVec->getElement(i);
369
inputY[i] = yVec->getElement(i);
370
/* Set sigma to be 0.1 at each point. This value is completely arbitrary.
371
* It won't change the fitting variables' values, but the deviation of each
372
* fitting variable (the value after "+/-") will be affected. */
373
sigma[i] = sigmaVec->getElement(i);
376
T = gsl_multifit_fdfsolver_lmsder;
377
s = gsl_multifit_fdfsolver_alloc (T, n, p);
378
gsl_multifit_fdfsolver_set (s, &f, &fittingVar.vector);
381
print_state12(iter, s);
385
status = gsl_multifit_fdfsolver_iterate (s);
388
printf ("status = %s\n", gsl_strerror (status));
389
print_state12(iter, s);
394
status = gsl_multifit_test_delta (s->dx, s->x, 1e-4, 1e-4);
396
while (status == GSL_CONTINUE && iter < 500);
398
gsl_multifit_covar (s->J, 0.0, covar);
400
#define FIT(i) gsl_vector_get(s->x, i)
401
#define ERR(i) sqrt(gsl_matrix_get(covar, i, i))
404
printf("\nElements in covariance matrix are:\n");
405
gsl_matrix_fprintf (stdout, covar, "%g");
407
printf("\nparameter #1: steepness = %.5f +/- %.5f\n", FIT(0), ERR(0));
408
printf("parameter #2: y-offset = %.5f +/- %.5f\n", FIT(1), ERR(1));
409
printf("parameter #3: x-shift = %.5f (fixed)\n", var3);
410
printf ("status = %s\n\n", gsl_strerror (status));
413
// Write the fitting result to a file (if output filename is available)
414
if (outputFile != '\0') {
415
FILE * f1 = fopen(outputFile, "w");
416
fprintf (f1, ";VB98\n");
417
fprintf (f1, ";REF1\n");
418
fprintf (f1, "; GSL nonlinear least-squares fitting\n");
419
string myTime, myDate;
420
maketimedate(myTime, myDate); // maketimedate() is defined in libvoxbo/vbutil.cpp
421
fprintf (f1, "; Date created: %s_%s\n", myTime.data(), myDate.data());
422
fprintf (f1, "; Initial guess: %f, %f\n", var1_init, var2_init);
423
fprintf (f1, "; Number of iteration: %d\n", (int)iter);
424
fprintf (f1, "; Status = %s\n", gsl_strerror (status));
425
fprintf (f1, "; Deviation: %f, %f\n\n", ERR(0), ERR(1));
426
fprintf (f1, "%f\n%f\n%f\n", FIT(0), FIT(1), var3);
430
// Set elements in fittingVec. If status is 0 (success),
431
// set elements one by one, Otherwise set all elements to zero.
433
fittingVec->setElement(0, 1.0);
434
fittingVec->setElement(1, iter);
435
fittingVec->setElement(2, FIT(0));
436
fittingVec->setElement(3, FIT(1));
437
fittingVec->setElement(4, var3);
438
fittingVec->setElement(5, ERR(0));
439
fittingVec->setElement(6, ERR(1));
440
fittingVec->setElement(7, 0);
443
fittingVec->setAll(0);
445
gsl_multifit_fdfsolver_free (s);
449
/******************************************************************************************
450
* wrapper function for fitting calculation. It accepts filenames for x, y, sigma values. *
451
* It also accepts filename for initial guess. *
452
******************************************************************************************/
453
VB_Vector * curvefit(const char *xFilename, const char *yFilename, const char *sigmaFilename,
454
double var3min, const char *initFile, const char *outputFile)
456
// Check xFilename format
457
VB_Vector * xVec = new VB_Vector();
458
string xString(xFilename);
459
if (xVec->ReadFile(xString)) {
460
printf("Invalid file format for X: %s\n", xFilename);
464
// Check yFilename format
465
VB_Vector * yVec = new VB_Vector();
466
string yString(yFilename);
467
if (yVec->ReadFile(yString)) {
468
printf("Invalid file format for Y: %s\n", yFilename);
472
// Check sigmaFilename format
473
VB_Vector * sigmaVec = new VB_Vector();
474
string sigmaString(sigmaFilename);
475
if (sigmaVec->ReadFile(sigmaString)) {
476
printf("Invalid file format for sigma: %s\n", sigmaFilename);
480
// Check initFile format
481
VB_Vector * initVec = new VB_Vector();
482
string initString(initFile);
483
if (initVec->ReadFile(initString)) {
484
printf("Invalid file format for initialization: %s\n", initFile);
488
/* Read the file which has initial values for three parameters */
489
double var1_init = initVec->getElement(0);
490
double var2_init = initVec->getElement(1);
491
double var3_init = initVec->getElement(2);
493
return curvefit(xVec, yVec, sigmaVec, var3min, var1_init, var2_init, var3_init, outputFile);
497
/******************************************************************************************
498
* Wrapper function for fitting calculation. It accepts filenames for x, y, sigma values. *
499
* input sigma values is a constant at each point. *
500
******************************************************************************************/
501
VB_Vector * curvefit(const char * xFilename, const char * yFilename, double inputSigma, double var3min,
502
double var1_init, double var2_init, double var3_init, const char *outputFile)
504
/* Read input file for X and Y values */
505
VB_Vector * xVec = new VB_Vector();
506
string xString(xFilename);
507
if (xVec->ReadFile(xString)) {
508
printf("Invalid file format for X: %s\n", xFilename);
512
VB_Vector * yVec = new VB_Vector(yFilename);
513
string yString(yFilename);
514
if (yVec->ReadFile(yString)) {
515
printf("Invalid file format for Y: %s\n", yFilename);
519
/* Read input file for sigma values */
520
int length = xVec->getLength();
521
VB_Vector * sigmaVec = new VB_Vector(length);
522
sigmaVec->setAll(inputSigma);
524
return curvefit(xVec, yVec, sigmaVec, var3min, var1_init, var2_init, var3_init, outputFile);
527
/******************************************************************************************
528
* Wrapper function for fitting calculation. It accepts filenames for x, y, sigma values. *
529
******************************************************************************************/
530
VB_Vector * curvefit(const char *xFilename, const char *yFilename, const char *sigmaFilename,
531
double var3min, double var1_init, double var2_init, double var3_init,
532
const char *outputFile)
534
// Check xFilename format
535
VB_Vector * xVec = new VB_Vector();
536
string xString(xFilename);
537
if (xVec->ReadFile(xString)) {
538
printf("Invalid file format for X: %s\n", xFilename);
542
// Check yFilename format
543
VB_Vector * yVec = new VB_Vector();
544
string yString(yFilename);
545
if (yVec->ReadFile(yString)) {
546
printf("Invalid file format for Y: %s\n", yFilename);
550
// Check sigmaFilename format
551
VB_Vector * sigmaVec = new VB_Vector();
552
string sigmaString(sigmaFilename);
553
if (sigmaVec->ReadFile(sigmaString)) {
554
printf("Invalid file format for sigma: %s\n", sigmaFilename);
558
return curvefit(xVec, yVec, sigmaVec, var3min, var1_init, var2_init, var3_init, outputFile);
561
/************************************************************************************************
562
* This is the generic function to get fitting for a input power spectrum vector.
563
* It has the following arguments:
564
* #1 (psVecIn): input vector which include power spectrum elements
565
* #2 (ignorePSIn): another vector which defines frequencies ignored (elements which is equal to 0)
566
* #3: input TR value in unit of ms (default is 2000)
567
* #4: input sigma value (default is 0.1)
568
* #5: initial guess for parameter 1 (default is 20.0)
569
* #6: initial guess for parameter 2 (default is 2.0)
570
* #7: initial guess for parameter 3 (default is -0.0001)
571
* #8: output filename (default is empty)
572
* #9: flag to show printout or not (default is yes)
573
************************************************************************************************/
574
VB_Vector * fitOneOverF(VB_Vector *psVecIn, VB_Vector *ignorePSin, double var3min, double TRin,
575
double sigma, double var1, double var2, double var3, const char * outputFile, bool printFlag)
577
int cmpStatus = cmpVec(psVecIn, ignorePSin);
578
if (cmpStatus == 2) {
579
printf("Error: The number of elements in ps vector doesn't match ignorePS vector.\n");
583
double TR = TRin / 1000.0;
584
double var1init = var1, var2init = var2, var3init = var3;
586
int numObs = psVecIn->getLength();
587
double duration = numObs * TR;
588
int sigLen = numObs / 2 - 4;
589
VB_Vector *xVec = new VB_Vector(sigLen);
590
VB_Vector *signal = new VB_Vector(sigLen);
591
VB_Vector *ignoreVec = new VB_Vector(sigLen);
593
for (int i = 0; i < sigLen; i++) {
594
double tmpVal = ignorePSin->getElement(i + 1);
595
ignoreVec->setElement(i, tmpVal);
596
xVec->setElement(i, ((double)i + 1.0)/ duration);
597
tmpVal = sqrt(psVecIn->getElement(i + 1));
598
signal->setElement(i, tmpVal);
600
var2init = signal->getElement(sigLen - 1);
603
size_t numNonZero = findNonZero(ignoreVec);
605
if (numNonZero == 0) {
606
printf("Error: no freqency found in ignorePS vector\n");
610
VB_Vector *subx = new VB_Vector(numNonZero);
611
VB_Vector *subsignal = new VB_Vector(numNonZero);
612
VB_Vector *sigmaVec = new VB_Vector(numNonZero);
613
sigmaVec->setAll(sigma);
615
for (size_t i = 0; i < signal->getLength(); i++) {
616
if (ignoreVec->getElement(i) != 0) {
617
subsignal->setElement(counter, signal->getElement(i));
618
subx->setElement(counter, xVec->getElement(i));
626
printf("Error: no freqency found in ignorePS vector\n");
630
VB_Vector * newVec = new VB_Vector(8);
631
newVec = curvefit(subx, subsignal, sigmaVec, var3min, var1init, var2init, var3init, outputFile, printFlag);
637
/******************************************************************************************
638
* Function used by generic fitOneOverF().
639
* Returns 0 if 1st vector has same number of elements as 2nd vector does;
640
* Returns 1 if 2nd vector's length is multiple of 1st vector's;
641
* Returns 2 otherwise
642
*******************************************************************************************/
643
int cmpVec(VB_Vector *shortVec, VB_Vector *longVec)
645
int length_s = shortVec->getLength();
646
int length_l = longVec->getLength();
648
if (length_s == length_l)
650
else if (length_l % length_s == 0)
656
/******************************************************************************************
657
* Function used by generic fitOneOverF().
658
* It returns number of non-zero elements in the input vector.
659
******************************************************************************************/
660
size_t findNonZero(VB_Vector *inputVec)
663
for (size_t i = 0; i < inputVec->getLength(); i++) {
664
if (inputVec->getElement(i) != 0)
673
/******************************************************************************************
674
* Wrapper function based on generic fitOneOverF().
675
* Only power spectrum vector is specified.
676
* All elements in ignorePS vector will be set to be 1 (no frequencies will be ignored.)
677
*******************************************************************************************/
678
VB_Vector * fitOneOverF(VB_Vector *psVec, double var3min, double TRin, double sigma, double var1,
679
double var2, double var3, const char * outputFile, bool printFlag)
681
VB_Vector *ignorePSin = new VB_Vector(psVec->getLength());
682
ignorePSin->setAll(1.0);
684
return fitOneOverF(psVec, ignorePSin, var3min, TRin, sigma, var1, var2, var3, outputFile, printFlag);
687
/******************************************************************************************
688
* Wrapper function based on generic fitOneOverF().
689
* It accepts a filename for ps vector.
690
*******************************************************************************************/
691
VB_Vector *fitOneOverF(const char *psFile, double var3min, double TRin, double sigma, double var1,
692
double var2, double var3, const char * outputFile, bool printFlag)
695
VB_Vector *psVec = new VB_Vector();
696
string psString(psFile);
697
if (psVec->ReadFile(psString)) {
698
printf("Invalid file format for power spectrum: %s\n", psFile);
702
VB_Vector *ignorePSin = new VB_Vector(psVec->getLength());
703
ignorePSin->setAll(1.0);
705
return fitOneOverF(psVec, ignorePSin, var3min, TRin, sigma, var1, var2, var3, outputFile, printFlag);
708
/******************************************************************************************
709
* Wrapper function based on generic fitOneOverF().
710
* It accepts a filename for reference function.
711
* Reference function will be used to build ignorePS vector
712
*******************************************************************************************/
713
VB_Vector * fitOneOverF(VB_Vector *psVec, const char *refFunc, double var3min, double TRin,
714
double sigma, double var1, double var2, double var3,
715
const char * outputFile, bool printFlag)
717
VB_Vector *refLocal = new VB_Vector();
719
int refStat = getCondVec(refFunc, condKey, refLocal);
721
printf("File not readable: %s\n", refFunc);
724
else if (refStat == -2) {
725
printf("Error: different number of keys in header and content: %s\n", refFunc);
728
else if (refStat == 1) {
729
printf("Error: different keys in header and content: %s\n", refFunc);
733
refLocal->meanCenter();
734
VB_Vector *psRef = new VB_Vector(refLocal->getLength());
735
refLocal->getPS(psRef);
737
VB_Vector *ignorePS = new VB_Vector(refLocal->getLength());
738
ignorePS->setAll(1.0);
739
double refPSmax = psRef->getMaxElement();
741
for (size_t i = 0; i < ignorePS->getLength(); i++) {
742
if (psRef->getElement(i) > refPSmax * 0.01)
743
ignorePS->setElement(i, 0);
747
return fitOneOverF(psVec, ignorePS, var3min, TRin, sigma, var1, var2, var3, outputFile, printFlag);
750
/******************************************************************************************
751
* Wrapper function based on generic fitOneOverF().
752
* It accepts a filename for reference function.
753
* Reference function will be used to build ignorePS vector
754
*******************************************************************************************/
755
VB_Vector * fitOneOverF(const char *psFile, const char *refFunc, double var3min, double TRin,
756
double sigma, double var1, double var2, double var3,
757
const char * outputFile, bool printFlag)
759
VB_Vector *psVec = new VB_Vector();
760
string psString(psFile);
761
if (psVec->ReadFile(psString)) {
762
printf("Invalid file format for power spectrum: %s\n", psFile);
766
VB_Vector *refLocal = new VB_Vector();
768
int refStat = getCondVec(refFunc, condKey, refLocal);
770
printf("File not readable: %s\n", refFunc);
773
else if (refStat == -2) {
774
printf("Error: different number of keys in header and content: %s\n", refFunc);
777
else if (refStat == 1) {
778
printf("Error: different keys in header and content: %s\n", refFunc);
782
refLocal->meanCenter();
783
VB_Vector *psRef = new VB_Vector(refLocal->getLength());
784
refLocal->getPS(psRef);
786
VB_Vector *ignorePS = new VB_Vector(refLocal->getLength());
787
ignorePS->setAll(1.0);
788
double refPSmax = psRef->getMaxElement();
790
for (size_t i = 0; i < ignorePS->getLength(); i++) {
791
if (psRef->getElement(i) > refPSmax * 0.01)
792
ignorePS->setElement(i, 0);
796
return fitOneOverF(psVec, ignorePS, var3min, TRin, sigma, var1, var2, var3, outputFile, printFlag);
799
/******************************************************************************************
800
* generic function to make time domain representation of the 1/f curve given a set of
801
* three parameters, number of images and TR (default is 2000 ms)
802
******************************************************************************************/
803
VB_Vector * makeOneOverF(int numImages, double var1, double var2, double var3, double TRin)
805
double TR = TRin / 1000.0;
806
int halfLen = numImages / 2;
809
VB_Vector *freqVec = new VB_Vector(numImages);
810
for (int i = 1; i <= halfLen; i++) {
811
xVal = (double) i / (TR * numImages);
812
yVal = 1.0 / (var1 * (xVal + var3)) + var2;
813
freqVec->setElement(i, yVal);
816
if (numImages % 2 == 0) {
818
for (int i = halfLen + 1; i < numImages; i++) {
819
yVal = freqVec->getElement(halfLen - j);
820
freqVec->setElement(i, yVal);
826
for (int i = halfLen + 2; i < numImages; i++) {
827
yVal = freqVec->getElement(halfLen - j);
828
freqVec->setElement(i, yVal);
831
yVal = freqVec->getElement(halfLen + 2);
832
freqVec->setElement(halfLen + 1, yVal);
835
freqVec->setElement(0, 0);
837
freqVec->scaleInPlace(1.0 / freqVec->getMaxElement());
838
VB_Vector *realPart = new VB_Vector(numImages);
839
VB_Vector *imagPart = new VB_Vector(numImages);
841
freqVec->ifft(realPart, imagPart);
847
/**************************************************************************
848
* wrapper function based on generic makeOneOverF()
849
* It accepts a filename which includes three fitting parameters.
850
**************************************************************************/
851
VB_Vector * makeOneOverF(int numImages, const char * paramFile, double TRin)
853
VB_Vector *initGuess = new VB_Vector(paramFile);
854
double var1= initGuess->getElement(0);
855
double var2= initGuess->getElement(1);
856
double var3= initGuess->getElement(2);
858
return makeOneOverF(numImages, var1, var2, var3, TRin);
added
removed
utils/fitOneOverF.cpp
