Skip to content

Commit

Permalink
Add opts for better control over RFI excision
Browse files Browse the repository at this point in the history
  • Loading branch information
@ymaan4
ymaan4 committed Feb 9, 2021
1 parent a147da5 commit 381422b
Show file tree
Hide file tree
Showing 5 changed files with 148 additions and 74 deletions.
14 changes: 7 additions & 7 deletions bin/crp_rficlean_fil.sh
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
#!/bin/bash
#=======================================================================
# Crude parallelization of rfiClean for processing filterbank files.
# Crude parallelization of RFIClean for processing filterbank files.
#
# Yogesh Maan, Dec. 2019
#=======================================================================
Expand Down Expand Up @@ -39,13 +39,13 @@ Nb=`echo "scale=0; $nsamples/($Np*$block) + 1" | bc`



echo "Starting parallel-rfiClean processing at: "
echo "Starting parallel-RFIClean processing at: "
date
echo ""

##-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
# prepare and execute the commands for all parts
mkdir -p rfiClean_ps
mkdir -p RFIClean_ps
combine_cmd="cat "
remove_cmd="rm "
count=0
Expand All @@ -57,15 +57,15 @@ for (( c=1; c<=$Np; c++ )); do
psfile=${rtag}_part${count}.ps
#-------------------------------------
if [ $count -eq 1 ]; then
cmd="rficlean -t $block $flag -o $outfile -ps rfiClean_ps/$psfile $infile -bst $bst -nbl $Nb &"
cmd="rficlean -t $block $flag -o $outfile -ps RFIClean_ps/$psfile $infile -bst $bst -nbl $Nb &"
elif [ $count -eq $Np ]; then
combine_cmd="${combine_cmd} ${tempout}"
remove_cmd="${remove_cmd} ${tempout}"
cmd="rficlean -t $block $flag -o $tempout -ps rfiClean_ps/$psfile $infile -bst $bst -headerless &"
cmd="rficlean -t $block $flag -o $tempout -ps RFIClean_ps/$psfile $infile -bst $bst -headerless &"
else
combine_cmd="${combine_cmd} ${tempout}"
remove_cmd="${remove_cmd} ${tempout}"
cmd="rficlean -t $block $flag -o $tempout -ps rfiClean_ps/$psfile $infile -bst $bst -nbl $Nb -headerless &"
cmd="rficlean -t $block $flag -o $tempout -ps RFIClean_ps/$psfile $infile -bst $bst -nbl $Nb -headerless &"
fi
echo $cmd
eval $cmd
Expand All @@ -83,7 +83,7 @@ echo $cmd
eval $cmd


echo "Finished parallel-rfiClean processing at: "
echo "Finished parallel-RFIClean processing at: "
date
echo ""

Expand Down
2 changes: 2 additions & 0 deletions include/rficlean.h
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
#include <stdlib.h>
#include <stdio.h>
#include <stdbool.h>
#include <math.h>
#include "header.h"
#include <fftw3.h>
Expand All @@ -14,6 +15,7 @@ float *fftstat, *chanstat, *predist, *xpredist, *postdist, *xpostdist, *finaldis
double *tfvar, *tfmean,meanvar,rmsvar ;
double *chandata, *mspec, *rspec, *vspec, *wspec, *wt;
long int *coff;
bool RFIx,rfiFDx,rfiTx,rfiSx,rfiMSx,rfiVSx,rfiSclip;

