1 #include "mrilib.h"
2
3 /*** bandpass functions: 30 Apr 2009 -- RWCox ***/
4
5 /*--------------------------------------------------------------------------*/
6
7 static int bpwrn = 1 ;
8
9 static int nfft_fixed = 0 ;
10
THD_bandpass_set_nfft(int n)11 int THD_bandpass_set_nfft( int n )
12 {
13 nfft_fixed = (n >= 16) ? csfft_nextup_even(n) : 0 ;
14 return nfft_fixed ;
15 }
16
17 /*--------------------------------------------------------------------------*/
18 /*! Check THD_bandpass_vectors() input parameters for OK-ness.
19 Returns 1 if OK, 0 if not OK. If verb!=0, prints an info message.
20 *//*------------------------------------------------------------------------*/
21
THD_bandpass_OK(int nx,float dt,float fbot,float ftop,int verb)22 int THD_bandpass_OK( int nx , float dt , float fbot , float ftop , int verb )
23 {
24 int nfft,nhalf , jbot,jtop ; float df , qbot,qtop ;
25
26 if( ftop > ICOR_MAX_FTOP ) return 1 ; /* 26 Feb 2010 */
27
28 if( nx < 9 ) return 0 ;
29 if( dt <= 0.0f ) dt = 1.0f ;
30 if( fbot < 0.0f ) fbot = 0.0f ;
31 if( ftop <= fbot || ftop > ICOR_MAX_FTOP ) ftop = ICOR_MAX_FTOP ;
32 if( bpwrn && dt > 60.0f ){
33 WARNING_message("Your bandpass timestep (%f) is high.\n"
34 " Make sure units are 'sec', not 'msec'.\n"
35 " This warning will not be repeated." ,
36 dt);
37 bpwrn = 0;
38 }
39
40 nfft = (nfft_fixed >= nx) ? nfft_fixed : csfft_nextup_even(nx) ;
41 nhalf = nfft/2 - 1 ;
42 df = 1.0f / (nfft * dt) ; /* freq step */
43 qbot = rintf(fbot/df) ;
44 qtop = rintf(ftop/df) ;
45 jbot = ( qbot < nhalf ) ? (int)qbot : 0 ; /* band bot index */
46 jtop = ( qtop < nhalf ) ? (int)qtop : nhalf ; /* band top index */
47 if( jbot+1 >= jtop ){
48 WARNING_message(
49 "bandpass: fbot=%g and ftop=%g too close ==> jbot=%d jtop=%d [nfft=%d dt=%g]",
50 fbot,ftop,jbot,jtop,nfft,dt) ;
51 WARNING_message( " -- ignoring bandpass filter" ) ;
52 jbot = 0 ; jtop = nhalf ;
53 }
54
55 if( verb )
56 ININFO_message(
57 "bandpass: ntime=%d nFFT=%d dt=%.6g dFreq=%.6g Nyquist=%.6g passband indexes=%d..%d",
58 nx, nfft, dt, df, (nfft/2)*df, jbot, jtop) ;
59
60 return 1 ;
61 }
62
63 /*--------------------------------------------------------------------------*/
64 /*! Return the number of degrees of freedom that would remain after
65 bandpassing (the dimension of the subspace the data will be
66 projected into).
67
68 Returns twice the length of the passband index range.
69 Returns 0 if life is bad.
70
71 based on THD_bandpass_OK 18 Mar 2015 [rickr]
72 *//*------------------------------------------------------------------------*/
73
THD_bandpass_remain_dim(int nx,float dt,float fbot,float ftop,int verb)74 int THD_bandpass_remain_dim(int nx, float dt, float fbot, float ftop, int verb)
75 {
76 int nfft,nhalf , jbot,jtop , dim ; float df , qbot,qtop ;
77
78 if( nx < 9 ) {
79 if( verb ) WARNING_message("length %d too short for bandpassing", nx);
80 return 0 ; }
81
82 if( dt <= 0.0f ) dt = 1.0f ;
83 if( fbot < 0.0f ) fbot = 0.0f ;
84 if( ftop <= fbot || ftop > ICOR_MAX_FTOP ) ftop = ICOR_MAX_FTOP ;
85 if( verb && dt > 60.0f ){
86 WARNING_message("Your bandpass timestep (%f) is high.\n"
87 " Make sure units are 'sec', not 'msec'.\n"
88 " This warning will not be repeated." ,
89 dt);
90 }
91
92 nfft = (nfft_fixed >= nx) ? nfft_fixed : csfft_nextup_even(nx) ;
93 nhalf = nfft/2 - 1 ;
94 df = 1.0f / (nfft * dt) ; /* freq step */
95 qbot = rintf(fbot/df) ;
96 qtop = rintf(ftop/df) ;
97 jbot = ( qbot < nhalf ) ? (int)qbot : 0 ; /* band bot index */
98 jtop = ( qtop < nhalf ) ? (int)qtop : nhalf ; /* band top index */
99 if( jbot+1 >= jtop ){
100 WARNING_message(
101 "bandpass: fbot=%g and ftop=%g too close ==> jbot=%d jtop=%d [nfft=%d dt=%g]",
102 fbot,ftop,jbot,jtop,nfft,dt) ;
103 WARNING_message( " -- ignoring bandpass filter" ) ;
104 jbot = 0 ; jtop = nhalf ;
105 }
106
107 dim = 2*(jtop-jbot+1); if( dim > nx ) dim = nx ;
108 return dim ;
109 }
110
111 /*--------------------------------------------------------------------------*/
112 /*! Bandpass a set of vectors, optionally removing some orts as well.
113 - Uses FFTs for the bandpass-ization, and least squares for the ort-ing.
114 - To do a highpass only, set ftop to something larger than the Nyquist
115 frequency (e.g., 999999.9). To do a lowpass only, set fbot to 0.0.
116 - However, the 0 and Nyquist frequencies are always removed.
117 - Return value is the number of linear dimensions projected out.
118 If 0 is returned, something funky happened.
119 *//*------------------------------------------------------------------------*/
120
THD_bandpass_vectors(int nlen,int nvec,float ** vec,float dt,float fbot,float ftop,int qdet,int nort,float ** ort)121 int THD_bandpass_vectors( int nlen , int nvec , float **vec ,
122 float dt , float fbot , float ftop ,
123 int qdet , int nort , float **ort )
124 {
125 int nfft,nby2,nhalf , iv, jbot,jtop , ndof=0 ; register int jj ;
126 float df , tapr , qbot,qtop ;
127 register float *xar, *yar=NULL ;
128 register complex *zar ; complex Zero={0.0f,0.0f} ;
129
130 ENTRY("THD_bandpass_vectors") ;
131
132 if( (ftop <= fbot || ftop > ICOR_MAX_FTOP) && qdet < 0 && (nort <= 0 || ort == NULL) )
133 RETURN(ndof) ; /* 26 Feb 2010: do nothing at all? */
134
135 if( nlen < 9 || nvec < 1 || vec == NULL ){ /* should not happen */
136 ERROR_message("bad filtering/bandpass data?"); RETURN(ndof);
137 }
138
139 if( nort >= nlen ){
140 ERROR_message("Too many bandpass orts?") ; RETURN(ndof);
141 }
142
143 if( bpwrn && dt > 60.0f ) {
144 WARNING_message("Your bandpass timestep (%f) is high.\n"
145 " Make sure units are 'sec', not 'msec'.\n"
146 " This warning will not be repeated." ,
147 dt);
148 bpwrn = 0;
149 }
150
151 if( dt <= 0.0f ) dt = 1.0f ;
152 if( fbot < 0.0f ) fbot = 0.0f ;
153 if( ftop <= fbot ) fbot = ICOR_MAX_FTOP ;
154
155 /** setup for FFT **/
156
157 nfft = (nfft_fixed >= nlen) ? nfft_fixed : csfft_nextup_even(nlen) ;
158 nby2 = nfft/2 ; nhalf = nby2-1 ;
159
160 df = 1.0f / (nfft * dt) ; /* frequency resolution */
161 qbot = rintf(fbot/df) ;
162 qtop = rintf(ftop/df) ;
163 jbot = ( qbot < nhalf ) ? (int)qbot : 0 ; /* band bot index */
164 jtop = ( qtop < nhalf ) ? (int)qtop : nhalf ; /* band top index */
165
166 if( ftop >= ICOR_MAX_FTOP ) jtop = nhalf ;
167
168 if( jtop >= nby2 ) jtop = nhalf ; /* can't go past Nyquist! */
169 if( jbot >= jtop+1 ){
170 jbot = 0 ; jtop = nhalf ;
171 }
172
173 /** quadratic detrending first? (should normally be used) **/
174
175 switch( qdet ){
176 case 2:
177 ndof += 2 ;
178 ININFO_message(" -- quadratic detrend") ;
179 for( iv=0 ; iv < nvec ; iv++ )
180 THD_quadratic_detrend( nlen, vec[iv], NULL,NULL,NULL ) ;
181 break ;
182
183 case 1:
184 ndof += 1 ;
185 ININFO_message(" -- linear detrend") ;
186 for( iv=0 ; iv < nvec ; iv++ )
187 THD_linear_detrend( nlen, vec[iv], NULL,NULL ) ;
188 break ;
189
190 case 0:
191 ININFO_message(" -- remove mean") ;
192 for( iv=0 ; iv < nvec ; iv++ )
193 THD_const_detrend( nlen, vec[iv], NULL ) ;
194 break ;
195 }
196
197 if( jbot > 0 || jtop < nhalf ){ /* 26 Feb 2010: do the FFTs */
198
199 zar = (complex *)malloc(sizeof(complex)*nfft) ; /* work array */
200 csfft_scale_inverse(1) ; /* scale inverse FFT by 1/nfft */
201
202 /** loop over vectors in pairs, FFT-ing and bandpassing **/
203
204 ndof += 2 ; /* for 0 and Nyquist freqs */
205
206 if( jbot >= 1 ) ndof += 2*jbot - 1 ; /* DOF for low freq */
207 ndof += 2*(nby2-jtop) - 1 ; /* DOF for high freq */
208
209 ININFO_message(" -- FFTs") ;
210 for( iv=0 ; iv < nvec ; iv+=2 ){
211
212 /* load a pair of vectors into zar to double up on FFTs of real data */
213
214 xar = vec[iv] ;
215 if( iv == nvec-1 ){ /* last one has nothing to pair with */
216 for( jj=0 ; jj < nlen ; jj++ ){ zar[jj].r = xar[jj] ; zar[jj].i = 0.0f ; }
217 } else {
218 yar = vec[iv+1] ;
219 for( jj=0 ; jj < nlen ; jj++ ){ zar[jj].r = xar[jj] ; zar[jj].i = yar[jj] ; }
220 }
221 for( jj=nlen ; jj < nfft ; jj++ ) zar[jj] = Zero ; /* zero fill */
222
223 csfft_cox( -1 , nfft , zar ) ; /*** the FFT ***/
224
225 /* delete unwanted frequencies */
226
227 zar[0] = zar[nby2] = Zero ;
228
229 tapr = (nort > 0 && ort != NULL) ? 0.05f : 0.5f ;
230
231 if( jbot >= 1 ){
232 zar[jbot].r *= tapr ; zar[jbot].i *= tapr ;
233 zar[nfft-jbot].r *= tapr ; zar[nfft-jbot].i *= tapr ;
234 for( jj=1 ; jj < jbot ; jj++ ) zar[jj] = zar[nfft-jj] = Zero ;
235 }
236
237 zar[jtop].r *= tapr ; zar[jtop].i *= tapr ;
238 zar[nfft-jtop].r *= tapr ; zar[nfft-jtop].i *= tapr ;
239 for( jj=jtop+1 ; jj < nby2 ; jj++ ) zar[jj] = zar[nfft-jj] = Zero ;
240
241 csfft_cox( 1 , nfft , zar ) ; /*** inverse FFT ***/
242
243 /* unload vector pair back into original data arrays */
244
245 if( iv == nvec-1 ){
246 for( jj=0 ; jj < nlen ; jj++ ) xar[jj] = zar[jj].r ;
247 } else {
248 for( jj=0 ; jj < nlen ; jj++ ){ xar[jj] = zar[jj].r ; yar[jj] = zar[jj].i ; }
249 }
250
251 } /* end of loop over vector pairs: bandpassing now done */
252
253 free(zar) ; csfft_scale_inverse(0) ;
254
255 } /* end of FFTs */
256
257 /*** remove orts? ***/
258
259 if( nort > 0 && ort != NULL ){
260 float **qort = (float **)malloc(sizeof(float *)*nort) ;
261 MRI_IMAGE *qim , *pim ; float *par, *qar , *rar , *pt,*qt ;
262 register float sum , xt ; register int kk ;
263
264 /* must bandpass copy of orts first -- via recursion
265 (so we don't re-introduce any of the removed frequencies) */
266
267 qim = mri_new( nlen , nort , MRI_float ) ; /** [Q] = nlen X nort **/
268 qar = MRI_FLOAT_PTR(qim) ; /** load with the orts **/
269 for( iv=0 ; iv < nort ; iv++ ){
270 qort[iv] = qar + iv*nlen ;
271 memcpy( qort[iv] , ort[iv] , sizeof(float)*nlen ) ;
272 }
273
274 (void)THD_bandpass_vectors( nlen, nort, qort, dt, fbot, ftop, qdet, 0,NULL ) ;
275 free(qort) ;
276
277 /* compute pseudo-inverse ([P] = inv{[Q]'[Q]}[Q]') of bandpassed orts */
278
279 pim = mri_matrix_psinv( qim , NULL , 1.e-8 ) ; /** [P] = nort X nlen **/
280 if( pim == NULL ){ /* should not happen */
281 mri_free(qim) ;
282 ERROR_message("can't remove bandpass orts?") ;
283 RETURN(ndof) ;
284 }
285
286 par = MRI_FLOAT_PTR(pim) ; /* nort X nlen matrix */
287
288 /* Project the bandpassed orts out of data vectors:
289 for each vector [y], replace it with [y] - [Q][P][y] ;
290 this is more efficient than computing the nlen X nlen projection
291 matrix [I]-[Q][P] and then applying to each vector, since
292 that would require about nlen*nlen flops per vector, whereas
293 the 2 matrix multiplications by [P] then [Q] will require
294 about 2*nlen*nort flops per vector, a clear win if nort < nlen/2. */
295
296 rar = (float *)malloc(sizeof(float)*nort) ; /* will hold [P][y] */
297 for( iv=0 ; iv < nvec ; iv++ ){
298 xar = vec[iv] ; /* [y] vector */
299
300 /* compute nort-vector [r] = [P][y] */
301
302 for( jj=0 ; jj < nort ; jj++ ) rar[jj] = 0.0f ; /* initialize to 0 */
303 for( kk=0 ; kk < nlen ; kk++ ){
304 pt = par + kk*nort ; xt = xar[kk] ;
305 for( jj=0 ; jj < nort ; jj++ ) rar[jj] += pt[jj]*xt ;
306 }
307
308 /* now subtract [Q][r] from [y] */
309
310 for( jj=0 ; jj < nort ; jj++ ){
311 qt = qar + jj*nlen ; xt = rar[jj] ;
312 for( kk=0 ; kk < nlen ; kk++ ) xar[kk] -= qt[kk]*xt ;
313 }
314 }
315
316 free(rar) ; mri_free(pim) ; mri_free(qim) ;
317
318 ndof += nort ;
319
320 } /* de-ortification is done */
321
322 /** done **/
323
324 if( nfft > nlen ){
325 double fac = ((double)nlen)/(double)nfft ;
326 ndof = (int)rint(fac*ndof) ;
327 }
328 RETURN(ndof) ;
329 }
330
331 /*--------------------------------------------------------------------------*/
332 /*! Bandpass a vectim. */
333
THD_bandpass_vectim(MRI_vectim * mrv,float dt,float fbot,float ftop,int qdet,int nort,float ** ort)334 int THD_bandpass_vectim( MRI_vectim *mrv ,
335 float dt , float fbot , float ftop ,
336 int qdet , int nort , float **ort )
337 {
338 float **vec ; int nlen , nvec , ndof , kk ;
339
340 ENTRY("THD_bandpass_vectim") ;
341
342 if( mrv == NULL ) RETURN(0) ;
343
344 nvec = mrv->nvec ; nlen = mrv->nvals ;
345 vec = (float **)malloc(sizeof(float *)*nvec) ;
346 for( kk=0 ; kk < nvec ; kk++ ) vec[kk] = VECTIM_PTR(mrv,kk) ;
347
348 ndof = THD_bandpass_vectors( nlen , nvec , vec ,
349 dt , fbot , ftop , qdet , nort , ort ) ;
350
351 free(vec) ; RETURN(ndof) ;
352 }
353
354 /*------------------- 08 Oct 2010: functions for despiking ----------------*/
355
356 #undef SWAP
357 #define SWAP(x,y) (temp=x,x=y,y=temp)
358
359 #undef SORT2
360 #define SORT2(a,b) if(a>b) SWAP(a,b)
361
362 /*--- fast median of 9 values ---*/
363
median9f(float * p)364 static INLINE float median9f(float *p)
365 {
366 register float temp ;
367 SORT2(p[1],p[2]) ; SORT2(p[4],p[5]) ; SORT2(p[7],p[8]) ;
368 SORT2(p[0],p[1]) ; SORT2(p[3],p[4]) ; SORT2(p[6],p[7]) ;
369 SORT2(p[1],p[2]) ; SORT2(p[4],p[5]) ; SORT2(p[7],p[8]) ;
370 SORT2(p[0],p[3]) ; SORT2(p[5],p[8]) ; SORT2(p[4],p[7]) ;
371 SORT2(p[3],p[6]) ; SORT2(p[1],p[4]) ; SORT2(p[2],p[5]) ;
372 SORT2(p[4],p[7]) ; SORT2(p[4],p[2]) ; SORT2(p[6],p[4]) ;
373 SORT2(p[4],p[2]) ; return(p[4]) ;
374 }
375 #undef SORT2
376 #undef SWAP
377
378 /*--- get the local median and MAD of values vec[j-4 .. j+4] ---*/
379
380 #undef mead9
381 #define mead9(j) \
382 { float qqq[9] ; int jj = (j)-4 ; \
383 if( jj < 0 ) jj = 0; else if( jj+8 >= num ) jj = num-9; \
384 qqq[0] = vec[jj+0]; qqq[1] = vec[jj+1]; qqq[2] = vec[jj+2]; \
385 qqq[3] = vec[jj+3]; qqq[4] = vec[jj+4]; qqq[5] = vec[jj+5]; \
386 qqq[6] = vec[jj+6]; qqq[7] = vec[jj+7]; qqq[8] = vec[jj+8]; \
387 med = median9f(qqq); qqq[0] = fabsf(qqq[0]-med); \
388 qqq[1] = fabsf(qqq[1]-med); qqq[2] = fabsf(qqq[2]-med); \
389 qqq[3] = fabsf(qqq[3]-med); qqq[4] = fabsf(qqq[4]-med); \
390 qqq[5] = fabsf(qqq[5]-med); qqq[6] = fabsf(qqq[6]-med); \
391 qqq[7] = fabsf(qqq[7]-med); qqq[8] = fabsf(qqq[8]-med); \
392 mad = median9f(qqq); }
393
394 /*-------------------------------------------------------------------------*/
395 /*! Remove spikes from a time series, in a very simplistic way.
396 Return value is the number of spikes that were squashed [RWCox].
397 *//*-----------------------------------------------------------------------*/
398
THD_despike9(int num,float * vec)399 int THD_despike9( int num , float *vec )
400 {
401 int ii , nsp ; float *zma,*zme , med,mad,val ;
402
403 if( num < 9 || vec == NULL ) return 0 ;
404 zme = (float *)malloc(sizeof(float)*num) ;
405 zma = (float *)malloc(sizeof(float)*num) ;
406
407 for( ii=0 ; ii < num ; ii++ ){
408 mead9(ii) ; zme[ii] = med ; zma[ii] = mad ;
409 }
410 mad = qmed_float(num,zma) ; free(zma) ;
411 if( mad <= 0.0f ){ free(zme); return 0; } /* should not happen */
412 mad *= 6.789f ; /* threshold value */
413
414 for( nsp=ii=0 ; ii < num ; ii++ )
415 if( fabsf(vec[ii]-zme[ii]) > mad ){ vec[ii] = zme[ii]; nsp++; }
416
417 free(zme) ; return nsp ;
418 }
419 #undef mead9
420
421 /*-------------------------------------------------------------------------*/
422 /* Return value: .i component = number of vectors with spikes
423 .j component = total number of spikes squashed.
424 *//*-----------------------------------------------------------------------*/
425
THD_vectim_despike9(MRI_vectim * mrv)426 int_pair THD_vectim_despike9( MRI_vectim *mrv )
427 {
428 int_pair pout = {0,0} ; int nv,ns , ss , kk ;
429
430 ENTRY("THD_vectim_despike9") ;
431
432 if( mrv == NULL || mrv->nvals < 9 ) RETURN(pout) ;
433
434 for( nv=ns=kk=0 ; kk < mrv->nvec ; kk++ ){
435 ss = THD_despike9( mrv->nvals , VECTIM_PTR(mrv,kk) ) ;
436 if( ss > 0 ){ nv++ ; ns += ss ; }
437 }
438
439 pout.i = nv ; pout.j = ns ; RETURN(pout) ;
440 }
441
442 /*-------------------------------------------------------------------------*/
443
THD_despike9_dataset(THD_3dim_dataset * inset,byte * mask)444 THD_3dim_dataset * THD_despike9_dataset( THD_3dim_dataset *inset , byte *mask )
445 {
446 THD_3dim_dataset *outset ;
447 MRI_vectim *mrv ;
448 int ii ;
449
450 ENTRY("THD_despike9_dataset") ;
451
452 if( !ISVALID_DSET(inset) || DSET_NVALS(inset) < 9 ) RETURN(NULL) ;
453
454 mrv = THD_dset_to_vectim(inset,mask,0) ; DSET_unload(inset) ;
455 if( mrv == NULL ) RETURN(NULL) ;
456
457 (void)THD_vectim_despike9(mrv) ;
458
459 outset = EDIT_empty_copy(inset) ;
460 for( ii=0 ; ii < DSET_NVALS(outset) ; ii++ )
461 EDIT_substitute_brick(outset,ii,MRI_float,NULL) ;
462
463 THD_vectim_to_dset(mrv,outset) ; VECTIM_destroy(mrv) ;
464 RETURN(outset) ;
465 }
466