00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006    
00007 /* This file is #included by the filtering program 3dFourier and the plugin plug_fourier */
00008 /* By T. Ross and K. Heimerl 8-99 */
00009 
00010 
00011 static void *My_Malloc( size_t size) {
00012         void *ptr=NULL;
00013         
00014         ptr = (void *)malloc(size);
00015         if (ptr == NULL) {
00016 #ifdef IN_FOURIER_PLUGIN
00017                 fprintf(stderr, "Fatal error in Fourier Filter Driver, malloc returned NULL");
00018                 exit(1);
00019 #else
00020                 Error_Exit("Fatal: malloc returned NULL"); 
00021 #endif
00022         }
00023         return ptr;
00024 }
00025 
00026 
00027 
00028 
00029 /**************************************************************************************
00030 This function zero pads a signal to a power of 2 * factor of 3 or 5 then FFT's the data. 
00031 The zero-padded data are FFT'd, and the resultant spectrum is multiplied by an ideal
00032 low or high pass window.  The signal is then inverse FFt'd and the zero padded portion
00033 of the signal is eliminated.
00034 **************************************************************************************/
00035 static char *filter(float *ORIG_SIG, float low_fc, float high_fc, int N, float period, int ignore, int retrend, int transform)
00036 {
00037 
00038         /**************************************************************************************
00039                         DECLARE and INITIALIZE VARIABLES
00040         **************************************************************************************/
00041         float sum, meanorig, mean_new, Fs, slope, inter;
00042         int  dummy, i,j,power,padded_N, M;
00043         int Ncutlo1,Ncutlo2,Ncuthi1,Ncuthi2;
00044         FILE *outfile;
00045         complex *NEW_SIG, *newfft, *win_fft, *win_low, *win_high;  /*typedef'd by B Cox */
00046         complex *restored_sig, *final_sig, sumcx, meanfinal;
00047 
00048 
00049 
00050 
00051         /*************************************************************************************
00052                 DETREND ORIGINAL DATA 
00053         **************************************************************************************/
00054 
00055         if (ignore>=N) {
00056                 if (!transform)
00057                         return "You cannot ignore all of the data";
00058                 else
00059                         high_fc = 1000;  /* fine, make it an all-stop */
00060         } else {
00061                 for (i=0; i<(N-ignore); i++)
00062                         ORIG_SIG[i] = ORIG_SIG[i+ignore];
00063                 N -= ignore;
00064         }
00065 
00066         /* Detrend the data, setting the end points to zero */
00067         slope = ((float)ORIG_SIG[N-1] - (float)ORIG_SIG[0] ) / ((float) N - 1.0);
00068         inter = (float)ORIG_SIG[0];
00069         for (i=0; i<N ;++i)
00070                 ORIG_SIG[i] -= (inter + slope*(float)i);
00071 
00072         /*************************************************************************************
00073                 CALCULATE LENGTH TO WHICH DATA MUST BE ZERO PADDED for csfft_cox()
00074         **************************************************************************************/
00075 
00076 
00077         padded_N=csfft_nextup_one35(N);
00078 
00079 
00080         NEW_SIG=(complex *)My_Malloc(padded_N*sizeof(complex));         /*Tom Ross's malloc*/
00081         newfft=(complex *)My_Malloc(padded_N*sizeof(complex));          /*e.g. newfft[7].r or newfft[7].i*/
00082         win_fft=(complex *)My_Malloc(padded_N*sizeof(complex));
00083         win_low=(complex *)My_Malloc(padded_N*sizeof(complex));
00084         win_high=(complex *)My_Malloc(padded_N*sizeof(complex));
00085         restored_sig=(complex *)My_Malloc(padded_N*sizeof(complex));
00086         final_sig=(complex *)My_Malloc(N*sizeof(complex));              
00087 
00088 
00089         Fs= 1.0/period;                                                 /* sampling frequency*/
00090 
00091         M=padded_N;
00092         
00093 
00094         /*************************************************************************************
00095                         PERFORM FFT ON NEW DATA 
00096         **************************************************************************************/
00097 
00098         for(i=0;i<N;++i)
00099         {
00100                 newfft[i].r=ORIG_SIG[i];
00101                 newfft[i].i=0;
00102         }
00103         
00104         for(i=N;i<M;++i)                        /*Zero pad to length padded_N*/
00105         {
00106                 newfft[i].r=0;
00107                 newfft[i].i=0;
00108         }
00109         /* M must be pwr of 2 OR pwr of 2 times a factor of 3 or 5*/
00110 
00111         csfft_cox(-1, M, newfft);
00112 
00113 
00114         /*************************************************************************************
00115                         GENERATE LOW PASS WINDOW
00116         **************************************************************************************/
00117         if((low_fc >= Fs) || (low_fc == 0.0)) {
00118                 if ((low_fc >= Fs) && (!transform)) {
00119                         return "Lowpass filter is an all-pass filter since Fc > sampling frequency\n";
00120                 }
00121                 for(i=0;i<M;++i)
00122                 {
00123                         win_low[i].r= 1;
00124                         win_low[i].i= 0;
00125                 }
00126         } 
00127         else if ( (low_fc > 0) && (low_fc < (Fs/(float)M) ) )
00128         {
00129                 if (!transform) {
00130                         return "Lowpass filter is a no-pass filter since Fc is too small.\n";
00131                 }
00132 
00133                 for(i=0;i<M;++i)
00134                 {
00135                         win_low[i].r=0;  
00136                         win_low[i].i=0;
00137                 }
00138         } 
00139         else
00140         { 
00141 
00142         
00143                 Ncutlo1=(int)((low_fc*M)/Fs);   /*index number of lower cutoff frequency*/
00144                 Ncutlo2=(M-Ncutlo1);            /*index number of upper cutoff frequency*/
00145 
00146                 win_low[0].r=1;                 /*keep DC offset but no phase change*/
00147                 win_low[0].i=0;
00148 
00149                 for(i=1;i<=Ncutlo1;++i)
00150                 {
00151                         win_low[i].r=1;
00152                         win_low[i].i=0;
00153                 }
00154                 
00155                 for(i=Ncutlo1+1;i<Ncutlo2;++i)
00156                 {
00157                         win_low[i].r=0;  
00158                         win_low[i].i=0;
00159                 }
00160 
00161                 for(i=Ncutlo2;i<M;++i)
00162                 {
00163                         win_low[i].r=1;
00164                         win_low[i].i=0;
00165                 }
00166         }
00167 
00168 
00169         /*************************************************************************************
00170                         GENERATE HIGH PASS WINDOW
00171         **************************************************************************************/
00172         
00173         if (high_fc < Fs/(float)M) 
00174         {
00175                 if ((!transform) && (high_fc!=0)){
00176                         return "Highpass filter is an all-pass filter since Fc is too small.\n";
00177                 }
00178                 for(i=0;i<M;++i)
00179                 {
00180                         win_high[i].r= 1;
00181                         win_high[i].i= 0;
00182                 }       
00183         }
00184         else if (high_fc >= Fs) 
00185         {
00186                 if (!transform) {
00187                         return "Highpass filter is an all-stop filter since Fc > sampling frequency.\n";
00188                 }
00189 
00190                 for(i=0;i<M;++i)
00191                 {
00192                         win_low[i].r=0;  
00193                         win_low[i].i=0;
00194                 }
00195         } 
00196         else
00197         { 
00198                  
00199                 Ncuthi1=(int)((high_fc*M)/Fs);          /*index number of lower cutoff frequency*/
00200                 Ncuthi2=(M-Ncuthi1);                    /*index number of upper cutoff frequency*/
00201 
00202                 
00203                 for(i=0;i<Ncuthi1;++i)
00204                 {
00205                         win_high[i].r=0;
00206                         win_high[i].i=0;
00207                 }
00208                 
00209                 for(i=Ncuthi1;i<(int)(M/2);++i)
00210                 {
00211                         win_high[i].r=1;
00212                         win_high[i].i=0;
00213                 }
00214 
00215 
00216                         win_high[M/2].r=1;              
00217                         win_high[M/2].i=0;
00218 
00219 
00220 
00221                 for(i=(int)(M/2)+1;i<=Ncuthi2;++i)
00222                 {
00223                         win_high[i].r=1;
00224                         win_high[i].i=0;
00225                 }
00226 
00227 
00228                 for(i=Ncuthi2+1;i<M;++i)
00229                 {
00230                         win_high[i].r=0;
00231                         win_high[i].i=0;
00232                 }
00233         }
00234 
00235 
00236         /**************************************************************************************
00237                         MULTIPLY FFT'D NEWSIGNAL BY WINDOWS
00238         **************************************************************************************/
00239 
00240         /* first multiply through by low pass window  */
00241 
00242         for(i=0;i<M;++i)
00243         {
00244                 win_fft[i].r=win_low[i].r*newfft[i].r-win_low[i].i*newfft[i].i;
00245                 win_fft[i].i=win_low[i].r*newfft[i].i+win_low[i].i*newfft[i].r;
00246         }
00247 
00248 
00249 
00250         /* then multiply result by high pass window  */
00251 
00252         for(i=0;i<M;++i)
00253         {
00254                 restored_sig[i].r=win_high[i].r*win_fft[i].r-win_high[i].i*win_fft[i].i;
00255                 restored_sig[i].i=win_high[i].r*win_fft[i].i+win_high[i].i*win_fft[i].r;
00256         }
00257 
00258 
00259         /**************************************************************************************
00260                         INVERSE FFT WINDOWED NEW SIGNAL 
00261         **************************************************************************************/
00262 
00263         csfft_cox(1, M, restored_sig);          /* Doesn't scale for inverse FFT  */
00264 
00265         
00266         for(i=0;i<M;++i)
00267         {
00268                 restored_sig[i].r=restored_sig[i].r/(float)M;
00269                 restored_sig[i].i=restored_sig[i].i/(float)M;
00270         }
00271 
00272                 
00273 
00274         /**************************************************************************************
00275                         ORIG_SIGNAL IS RETURNED TO AFNI, FREE ALLOCATED MEMORY
00276         **************************************************************************************/
00277 
00278         /*Note that data points N-M are the zero padded portion - we don't keep those data*/
00279 
00280         for (i=0;i<N; ++i)
00281         {
00282                 ORIG_SIG[i]=restored_sig[i].r;
00283 
00284         }
00285 
00286 
00287         /* Retrend the data */
00288         if (retrend) 
00289                 for (i=0; i<N ;++i)
00290                         ORIG_SIG[i] += (inter + slope*(float)i);
00291         
00292         if (ignore) {
00293                 for (i=(N+ignore-1); i>=ignore; i--)
00294                         ORIG_SIG[i]=ORIG_SIG[i-ignore];
00295                 for (i=0; i<ignore; i++)
00296                         ORIG_SIG[i]=ORIG_SIG[ignore];
00297         }
00298 
00299 
00300         free(NEW_SIG);          
00301         free(newfft);   
00302         free(win_fft);
00303         free(win_low);
00304         free(win_high);
00305         free(restored_sig);
00306         free(final_sig);                
00307 
00308 
00309 
00310 
00311 
00312         /*************************************************************************************
00313                                         the end of the filter
00314         **************************************************************************************/
00315         return NULL;  /* No errors */
00316 }
00317 
00318 
00319 
00320 #ifdef IN_FOURIER_PLUGIN
00321 static char *Fourier_Filter_Driver(PLUGIN_interface *plint, THD_3dim_dataset *input, float low_fc, float high_fc, int ignore, int autocorr, int retrend, char *output_prefix)  {
00322 #else
00323 static char *Fourier_Filter_Driver(THD_3dim_dataset *input, float low_fc, float high_fc, int ignore, int autocorr, int retrend, char *output_prefix)  {
00324 #endif
00325         int i, j;
00326         int ntimes, nvox;
00327         float **out_data, *out_temp, *scale, *input_data_f, period;
00328         byte *input_data_b;
00329         short *input_data_s;
00330         THD_3dim_dataset *output=NULL;
00331         char *err;
00332 
00333         /* should be a valid 3d+time input */   
00334         DSET_load(input);
00335         
00336         ntimes = DSET_NUM_TIMES(input);
00337         nvox = DSET_NVOX(input);
00338         
00339         /* Create a float array for the output */
00340         out_data = (float **)My_Malloc(ntimes * sizeof(float *));
00341         for (i=0; i<ntimes; out_data[i++] = (float *)My_Malloc(nvox * sizeof(float)));
00342         
00343         /* Create the tempory float array */
00344         out_temp = (float *)My_Malloc(ntimes*sizeof(float));
00345         
00346         /* Get the scale factors for later */
00347         scale = (float *)My_Malloc(ntimes*sizeof(float));
00348         for (i=0; i<ntimes; i++) {
00349                 if (!DSET_BRICK_FACTOR(input,i))
00350                         scale[i] = 1.0;
00351                 else
00352                         scale[i] = DSET_BRICK_FACTOR(input,i);
00353         }
00354         
00355         /* get the sampling period */
00356         period = DSET_TIMESTEP(input);
00357         switch (DSET_TIMEUNITS(input)) {
00358                 case UNITS_MSEC_TYPE: {period/=1000; break; }
00359                 case UNITS_SEC_TYPE:   break;
00360                 case UNITS_HZ_TYPE:   return ("FIlter_Driver Error:\nHmm, you would think I could handle a 3d+time that is already\nfrequency, but Im lame");
00361                 default: return("Filter_Driver Error:\nI dont understand the units of the 3d+time");
00362         }
00363 
00364 #ifdef IN_FOURIER_PLUGIN
00365         PLUTO_popup_meter(plint);               
00366 #endif                  
00367                         
00368         /* Loop over voxels, pull out the time series and filter */
00369         for (i=0; i< nvox; i++) {
00370                 for (j=0; j<ntimes; j++) {
00371                         switch(DSET_BRICK_TYPE(input,j)) {
00372                                 case MRI_byte: {
00373                                         input_data_b = (byte *)DSET_ARRAY(input,j);
00374                                         out_temp[j] = scale[j] * (float)input_data_b[i];
00375                                         break;
00376                                 }
00377 
00378                                 case MRI_short: {
00379                                         input_data_s = (short *)DSET_ARRAY(input,j);
00380                                         out_temp[j] = scale[j] * (float)input_data_s[i];
00381                                         break;
00382                                 }
00383 
00384                                 case MRI_float: {
00385                                         input_data_f = (float *)DSET_ARRAY(input,j);
00386                                         out_temp[j] = scale[j] * input_data_f[i];
00387                                         break;
00388                                 }
00389                                 
00390                                 default : {
00391                                         return("FIlter_Driver Error:\nInvalid data type for one of the sub-bricks");
00392                                 }
00393                         }
00394                 }
00395                                         
00396                 err = filter(out_temp, low_fc, high_fc, ntimes, period, ignore, retrend, FALSE);
00397                 if (err != NULL)
00398                         return err;
00399                 
00400                 for(j=0; j<ntimes; j++)
00401                         out_data[j][i] = out_temp[j];
00402 #ifdef IN_FOURIER_PLUGIN                
00403                 PLUTO_set_meter(plint, (int)(100.0*((float)i/(float)nvox)));
00404 #endif
00405                         
00406         }
00407         
00408         /* create new dataset and convert, etc. */
00409         output = EDIT_empty_copy(input);
00410         
00411         j=EDIT_dset_items(output,
00412                 ADN_prefix, output_prefix,
00413                 ADN_none);
00414         
00415         for (j=0; j<ntimes; j++) {
00416                 switch(DSET_BRICK_TYPE(input,j)) {
00417                         case MRI_byte: {
00418                                 input_data_b = (byte *)My_Malloc(nvox*sizeof(byte));
00419                                 for (i=0; i<nvox; i++) 
00420                                         input_data_b[i] = (byte)(out_data[j][i] / scale[j]);
00421                                 EDIT_substitute_brick(output, j, MRI_byte, (byte *)input_data_b); 
00422                                 break;
00423                         }       
00424                         case MRI_short: {
00425                                 input_data_s = (short *)My_Malloc(nvox*sizeof(short));
00426                                 for (i=0; i<nvox; i++) 
00427                                         input_data_s[i] = (short)(out_data[j][i] / scale[j]);
00428                                 EDIT_substitute_brick(output, j, MRI_short, (short *)input_data_s); 
00429                                 break;
00430                         }       
00431                         case MRI_float: {
00432                                 input_data_f = (float *)My_Malloc(nvox*sizeof(float));
00433                                 for (i=0; i<nvox; i++) 
00434                                         input_data_f[i] = (float)(out_data[j][i] / scale[j]);
00435                                 EDIT_substitute_brick(output, j, MRI_float, (float *)input_data_f); 
00436                                 break;
00437                         }       
00438                 }
00439 #ifdef IN_FOURIER_PLUGIN                
00440                 PLUTO_set_meter(plint, (int)(100.0*((float)j/(float)ntimes)));
00441 #endif
00442 }
00443         
00444         /* Write out the new brick at let the memory be free */
00445 #ifdef IN_FOURIER_PLUGIN
00446         PLUTO_add_dset(plint, output, DSET_ACTION_MAKE_CURRENT);
00447 #else
00448         DSET_write(output);
00449 #endif
00450         DSET_unload(input);
00451         DSET_unload(output);
00452         for (i=0; i<ntimes; free(out_data[i++]));
00453         free (out_data);
00454         free (scale);
00455         free (out_temp);
00456 #ifdef IN_FOURIER_PLUGIN
00457         PLUTO_popdown_meter(plint);             
00458 #endif                  
00459 
00460         return NULL;
00461 }
00462 
00463 
00464 
00465 
00466