fftw_complex *in;
fftw_complex *out;
Expand Down
89 changes: 53 additions & 36 deletions src/cleanit.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -371,11 +371,19 @@ void cleanit(float *data, int nchans, long int nadd)
isame = all_samef(chandata,nadd);
}
if( isame == 0){
fftclean(in,out,nadd,inc); // clean some periodic RFIs
if(rfiFDx){
fftclean(in,out,nadd,inc); // clean some periodic RFIs
}
for (ii = 0; ii<nadd; ii++) chandata[ii] = in[ii][0];
tsclip(chandata,nadd,sthresh); // clean some spiky RFIs
mspec[channum]=chandata[nadd];
rspec[channum]=chandata[nadd+1];
if(rfiTx){
tsclip(chandata,nadd,sthresh); // clean some spiky RFIs
mspec[channum]=chandata[nadd];
rspec[channum]=chandata[nadd+1];
} else {
robust_meanrms(chandata,nadd);
mspec[channum]=chandata[nadd];
rspec[channum]=chandata[nadd+1];
}
}
else {
mspec[channum]=chandata[0];
Expand All @@ -388,10 +396,14 @@ void cleanit(float *data, int nchans, long int nadd)


for (i=0; i<nchans; i++) wspec[i]=+1.0;
spfind(mspec,nchans,rthresh,wspec);
spfind(vspec,nchans,rthresh,wspec);
tsfind(mspec,nchans,rthresh,wspec);
tsfind(vspec,nchans,rthresh,wspec);
if(rfiMSx){
spfind(mspec,nchans,rthresh,wspec);
tsfind(mspec,nchans,rthresh,wspec);
}
if(rfiVSx){
spfind(vspec,nchans,rthresh,wspec);
tsfind(vspec,nchans,rthresh,wspec);
}
if(pcl>0){
kk = (int)(0.9*nchans);
for (i=kk;i<nchans;i++) wspec[i]=0.0;
Expand Down Expand Up @@ -437,28 +449,32 @@ void cleanit(float *data, int nchans, long int nadd)


// Try clipping some channels in individual samples
for (t=0; t<nadd; t++){
nxc = t*nchans;
for (c=0; c<nchans; c++) chandata[c] = data[nxc+c];
spclip(chandata,nchans,clipthresh);
for (c=0; c<nchans; c++) data[nxc+c] = chandata[c];
if(rfiSclip){
for (t=0; t<nadd; t++){
nxc = t*nchans;
for (c=0; c<nchans; c++) chandata[c] = data[nxc+c];
spclip(chandata,nchans,clipthresh);
for (c=0; c<nchans; c++) data[nxc+c] = chandata[c];
}
}
// Now some timeseries cleaning
an = (double)nadd;
for (t=0; t<nadd; t++){
nxc = t*nchans;
chandata[t]=0.0;
for (c=0; c<nchans; c++) chandata[t] = chandata[t]+data[nxc+c];
chandata[t]=chandata[t]/an;
}
for (i=0; i<nadd; i++) wt[i]=+1.0;
tsfind(chandata,nadd,sthresh,wt);
spfind(chandata,nadd,sthresh,wt);
for (t=0; t<nadd; t++){
nxc = t*nchans;
if(wt[t] < 0.0){
for (c=0; c<nchans; c++) data[nxc+c] = mspec[c];
}
if(rfiTx){
an = (double)nadd;
for (t=0; t<nadd; t++){
nxc = t*nchans;
chandata[t]=0.0;
for (c=0; c<nchans; c++) chandata[t] = chandata[t]+data[nxc+c];
chandata[t]=chandata[t]/an;
}
for (i=0; i<nadd; i++) wt[i]=+1.0;
tsfind(chandata,nadd,sthresh,wt);
spfind(chandata,nadd,sthresh,wt);
for (t=0; t<nadd; t++){
nxc = t*nchans;
if(wt[t] < 0.0){
for (c=0; c<nchans; c++) data[nxc+c] = mspec[c];
}
}
}

// sometimes gpt results might need additional checks
Expand Down Expand Up @@ -536,14 +552,15 @@ void cleanit(float *data, int nchans, long int nadd)
}



if(iflip==1){ // flip the band
for (t=0; t<nadd; t++){
nxc = t*nchans;
for (c=0; c<nchans; c++) mspec[c] = data[nxc+c];
for (c=0; c<nchans; c++) data[nxc+c] = mspec[nchans-c-1];
}
}
/* ***Now the following part is in rficlean_data.c ***
// if(iflip==1){ // flip the band
// for (t=0; t<nadd; t++){
// nxc = t*nchans;
// for (c=0; c<nchans; c++) mspec[c] = data[nxc+c];
// for (c=0; c<nchans; c++) data[nxc+c] = mspec[nchans-c-1];
// }
// }
***/

}
//============================================================================
89 changes: 64 additions & 25 deletions src/rficlean.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,35 +42,53 @@ int help_required(char *string)
void rficlean_help()
{
puts("");
puts("rficlean - clean some periodic and/or narrow-band RFIs in filterbank data\n");
puts("rficlean - clean some periodic and/or narrow-band and/or spiky/bursty RFIs in filterbank data\n");
puts("usage: rficlean -{options} {input-filename} \n");
puts("options:\n");
puts("");
puts("General options:\n");
puts(" <filename> - filterbank data file (def=stdin)");
puts("-t <numsamps> - number of time samples in a block (for periodicity search; def=4096)");
puts("-t <numsamps> - number of time samples in a block (def=4096)");
puts("-ft <thres.> - specify threshold for FT-cleaning (def=6.0)");
puts("-st <thres.> - specify threshold for timeseries-cleaning (def=10.0)");
puts("-rt <thres.> - specify threshold for chan-diff cleaning (def=4.0)");
puts("-n <numbits> - specify output number of bits (def=input)");
puts("-white - whitten each spectrum before RFI excision");
puts("-zeordm - remove 0-dm powers before writing out data");
puts("-pcl - assume that input data may have some parts replaced");
puts(" by 0s or mean/median by some other RFI excision program");
puts("-nharm <nh> - Number of fundamental+harmonics to be removed");
puts("-nharm <nh> - Number of fundamental+harmonics to be excised");
puts(" (def=1, i.e., only fundamental)");
puts("-T <nsamp> - number of time samples to avg before output (def=0)");
puts("-gm <filename> - For gmrt raw-data, name of the file that contains header information");
puts("-gmtstamp <filename> - For gmrt data, name of the file that contains time-stamp information");
puts("-st <thres.> - specify threshold for timeseries cleaning (def=10.0)");
puts("-rt <thres.> - specify threshold for channel cleaning (def=4.0)");
puts("-clt <thres.> - specify threshold for channel clipping (def=5.0)");
puts("-ps <filename> - File-name of the output diagnostic plot (def=rficlean_output.ps)");
puts("-o <filename> - specify output filename (def=stdout)");
puts("");
puts("To safeguard a known periodic signal:");
puts("-psrf <F0> - fundamental freq. (Hz) of pulsar to be protected (def=none)");
puts("-psrfdf <dF> - Delta-frequency (on either side) of the fundamental & harmonic");
puts("-psrfbins <Nb> - *Nbins (on either side) of the fundamental & harmonic");
puts(" (*specify either psrfbins or psrfdf!)");
puts(" frequencies to be protected (def=2)");
puts("-psrfdf* <dF> - Delta-frequency (on either side) of the fundamental & harmonic");
puts(" frequencies to be protected (def=default psrfbins below)");
puts("-psrfbins* <Nb> - *Nbins (on either side) of the fundamental & harmonic");
puts(" frequencies to be protected (def=8)");
puts(" (*specify either psrfbins or psrfdf!)");
puts("");
puts("Add-on options:");
puts("-n <numbits> - specify output number of bits (def=input)");
puts("-zeordm - remove 0-dm powers before writing out data");
puts("-pcl - assume that input data may have some parts replaced");
puts(" by 0s or mean/median by some other RFI excision program");
puts("-T <nsamp> - number of time samples to avg before output (def=0)");
puts("-headerless - do not broadcast resulting header (def=broadcast)");
puts("-bst <bstart> - Starting block number to be processed (def=1, i.e. start of file)");
puts("-nbl <nblocks> - No. of blocks to be processed (def=till EoF)");
puts("");
puts("For GMRT data in native format:");
puts("-gm <filename> - For gmrt raw-data, name of the file that contains header information");
puts("-gmtstamp <filename> - For gmrt data, name of the file that contains time-stamp information");
puts("");
puts("Options to control which RFI excision methods to skip (if at all):");
puts("-noRFIx - Do not excise any RFI ");
puts("-noFDx - Do not excise periodic RFI in the Fourier domain");
puts("-noTx - Do not excise spiky RFI from channel/band-averaged time-series");
puts("-noSx - Do not excise narrow-band RFI from spectra");
puts("-noMSx - Do not excise narrow-band RFI using mean spectra");
puts("-noVSx - Do not excise narrow-band RFI using variance spectra");
puts("-noSclip - Do not clip above clip-threshold in individual spectra");
puts("");
}

int file_exists(char *filename)
Expand Down Expand Up @@ -99,7 +117,7 @@ void print_version(char *program, char *argument)
{
if ( (strings_equal(argument,"-v")) ||
(strings_equal(argument,"--version"))) {
printf("PROGRAM: %s is part of rfiClean version: %.1f\n",program,RFICLEAN_VERSION);
printf("PROGRAM: %s is part of RFIClean version: %.1f\n",program,RFICLEAN_VERSION);
exit(0);
}
}
Expand All @@ -112,6 +130,7 @@ void main (int argc, char *argv[])

/* set up default global variables */
obits=headerless=naddt=nsamp=0;
RFIx=rfiFDx=rfiTx=rfiSx=rfiMSx=rfiVSx=rfiSclip=true;
ibits=0;
fthresh = 6.0;
forcefthresh = 1000.0;
Expand Down Expand Up @@ -146,6 +165,22 @@ void main (int argc, char *argv[])
if (strings_equal(argv[i],"-t")) {
i++;
naddt=atoi(argv[i]);
} else if (strings_equal(argv[i],"-noRFIx")) {
RFIx = false;
} else if (strings_equal(argv[i],"-noFDx")) {
rfiFDx = false;
} else if (strings_equal(argv[i],"-noTx")) {
rfiTx = false;
} else if (strings_equal(argv[i],"-noSx")) {
rfiSx = false;
rfiMSx = false;
rfiVSx = false;
} else if (strings_equal(argv[i],"-noMSx")) {
rfiMSx = false;
} else if (strings_equal(argv[i],"-noVSx")) {
rfiVSx = false;
} else if (strings_equal(argv[i],"-noSclip")) {
rfiSclip = false;
} else if (strings_equal(argv[i],"-normf")) {
fnorm=1;
} else if (strings_equal(argv[i],"-normt")) {
Expand Down Expand Up @@ -186,6 +221,9 @@ void main (int argc, char *argv[])
} else if (strings_equal(argv[i],"-st")) {
i++;
sthresh=atof(argv[i]);
} else if (strings_equal(argv[i],"-clt")) {
i++;
clipthresh=atof(argv[i]);
} else if (strings_equal(argv[i],"-n")) {
i++;
obits=atoi(argv[i]);
Expand Down Expand Up @@ -238,24 +276,25 @@ void main (int argc, char *argv[])
rficlean_help();
exit(0);
}
if(!rfiVSx && !rfiMSx) rfiSx = false;

// some sanity checks
if (psrf<1000000.0 && fthresh<4.0 && forcefthresh>0) fthresh=4.0;
if (psrfdf>0.0 && psrfbins > 0){
printf ("Both psrfdf and psrfbins are specified!!");
printf ("Using the delta-F information from psrfdf.");
printf ("\nBoth psrfdf and psrfbins are specified!!\n");
printf ("Using the delta-F information from psrfdf.\n");
psrfbins = -1;
}
if (psrfdf<0.0 && psrfbins<0){
printf ("Nether psrfdf nor psrfbins is specified!!");
printf ("Using psrfbins=8");
if (psrf<100000.0 && psrfdf<0.0 && psrfbins<0){
printf ("\nNether psrfdf nor psrfbins is specified!!\n");
printf ("Using psrfbins=8\n");
psrfbins = 8;
}

/* read in the header to establish what the input data are... */
if (nifs>1){
printf ("Cannot process data with IFs more than 1.");
printf ("Halting !!");
printf ("Cannot process data with IFs more than 1.\n");
printf ("Halting !!\n");
exit(0);
}
if (gm>0) {
Expand Down
28 changes: 22 additions & 6 deletions src/rficlean_data.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -152,7 +152,7 @@ void rficlean_data(FILE *input, FILE *output)
printf (" All buffers prepared! \n\n");

/* main loop */
printf (" Now rfiClean-ing (and 0-DM cleaning & downsampling, if asked for), and writing out the data...\n \n");
printf (" Now RFIClean-ing (and 0-DM cleaning & downsampling, if asked for), and writing out the data...\n \n");
istart = 0;
iblock = 0;
while ((ns=read_block(input,nbits,fblock,nsblk,byte_offset))>0 && iblock<nblocks) {
Expand All @@ -163,7 +163,18 @@ void rficlean_data(FILE *input, FILE *output)
jj = naddt*nchans;
for (j=ns;j<jj;j++) fblock[j] = 0.0;
}
cleanit(fblock,nchans,naddt);
if(RFIx){
cleanit(fblock,nchans,naddt);
}

// flip the band, if needed
if(iflip==1){
for (k=0; k<naddt; k++){
jj = k*nchans;
for (j=0; j<nchans; j++) mspec[j] = fblock[jj+j];
for (j=0; j<nchans; j++) fblock[jj+j] = mspec[nchans-j-1];
}
}
//-------------------------------------------------
// compute post-cleaning 0-DM tseries
for (j=0;j<n0;j++) {
Expand Down Expand Up @@ -246,11 +257,16 @@ void rficlean_data(FILE *input, FILE *output)
//sprintf(string,"time:%.1fs",realtime);
}

printf (" Data rfiClean-ed and written out! \n\n");
printf (" Now making diagnostic plots ... ");
if(RFIx){
printf (" Data RFIClean-ed and written out! \n\n");
printf (" Now making diagnostic plots ... ");

plot_data(plotdevice,wpout);
printf (" Done! \n ");
plot_data(plotdevice,wpout);
printf (" Done! \n ");
} else {
printf (" Data not RFIClean-ed,\n");
printf (" other desired operation, if any, performed and data written out! \n\n");
}

fclose(input);
fclose(output);
Expand Down

0 comments on commit 381422b

Please sign in to comment.