1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 /*---------------------------------------------------------------------------*/
8 /*
9 Plugin to perform wavelet analysis of time series data.
10
11 File: plug_wavelet.c
12 Author: B. Douglas Ward
13 Date: 28 March 2000
14
15 Mod: Set MAX_NAME_LENGTH equal to THD_MAX_NAME.
16 Date: 02 December 2002
17 */
18
19
20 /*---------------------------------------------------------------------------*/
21 /*
22 Software identification.
23 */
24
25 #define PROGRAM_NAME "plug_wavelets" /* name of this program */
26 #define PROGRAM_AUTHOR "B. Douglas Ward" /* program author */
27 #define PROGRAM_INITIAL "28 March 2000" /* initial program release date */
28 #define PROGRAM_LATEST "02 December 2002" /* latest program revision date*/
29
30 /*---------------------------------------------------------------------------*/
31 /*
32 Include header files and source code.
33 */
34
35 #include "afni.h"
36 #include "Wavelets.h"
37 #include "Wavelets.c"
38
39 /*---------------------------------------------------------------------------*/
40 /*
41 Global constants.
42 */
43
44 #define MAX_NAME_LENGTH THD_MAX_NAME /* max. string length for file names */
45 #define MAX_FILTERS 20 /* max. number of blocking filters */
46 #define MAX_BAND 20 /* max. frequency band */
47
48 /*---------------------------------------------------------------------------*/
49
50 #ifndef ALLOW_PLUGINS
51 # error "Plugins not properly set up -- see machdep.h"
52 #endif
53
54
55 /*------------- string to 'help' the user -------------*/
56
57 static char helpstring[] =
58 " Purpose: Wavelet Analysis of FMRI time series data. \n"
59 " \n"
60 " Control: Wavelet = Type of wavelet to be used in analysis \n"
61 " NFirst = Number of first time series point to use. \n"
62 " NLast = Number of last time series point to use. \n"
63 " \n"
64 " Filter: Type = Type of filtering to apply in this \n"
65 " time-frequency window: \n"
66 " Stop = set wavelet coefficients to zero \n"
67 " Baseline = assign wavelet coefficients to baseline model \n"
68 " Signal = assign wavelet coefficients to signal model \n"
69 " \n"
70 " Band = Frequency band at which to apply filtering. \n"
71 " Min TR = Minimum value for time window (in TR) \n"
72 " Max TR = Maximum value for time window (in TR) \n"
73 " \n"
74 "For more information, see 'Wavelet Analysis of FMRI Time Series Data' \n"
75 "which is contained in file 3dWavelets.ps of the AFNI distribution. \n"
76 ;
77
78
79 /*------- Strings for filter type ------*/
80 #define NFILTER 3
81 static char * filter_strings[NFILTER] = {"Stop", "Baseline", "Signal"};
82
83
84 /*--------------- prototypes for internal routines ---------------*/
85
86 char * WA_main( PLUGIN_interface * ) ; /* the entry point */
87
88 void WA_fwt (int nt, double to, double dt, float * vec, char ** label);
89 void WA_fit (int nt, double to, double dt, float * vec, char ** label);
90 void WA_sgnl (int nt, double to, double dt, float * vec, char ** label);
91 void WA_err (int nt, double to, double dt, float * vec, char ** label);
92
93
94
95 /*---------------- global data -------------------*/
96
97 static PLUGIN_interface * global_plint = NULL;
98
99 static int plug_wavelet_type=0; /* type of wavelet to use in analysis */
100 static int plug_NFirst=0; /* first image from input 3d+time dataset to use */
101 static int plug_NLast=32767; /* last image from input 3d+time dataset to use */
102 static int plug_initialize=0; /* flag for perform initialization */
103 static int plug_prev_nt=0; /* previous time series length */
104 static int plug_filter_type=0; /* type of filter to use in analysis */
105
106 static int plug_num_stop_filters = 0; /* number of stop filters */
107 static int plug_stop_band[MAX_FILTERS]; /* freq band for stop filter */
108 static int plug_stop_mintr[MAX_FILTERS]; /* min time for stop filter */
109 static int plug_stop_maxtr[MAX_FILTERS]; /* max time for stop filter */
110 static float * plug_stop_filter = NULL; /* select wavelet coefs. to stop */
111
112 static int plug_num_base_filters = 0; /* number of base filters */
113 static int plug_base_band[MAX_FILTERS]; /* freq band for base filter */
114 static int plug_base_mintr[MAX_FILTERS]; /* min time for base filter */
115 static int plug_base_maxtr[MAX_FILTERS]; /* max time for base filter */
116 static float * plug_base_filter = NULL; /* select wavelet coefs.
117 for baseline */
118
119 static int plug_num_sgnl_filters = 0; /* number of signal filters */
120 static int plug_sgnl_band[MAX_FILTERS]; /* freq band for signal filter */
121 static int plug_sgnl_mintr[MAX_FILTERS]; /* min time for signal filter */
122 static int plug_sgnl_maxtr[MAX_FILTERS]; /* max time for signal filter */
123 static float * plug_sgnl_filter = NULL; /* select wavelet coefs. for signal */
124
125
126 /*---------------------------------------------------------------------------*/
127 /*
128 Set up the interface to the user.
129 */
130
131
132 DEFINE_PLUGIN_PROTOTYPE
133
PLUGIN_init(int ncall)134 PLUGIN_interface * PLUGIN_init( int ncall )
135 {
136 PLUGIN_interface * plint ; /* will be the output of this routine */
137 int is; /* filter index */
138
139 CHECK_IF_ALLOWED("WAVELETS","Wavelets") ; /* 30 Sep 2016 */
140
141 if( ncall > 0 ) return NULL ; /* generate interface for ncall 0 */
142
143
144 /***** do interface #0 *****/
145
146 /*---------------- set titles and call point ----------------*/
147
148 plint = PLUTO_new_interface ("Wavelets" ,
149 "Wavelet Analysis of Time Series Data" ,
150 helpstring, PLUGIN_CALL_VIA_MENU, WA_main);
151
152 global_plint = plint ; /* make global copy */
153
154 PLUTO_add_hint (plint,
155 "Control Wavelet Analysis Functions");
156
157 PLUTO_set_sequence( plint , "A:funcs:fitting" ) ;
158
159
160 /*----- Initialize Global Variables -----*/
161 for (is =0; is < MAX_FILTERS; is++)
162 {
163 plug_stop_band[is] = 0;
164 plug_stop_mintr[is] = 0.0;
165 plug_stop_maxtr[is] = 0.0;
166 plug_base_band[is] = 0;
167 plug_base_mintr[is] = 0.0;
168 plug_base_maxtr[is] = 0.0;
169 plug_sgnl_band[is] = 0;
170 plug_sgnl_mintr[is] = 0.0;
171 plug_sgnl_maxtr[is] = 0.0;
172 }
173
174
175 /*----- Parameters -----*/
176 PLUTO_add_option (plint, "Control", "Control", TRUE);
177 PLUTO_add_string (plint, "Wavelet", MAX_WAVELET_TYPE, WAVELET_TYPE_name,
178 plug_wavelet_type);
179 PLUTO_add_number (plint, "NFirst", 0, 32767, 0, plug_NFirst, TRUE);
180 PLUTO_add_number (plint, "NLast", 0, 32767, 0, plug_NLast, TRUE);
181
182
183 /*----- Stop Filters -----*/
184 for (is = 0; is < MAX_FILTERS; is++)
185 {
186 PLUTO_add_option (plint, "Filter", "Filter", FALSE);
187 PLUTO_add_string (plint, "Type", NFILTER, filter_strings,
188 plug_filter_type);
189 PLUTO_add_number (plint, "Band", -1, MAX_BAND, 0, 0, TRUE);
190 PLUTO_add_number (plint, "Min TR", 0, 10000, 0, 0, TRUE);
191 PLUTO_add_number (plint, "Max TR", 0, 10000, 0, 0, TRUE);
192 }
193
194
195 /*--------- done with interface setup ---------*/
196 PLUTO_register_1D_funcstr ("WA_FWT", WA_fwt);
197 PLUTO_register_1D_funcstr ("WA_Fit", WA_fit);
198 PLUTO_register_1D_funcstr ("WA_Sgnl", WA_sgnl);
199 PLUTO_register_1D_funcstr ("WA_Err", WA_err);
200
201
202 return plint ;
203 }
204
205
206 /*---------------------------------------------------------------------------*/
207 /*
208 Main routine for this plugin (will be called from AFNI).
209 If the return string is not NULL, some error transpired, and
210 AFNI will popup the return string in a message box.
211 */
212
WA_main(PLUGIN_interface * plint)213 char * WA_main( PLUGIN_interface * plint )
214 {
215 char * str; /* input string */
216 int is; /* filter index */
217
218
219 /*----- reset flag for successful initialization -----*/
220 plug_initialize = 0;
221
222
223 /*--------- go to Control input line ---------*/
224 PLUTO_next_option (plint);
225 str = PLUTO_get_string (plint);
226 plug_wavelet_type = PLUTO_string_index (str, MAX_WAVELET_TYPE,
227 WAVELET_TYPE_name);
228 plug_NFirst = PLUTO_get_number (plint);
229 plug_NLast = PLUTO_get_number (plint);
230
231
232 /*------ read input line(s) -----*/
233 plug_num_stop_filters = 0;
234 plug_num_base_filters = 0;
235 plug_num_sgnl_filters = 0;
236
237 do
238 {
239 str = PLUTO_get_optiontag(plint);
240 if (str == NULL) break;
241 if (strcmp (str, "Filter") != 0)
242 return "************************\n"
243 "Illegal optiontag found!\n"
244 "************************";
245
246
247 /*----- Read Filter Specification Line -----*/
248 str = PLUTO_get_string (plint);
249 plug_filter_type = PLUTO_string_index (str, NFILTER, filter_strings);
250
251 switch (plug_filter_type)
252 {
253 case 0:
254 {
255 plug_stop_band[plug_num_stop_filters] = PLUTO_get_number(plint);
256 plug_stop_mintr[plug_num_stop_filters]
257 = PLUTO_get_number(plint);
258 plug_stop_maxtr[plug_num_stop_filters]
259 = PLUTO_get_number(plint);
260
261 if (plug_stop_mintr[plug_num_stop_filters] >
262 plug_stop_maxtr[plug_num_stop_filters])
263 return "*************************\n"
264 "Require Min TR <= Max TR \n"
265 "*************************" ;
266
267 plug_num_stop_filters++;
268 break;
269 }
270
271 case 1:
272 {
273 plug_base_band[plug_num_base_filters] = PLUTO_get_number(plint);
274 plug_base_mintr[plug_num_base_filters]
275 = PLUTO_get_number(plint);
276 plug_base_maxtr[plug_num_base_filters]
277 = PLUTO_get_number(plint);
278
279 if (plug_base_mintr[plug_num_base_filters] >
280 plug_base_maxtr[plug_num_base_filters])
281 return "*************************\n"
282 "Require Min TR <= Max TR \n"
283 "*************************" ;
284
285 plug_num_base_filters++;
286 break;
287 }
288
289 case 2:
290 {
291 plug_sgnl_band[plug_num_sgnl_filters]=PLUTO_get_number(plint);
292 plug_sgnl_mintr[plug_num_sgnl_filters]
293 = PLUTO_get_number(plint);
294 plug_sgnl_maxtr[plug_num_sgnl_filters]
295 = PLUTO_get_number(plint);
296
297 if (plug_sgnl_mintr[plug_num_sgnl_filters] >
298 plug_sgnl_maxtr[plug_num_sgnl_filters])
299 return "*************************\n"
300 "Require Min TR <= Max TR \n"
301 "*************************" ;
302
303 plug_num_sgnl_filters++;
304 break;
305 }
306
307 }
308 }
309 while (1);
310
311
312 /*----- Identify software -----*/
313 printf ("\n\n");
314 printf ("Program: %s \n", PROGRAM_NAME);
315 printf ("Author: %s \n", PROGRAM_AUTHOR);
316 printf ("Initial Release: %s \n", PROGRAM_INITIAL);
317 printf ("Latest Revision: %s \n", PROGRAM_LATEST);
318 printf ("\n");
319
320 /*----- show current input options -----*/
321 printf ("\nControls: \n");
322 printf ("Wavelet Type = %10s \n", WAVELET_TYPE_name[plug_wavelet_type]);
323 printf ("NFirst = %10d \n", plug_NFirst);
324 printf ("NLast = %10d \n", plug_NLast);
325
326 for (is = 0; is < plug_num_stop_filters; is++)
327 {
328 printf ("\nStop Filter: Band = %4d ", plug_stop_band[is]);
329 printf ("Min. TR = %4d Max. TR = %4d \n",
330 plug_stop_mintr[is], plug_stop_maxtr[is]);
331 }
332
333 for (is = 0; is < plug_num_base_filters; is++)
334 {
335 printf ("\nBaseline Filter: Band = %4d ", plug_base_band[is]);
336 printf ("Min. TR = %4d Max. TR = %4d \n",
337 plug_base_mintr[is], plug_base_maxtr[is]);
338 }
339
340 for (is = 0; is < plug_num_sgnl_filters; is++)
341 {
342 printf ("\nSignal Filter: Band = %4d ", plug_sgnl_band[is]);
343 printf ("Min. TR = %4d Max. TR = %4d \n",
344 plug_sgnl_mintr[is], plug_sgnl_maxtr[is]);
345 }
346
347
348 /*--- nothing left to do until data arrives ---*/
349 plug_initialize = 1 ; /* successful initialization */
350 plug_prev_nt = 0; /* previous time series length */
351
352 return NULL ;
353 }
354
355
356 /*---------------------------------------------------------------------------*/
357 /*
358 Perform the wavelet analysis.
359 */
360
calculate_results(int nt,float * vec,int * NFirst,int * NLast,char ** label,float ** coefts,float ** fitts,float ** sgnlts,float ** errts)361 int calculate_results
362 (
363 int nt, /* number of time points */
364 float * vec, /* input measured time series data */
365 int * NFirst, /* first image from input 3d+time dataset to use */
366 int * NLast, /* last image from input 3d+time dataset to use */
367 char ** label, /* string containing statistical summary of results */
368 float ** coefts, /* forward wavelet transform coefficients */
369 float ** fitts, /* full model fit to the time series data */
370 float ** sgnlts, /* signal model fit to the time series data */
371 float ** errts /* residual error time series */
372 )
373
374 {
375 float * ts_array; /* array of measured data for one voxel */
376 float * coef; /* regression parameters */
377 float sse_base; /* baseline model error sum of squares */
378 float sse_full; /* full model error sum of squares */
379 float ffull; /* full model F-statistic */
380 float rfull; /* full model R^2 stat. */
381
382 int N; /* number of data points */
383 int f; /* number of parameters removed by the filter */
384 int p; /* number of parameters in the full model */
385 int q; /* number of parameters in the baseline model */
386 int n; /* data point index */
387 int i; /* data point index */
388 int ok = 1;
389
390
391 /*----- Check initialization flag -----*/
392 if (! plug_initialize) return (0);
393
394
395 /*----- Initialize local variables -----*/
396 *NFirst = plug_NFirst;
397
398 *NLast = plug_NLast;
399 if (*NLast > nt-1) *NLast = nt-1;
400
401 N = *NLast - *NFirst + 1;
402 N = powerof2(my_log2(N));
403 *NLast = N + *NFirst - 1;
404
405
406 plug_stop_filter =
407 FWT_1d_stop_filter (plug_num_stop_filters, plug_stop_band,
408 plug_stop_mintr, plug_stop_maxtr, *NFirst, N);
409
410 plug_base_filter =
411 FWT_1d_pass_filter (plug_num_base_filters, plug_base_band,
412 plug_base_mintr, plug_base_maxtr, *NFirst, N);
413
414 plug_sgnl_filter =
415 FWT_1d_pass_filter (plug_num_sgnl_filters, plug_sgnl_band,
416 plug_sgnl_mintr, plug_sgnl_maxtr, *NFirst, N);
417
418 f = 0;
419 for (i = 0; i < N; i++)
420 if (plug_stop_filter[i] == 0.0)
421 {
422 f++;
423 plug_base_filter[i] = 0.0;
424 plug_sgnl_filter[i] = 0.0;
425 }
426
427 q = 0;
428 for (i = 0; i < N; i++)
429 if (plug_base_filter[i] == 1.0)
430 {
431 q++;
432 plug_sgnl_filter[i] = 1.0;
433 }
434
435 p = 0;
436 for (i = 0; i < N; i++)
437 if (plug_sgnl_filter[i] == 1.0)
438 {
439 p++;
440 }
441
442
443 /*----- Allocate memory for fitted time series and residuals -----*/
444 coef = (float *) malloc (sizeof(float) * p); MTEST (coef);
445 *coefts = (float *) malloc (sizeof(float) * N); MTEST (coefts);
446 *fitts = (float *) malloc (sizeof(float) * N); MTEST (fitts);
447 *sgnlts = (float *) malloc (sizeof(float) * N); MTEST (sgnlts);
448 *errts = (float *) malloc (sizeof(float) * N); MTEST (errts);
449
450
451 /*----- Extract data for this voxel -----*/
452 ts_array = vec + *NFirst;
453
454
455 /*----- Perform the wavelet analysis for this voxel-----*/
456 wavelet_analysis (plug_wavelet_type,
457 f, plug_stop_filter,
458 q, plug_base_filter,
459 p, plug_sgnl_filter,
460 N, ts_array, coef,
461 &sse_base, &sse_full, &ffull, &rfull,
462 *coefts, *fitts, *sgnlts, *errts);
463
464
465 /*----- Report results for this voxel -----*/
466 printf ("\nResults for Voxel: \n");
467 report_results (N, *NFirst, f, p, q,
468 plug_base_filter, plug_sgnl_filter,
469 coef, sse_base, sse_full, ffull, rfull,
470 label);
471 printf ("%s \n", *label);
472
473 plug_prev_nt = nt;
474
475
476 /*----- Release memory -----*/
477 free (plug_stop_filter); plug_stop_filter = NULL;
478 free (plug_base_filter); plug_base_filter = NULL;
479 free (plug_sgnl_filter); plug_sgnl_filter = NULL;
480 free (coef); coef = NULL;
481
482
483 return (ok);
484 }
485
486
487 /*---------------------------------------------------------------------------*/
488 /*
489 Calculate the forward wavelet transform coefficients.
490 */
491
WA_fwt(int nt,double to,double dt,float * vec,char ** label)492 void WA_fwt (int nt, double to, double dt, float * vec, char ** label)
493
494 {
495 int NFirst; /* first image from input 3d+time dataset to use */
496 int NLast; /* last image from input 3d+time dataset to use */
497 int n; /* time index */
498 int ok; /* RwcBoolean for successful calculation */
499 float * coefts = NULL; /* forward wavelet transform coefficients */
500 float * fitts = NULL; /* wavelet filtered time series */
501 float * sgnlts = NULL; /* signal model fitted time series */
502 float * errts = NULL; /* residual error time series */
503
504
505 /*----- Calculate the wavelet filtered time series data -----*/
506 ok = calculate_results (nt, vec, &NFirst, &NLast, label,
507 &coefts, &fitts, &sgnlts, &errts);
508 if (!ok)
509 {
510 for (n = 0; n < nt; n++) vec[n] = 0.0;
511 return;
512 }
513
514
515 /*----- Store the filtered time series data back into the array -----*/
516 for (n = NFirst; n <= NLast; n++)
517 {
518 vec[n] = coefts[n-NFirst];
519 }
520
521 for (n = 0; n < NFirst; n++)
522 vec[n] = 0.0;
523
524 for (n = NLast+1; n < nt; n++)
525 vec[n] = 0.0;
526
527
528 /*----- Deallocate memory -----*/
529 free (coefts); coefts = NULL;
530 free (fitts); fitts = NULL;
531 free (sgnlts); sgnlts = NULL;
532 free (errts); errts = NULL;
533
534
535 return;
536 }
537
538
539 /*---------------------------------------------------------------------------*/
540 /*
541 Calculate the full model wavelet fit to the time series data.
542 */
543
WA_fit(int nt,double to,double dt,float * vec,char ** label)544 void WA_fit (int nt, double to, double dt, float * vec, char ** label)
545
546 {
547 int NFirst; /* first image from input 3d+time dataset to use */
548 int NLast; /* last image from input 3d+time dataset to use */
549 int n; /* time index */
550 int ok; /* RwcBoolean for successful calculation */
551 float * coefts = NULL; /* forward wavelet transform coefficients */
552 float * fitts = NULL; /* wavelet filtered time series */
553 float * sgnlts = NULL; /* signal model fitted time series */
554 float * errts = NULL; /* residual error time series */
555
556
557 /*----- Calculate the wavelet filtered time series data -----*/
558 ok = calculate_results (nt, vec, &NFirst, &NLast, label,
559 &coefts, &fitts, &sgnlts, &errts);
560 if (!ok)
561 {
562 for (n = 0; n < nt; n++) vec[n] = 0.0;
563 return;
564 }
565
566
567 /*----- Store the filtered time series data back into the array -----*/
568 for (n = NFirst; n <= NLast; n++)
569 {
570 vec[n] = fitts[n-NFirst];
571 }
572
573 for (n = 0; n < NFirst; n++)
574 vec[n] = vec[NFirst];
575
576 for (n = NLast+1; n < nt; n++)
577 vec[n] = vec[NLast];
578
579
580 /*----- Deallocate memory -----*/
581 free (coefts); coefts = NULL;
582 free (fitts); fitts = NULL;
583 free (sgnlts); sgnlts = NULL;
584 free (errts); errts = NULL;
585
586
587 return;
588 }
589
590
591 /*---------------------------------------------------------------------------*/
592 /*
593 Calculate the signal model wavelet fit to the time series data.
594 */
595
WA_sgnl(int nt,double to,double dt,float * vec,char ** label)596 void WA_sgnl (int nt, double to, double dt, float * vec, char ** label)
597
598 {
599 int NFirst; /* first image from input 3d+time dataset to use */
600 int NLast; /* last image from input 3d+time dataset to use */
601 int n; /* time index */
602 int ok; /* RwcBoolean for successful calculation */
603 float * coefts = NULL; /* forward wavelet transform coefficients */
604 float * fitts = NULL; /* wavelet filtered time series */
605 float * sgnlts = NULL; /* signal model fitted time series */
606 float * errts = NULL; /* residual error time series */
607
608
609 /*----- Calculate the wavelet filtered time series data -----*/
610 ok = calculate_results (nt, vec, &NFirst, &NLast, label,
611 &coefts, &fitts, &sgnlts, &errts);
612 if (!ok)
613 {
614 for (n = 0; n < nt; n++) vec[n] = 0.0;
615 return;
616 }
617
618
619 /*----- Store the filtered time series residuals back into the array -----*/
620 for (n = NFirst; n <= NLast; n++)
621 {
622 vec[n] = sgnlts[n-NFirst];
623 }
624
625 for (n = 0; n < NFirst; n++)
626 vec[n] = 0.0;
627
628 for (n = NLast+1; n < nt; n++)
629 vec[n] = 0.0;
630
631
632 /*----- Deallocate memory -----*/
633 free (coefts); coefts = NULL;
634 free (fitts); fitts = NULL;
635 free (sgnlts); sgnlts = NULL;
636 free (errts); errts = NULL;
637
638
639 return;
640 }
641
642
643 /*---------------------------------------------------------------------------*/
644 /*
645 Calculate the residual error from the full model wavelet fit.
646 */
647
WA_err(int nt,double to,double dt,float * vec,char ** label)648 void WA_err (int nt, double to, double dt, float * vec, char ** label)
649
650 {
651 int NFirst; /* first image from input 3d+time dataset to use */
652 int NLast; /* last image from input 3d+time dataset to use */
653 int n; /* time index */
654 int ok; /* RwcBoolean for successful calculation */
655 float * coefts = NULL; /* forward wavelet transform coefficients */
656 float * fitts = NULL; /* wavelet filtered time series */
657 float * sgnlts = NULL; /* signal model fitted time series */
658 float * errts = NULL; /* residual error time series */
659
660
661 /*----- Calculate the wavelet filtered time series data -----*/
662 ok = calculate_results (nt, vec, &NFirst, &NLast, label,
663 &coefts, &fitts, &sgnlts, &errts);
664 if (!ok)
665 {
666 for (n = 0; n < nt; n++) vec[n] = 0.0;
667 return;
668 }
669
670
671 /*----- Store the filtered time series residuals back into the array -----*/
672 for (n = NFirst; n <= NLast; n++)
673 {
674 vec[n] = errts[n-NFirst];
675 }
676
677 for (n = 0; n < NFirst; n++)
678 vec[n] = 0.0;
679
680 for (n = NLast+1; n < nt; n++)
681 vec[n] = 0.0;
682
683
684 /*----- Deallocate memory -----*/
685 free (coefts); coefts = NULL;
686 free (fitts); fitts = NULL;
687 free (sgnlts); sgnlts = NULL;
688 free (errts); errts = NULL;
689
690
691 return;
692 }
693
694
695 /*---------------------------------------------------------------------------*/
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710