Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP

Comparing changes

Choose two branches to see what's changed or to start a new pull request. If you need to, you can also compare across forks.

Open a pull request

Create a new pull request by comparing changes across two branches. If you need to, you can also compare across forks.
...
Checking mergeability… Don't worry, you can still create the pull request.
  • 5 commits
  • 5 files changed
  • 0 commit comments
  • 1 contributor
View
12 src/backend_common.c
@@ -236,7 +236,7 @@ void print_spectra_info(struct spectra_info *s)
printf(" bytes per subint = %d\n", s->bytes_per_subint);
printf(" samples per subint = %d\n", s->samples_per_subint);
printf(" zero offset = %-17.15g\n", s->zero_offset);
- printf(" Invert the band? = %s\n", s->apply_flipband ? "True" : "False");
+ printf(" Invert the band? = %s\n", (s->apply_flipband>0) ? "True" : "False");
if (s->header_offset!=NULL)
printf(" bytes in file header = %d\n", s->header_offset[0]);
if (s->datatype==PSRFITS) {
@@ -262,7 +262,7 @@ void print_spectra_info(struct spectra_info *s)
void print_spectra_info_summary(struct spectra_info *s)
// Print the basic details of the files that are being processed
{
- int ii;
+ int ii, nn;
printf(" Number of files = %d\n", s->num_files);
if (s->num_polns>=2 && !s->summed_polns)
printf(" Num of polns = %d\n", s->num_polns);
@@ -279,8 +279,9 @@ void print_spectra_info_summary(struct spectra_info *s)
printf(" Clipping sigma = %.3f\n", s->clip_sigma);
if (s->zero_offset!=0.0)
printf(" zero offset = %-17.15g\n", s->zero_offset);
- printf(" Invert the band? = %s\n", s->apply_flipband ? "True" : "False");
+ printf(" Invert the band? = %s\n", (s->apply_flipband>0) ? "True" : "False");
printf(" Byteswap? = %s\n", s->flip_bytes ? "True" : "False");
+ printf(" Remove zeroDM? = %s\n", s->remove_zerodm ? "True" : "False");
if (s->datatype==PSRFITS) {
printf(" Apply scaling? = %s\n", s->apply_scale ? "True" : "False");
printf(" Apply offsets? = %s\n", s->apply_offset ? "True" : "False");
@@ -288,7 +289,10 @@ void print_spectra_info_summary(struct spectra_info *s)
}
printf("\nFile Samples Padding Start MJD\n");
printf("---- ---------- ---------- --------------------\n");
- for (ii = 0; ii < s->num_files; ii++)
+ if (s->datatype==SUBBAND || s->datatype==DAT ||
+ s->datatype==EVENTS || s->datatype==SDAT) nn = 1;
+ else nn = s->num_files;
+ for (ii = 0; ii < nn; ii++)
printf("%-4d %10lld %10lld %19.14Lf\n", ii + 1,
s->num_spec[ii], s->num_pad[ii], s->start_MJD[ii]);
printf("\n");
View
39 src/clipping.c
@@ -3,8 +3,8 @@
#define BLOCKSTOSKIP 100
/* NEW Clipping Routine (uses channel running averages) */
-int new_clip_times(float *rawdata, int ptsperblk, int numchan,
- float clip_sigma, float *good_chan_levels)
+int clip_times(float *rawdata, int ptsperblk, int numchan,
+ float clip_sigma, float *good_chan_levels)
// Perform time-domain clipping of rawdata. This is a 2D array with
// ptsperblk*numchan points, each of which is a float. The clipping
// is done at clip_sigma sigma above/below the running mean. The
@@ -15,7 +15,7 @@ int new_clip_times(float *rawdata, int ptsperblk, int numchan,
static float running_avg = 0.0, running_std = 0.0;
static int blocksread = 0, firsttime = 1;
static long long current_point = 0;
- float *zero_dm_block, *median_temp;
+ float *zero_dm_block, *ftmp;
double *chan_avg_temp;
float current_med, trigger;
double current_avg = 0.0, current_std = 0.0;
@@ -28,17 +28,22 @@ int new_clip_times(float *rawdata, int ptsperblk, int numchan,
}
chan_avg_temp = gen_dvect(numchan);
zero_dm_block = gen_fvect(ptsperblk);
- median_temp = gen_fvect(ptsperblk);
+ ftmp = gen_fvect(ptsperblk);
/* Calculate the zero DM time series */
for (ii = 0; ii < ptsperblk; ii++) {
zero_dm_block[ii] = 0.0;
powptr = rawdata + ii * numchan;
for (jj = 0; jj < numchan; jj++)
- zero_dm_block[ii] += *powptr++;
- median_temp[ii] = zero_dm_block[ii];
+ zero_dm_block[ii] += rawdata[ii * numchan + jj];
+ //zero_dm_block[ii] += *powptr++;
+ ftmp[ii] = zero_dm_block[ii];
+ if (ii < 10) printf("%d %f", ii, ftmp[ii]);
}
- current_med = median(median_temp, ptsperblk);
+ printf("\n");
+ avg_var(ftmp, ptsperblk, &current_avg, &current_std);
+ current_std = sqrt(current_std);
+ current_med = median(ftmp, ptsperblk);
/* Calculate the current standard deviation and mean */
/* but only for data points that are within a certain */
@@ -48,20 +53,22 @@ int new_clip_times(float *rawdata, int ptsperblk, int numchan,
float lo_cutoff, hi_cutoff;
int numgoodpts = 0;
- lo_cutoff = 0.7 * current_med;
- hi_cutoff = 1.3 * current_med;
+ lo_cutoff = current_med - 3.0 * current_std;
+ hi_cutoff = current_med + 3.0 * current_std;;
for (jj = 0; jj < numchan; jj++)
chan_avg_temp[jj] = 0.0;
/* Find the "good" points */
for (ii = 0; ii < ptsperblk; ii++) {
if (zero_dm_block[ii] > lo_cutoff && zero_dm_block[ii] < hi_cutoff) {
- median_temp[numgoodpts] = zero_dm_block[ii];
+ ftmp[numgoodpts] = zero_dm_block[ii];
powptr = rawdata + ii * numchan;
for (jj = 0; jj < numchan; jj++)
chan_avg_temp[jj] += *powptr++;
numgoodpts++;
}
}
+//printf("med = %f std = %f numgoodpts = %d\n", current_med, current_std, numgoodpts);
+
/* Calculate the current average and stddev */
if (numgoodpts < 1) {
current_avg = running_avg;
@@ -69,7 +76,7 @@ int new_clip_times(float *rawdata, int ptsperblk, int numchan,
for (jj = 0; jj < numchan; jj++)
chan_avg_temp[jj] = chan_running_avg[jj];
} else {
- avg_var(median_temp, numgoodpts, &current_avg, &current_std);
+ avg_var(ftmp, numgoodpts, &current_avg, &current_std);
current_std = sqrt(current_std);
for (jj = 0; jj < numchan; jj++)
chan_avg_temp[jj] /= numgoodpts;
@@ -87,8 +94,8 @@ int new_clip_times(float *rawdata, int ptsperblk, int numchan,
running_std = current_std;
for (ii = 0; ii < numchan; ii++)
chan_running_avg[ii] = chan_avg_temp[ii];
- if (running_avg == 0.0 || current_avg == 0.0)
- printf("BAD RFI IN BLOCK#1!!!\n\n");
+ if (current_avg == 0.0)
+ printf("Warning: problem with clipping in first block!!!\n\n");
}
/* See if any points need clipping */
@@ -124,15 +131,15 @@ int new_clip_times(float *rawdata, int ptsperblk, int numchan,
vect_free(chan_avg_temp);
vect_free(zero_dm_block);
- vect_free(median_temp);
+ vect_free(ftmp);
return clipped;
}
/* OLD Clipping Routine (uses channel medians) */
-int clip_times(float *rawdata, int ptsperblk, int numchan,
- float clip_sigma, float *good_chan_levels)
+int old_clip_times(float *rawdata, int ptsperblk, int numchan,
+ float clip_sigma, float *good_chan_levels)
// Perform time-domain clipping of rawdata. This is a 2D array with
// ptsperblk*numchan points, each of which is a float. The clipping
// is done at clip_sigma sigma above/below the running mean. The
View
5 src/prepfold.c
@@ -194,7 +194,7 @@ int main(int argc, char *argv[])
printf("Reading %s data from 1 file:\n", description);
for (ii = 0; ii < s.num_files; ii++) {
printf(" '%s'\n", cmd->argv[ii]);
- if (insubs) s.files[ii] = chkfopen(s.filenames[0], "rb");
+ if (insubs) s.files[ii] = chkfopen(s.filenames[ii], "rb");
}
printf("\n");
if (RAWDATA) {
@@ -207,7 +207,7 @@ int main(int argc, char *argv[])
} else { // insubs
cmd->nsub = s.num_files;
s.N = chkfilelen(s.files[0], sizeof(short));
- ptsperrec = SUBSBLOCKLEN;
+ s.spectra_per_subint = ptsperrec = SUBSBLOCKLEN;
numrec = s.N / ptsperrec;
s.padvals = gen_fvect(s.num_files);
for (ii = 0 ; ii < s.num_files ; ii++)
@@ -252,6 +252,7 @@ int main(int argc, char *argv[])
s.hi_freq = s.lo_freq + (s.num_channels - 1.0) * s.df;
s.BW = s.num_channels * s.df;
s.fctr = s.lo_freq - 0.5 * s.df + 0.5 * s.BW;
+ print_spectra_info_summary(&s);
} else {
printf("\nThe input files (%s) must be subbands! (i.e. *.sub##)\n\n",
cmd->argv[0]);
View
9 src/prepsubband.c
@@ -143,7 +143,7 @@ int main(int argc, char *argv[])
printf("Reading %s data from 1 file:\n", description);
for (ii = 0; ii < s.num_files; ii++) {
printf(" '%s'\n", cmd->argv[ii]);
- if (insubs) s.files[ii] = chkfopen(s.filenames[0], "rb");
+ if (insubs) s.files[ii] = chkfopen(s.filenames[ii], "rb");
}
printf("\n");
if (RAWDATA) {
@@ -196,6 +196,7 @@ int main(int argc, char *argv[])
s.BW = s.num_channels * s.df;
s.fctr = s.lo_freq - 0.5 * s.df + 0.5 * s.BW;
s.spectra_per_subint = SUBSBLOCKLEN;
+ print_spectra_info_summary(&s);
} else {
printf("\nThe input files (%s) must be subbands! (i.e. *.sub##)\n\n",
cmd->argv[0]);
@@ -209,7 +210,7 @@ int main(int argc, char *argv[])
datafilenm = (char *) calloc(strlen(cmd->outfile) + 20, 1);
if (!cmd->subP) {
- printf("\nWriting output data to:\n");
+ printf("Writing output data to:\n");
outfiles = (FILE **) malloc(cmd->numdms * sizeof(FILE *));
dms = gen_dvect(cmd->numdms);
for (ii = 0; ii < cmd->numdms; ii++) {
@@ -230,7 +231,7 @@ int main(int argc, char *argv[])
" Setting the '-nobary' flag automatically.\n");
cmd->nobaryP = 1;
}
- printf("\nWriting subbands to:\n");
+ printf("Writing subbands to:\n");
cmd->numdms = 1;
dms = gen_dvect(cmd->numdms);
dms[0] = cmd->subdm;
@@ -317,7 +318,7 @@ int main(int argc, char *argv[])
/* Allocate our data array and start getting data */
- printf("De-dispersing using:\n");
+ printf("\nDe-dispersing using:\n");
printf(" Subbands = %d\n", cmd->nsub);
printf(" Average DM = %.7g\n", avgdm);
if (cmd->downsamp > 1) {
View
97 src/sigproc_fb.c
@@ -481,6 +481,8 @@ int get_filterbank_rawblock(float *fdata, struct spectra_info *s, int *padding)
{
int numread, numtopad = 0, numtoread;
float *fdataptr = fdata;
+ static float *bandpass = NULL;
+ static int firsttime = 1;
// If we only made a partial read last time, adjust our pointer
// offsetting into the output floating-point array
@@ -499,13 +501,11 @@ int get_filterbank_rawblock(float *fdata, struct spectra_info *s, int *padding)
}
convert_filterbank_block(fdataptr, cdatabuffer, numread, s);
- // Need to add -invert and zero-DMing
-
if (numread==numtoread) { // Got all we needed
*padding = 0;
numbuffered = 0;
currentblock++;
- return 1;
+ goto return_block;
} else { // Didn't get all the data we needed
numbuffered += numread;
if (feof(s->files[currentfile])) { // End of file?
@@ -531,7 +531,7 @@ int get_filterbank_rawblock(float *fdata, struct spectra_info *s, int *padding)
numpadded = 0;
currentfile++;
}
- return 1;
+ goto return_block;
} else { // Need < 1 block (or remaining block) of padding
int ii;
for (ii = 0; ii < numtopad; ii++)
@@ -551,6 +551,95 @@ int get_filterbank_rawblock(float *fdata, struct spectra_info *s, int *padding)
exit(1);
}
}
+
+return_block:
+ // Apply the corrections that need a full block
+
+ // Perform Zero-DMing if requested
+ if (s->remove_zerodm) {
+ int ii, jj, offset;
+ float zerodm, ftmp, padvalsum = 0.0, *chanwts;
+
+ // Make a static copy of the bandpass to use as DC offsets for each
+ // channel. This is necessary when masking as the newpadvals are
+ // used to mask. So our post-zero-DM-subtracted band had better look
+ // kinda like newpadvals. We want it static so that we don't get a
+ // different DC offset for each row of the FITS file (which will put
+ // a nice periodicity in the data).
+ if (firsttime) {
+ float bpsum = 0.0;
+ bandpass = gen_fvect(s->num_channels);
+ for (ii = 0; ii < s->num_channels; ii++) {
+ bandpass[ii] = s->padvals[ii];
+ bpsum += bandpass[ii];
+ }
+ if (bpsum==0.0) { // i.e. no padding is set
+ float favg;
+ printf("\nUsing first block channel averages for zeroDM bandpass.\n");
+ printf("Would be better to use statistics from an rfifind mask...\n\n");
+ for (ii = 0; ii < s->num_channels; ii++) {
+ favg = 0.0;
+ for (jj = 0; jj < s->spectra_per_subint; jj++)
+ favg += fdata[ii + jj * s->num_channels];
+ bandpass[ii] = favg / s->spectra_per_subint;
+ }
+ }
+ }
+
+ // Determine what weights to use for each channel for
+ // the zero-DM subtraction
+ chanwts = gen_fvect(s->num_channels);
+ if (s->clip_sigma > 0.0) {
+ // Determine the sum of newpadvals for weighting of the channels
+ for (ii = 0; ii < s->num_channels; ii++)
+ padvalsum += s->padvals[ii];
+ // Set the channel weights for this block
+ if (padvalsum!=0.0) {
+ for (ii = 0; ii < s->num_channels; ii++)
+ chanwts[ii] = (float) s->padvals[ii] / padvalsum;
+ } else {
+ for (ii = 0; ii < s->num_channels; ii++)
+ chanwts[ii] = 1.0 / s->num_channels;
+ }
+ } else {
+ for (ii = 0; ii < s->num_channels; ii++)
+ chanwts[ii] = 1.0 / s->num_channels;
+ }
+
+ // Now loop over all the spectra...
+ for (ii = 0; ii < s->spectra_per_subint; ii++) {
+ // Determine the DM=0 total power for this spectra
+ zerodm = 0.0;
+ offset = ii * s->num_channels;
+ for (jj = 0; jj < s->num_channels; jj++)
+ zerodm += fdata[offset+jj];
+ // Subtract off the correct amount of power from each point
+ for (jj = 0; jj < s->num_channels; jj++) {
+ // Put the constant bandpass back in since we are subtracting
+ // comparable sized numbers and don't (usually) want power < 0
+ ftmp = fdata[offset+jj] - chanwts[jj] * zerodm; // + bandpass[jj];
+ fdata[offset+jj] = ftmp;
+ }
+ }
+ vect_free(chanwts);
+ firsttime = 0;
+ }
+
+ // Invert the band if explicitly requested
+ if (s->apply_flipband==1) {
+ float ftmp;
+ int ii, jj, offset;
+ for (jj = 0 ; jj < s->spectra_per_subint ; jj++) {
+ offset = jj * s->num_channels;
+ for (ii = 0 ; ii < s->num_channels/2 ; ii++) {
+ ftmp = fdata[offset+ii];
+ fdata[offset+ii] = fdata[offset+s->num_channels-ii-1];
+ fdata[offset+s->num_channels-ii-1] = ftmp;
+ }
+ }
+ }
+
+ return 1;
}

No commit comments for this range

Something went wrong with that request. Please try again.