Doxygen Source Code Documentation
fourier_filter.c File Reference
Go to the source code of this file.
Functions | |
| void * | My_Malloc (size_t size) |
| char * | filter (float *ORIG_SIG, float low_fc, float high_fc, int N, float period, int ignore, int retrend, int transform) |
| char * | Fourier_Filter_Driver (THD_3dim_dataset *input, float low_fc, float high_fc, int ignore, int autocorr, int retrend, char *output_prefix) |
Function Documentation
|
||||||||||||||||||||||||||||||||||||
|
Definition at line 35 of file fourier_filter.c. References csfft_cox(), csfft_nextup_one35(), dummy, free, complex::i, i, My_Malloc(), and complex::r. Referenced by fourier_1D(), and Fourier_Filter_Driver().
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 }
|
|
||||||||||||||||||||||||||||||||
|
Definition at line 323 of file fourier_filter.c. References ADN_none, ADN_prefix, DSET_ARRAY, DSET_BRICK_FACTOR, DSET_BRICK_TYPE, DSET_load, DSET_NUM_TIMES, DSET_NVOX, DSET_TIMESTEP, DSET_TIMEUNITS, DSET_unload, DSET_write, EDIT_dset_items(), EDIT_empty_copy(), EDIT_substitute_brick(), filter(), free, i, My_Malloc(), PLUTO_add_dset(), PLUTO_popdown_meter(), PLUTO_popup_meter(), PLUTO_set_meter(), scale, UNITS_HZ_TYPE, UNITS_MSEC_TYPE, and UNITS_SEC_TYPE. Referenced by Fourier_Main(), and main().
00323 {
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 }
|
|
|
Definition at line 11 of file fourier_filter.c. References Error_Exit(), and malloc. Referenced by filter(), Fourier_Filter_Driver(), and main().
00011 {
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 }
|
