Last active
January 12, 2022 17:36
-
-
Save yongrenjie/79b7ba39fcfe548af624dcb50b69246e to your computer and use it in GitHub Desktop.
Substantially faster version of pshift AU programme for 2D data
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
/* min_pshift | |
* ---------- | |
* | |
* Minimal version of pshift. Only deals with 2D data. Lots of improvements to | |
* code, so runs much faster especially when the entire 2D dataset doesn't fit | |
* into memory. | |
* | |
* Jonathan Yong, University of Oxford | |
* 31 Oct 2020 */ | |
// Maximum size of TD that we can handle. | |
#define MAX_TD 131072 | |
// Macro to round doubles or floats to integer. | |
#define ROUND(x) ((x) < 0 ? (int) ((x) - 0.5) : (int) ((x) + 0.5)) | |
// New expno to store the processed data in. | |
int new_expno = expno + 1000; | |
GETCURDATA | |
/* The GRPDLY parameter tells us how much to circularly shift the FID by (i.e. | |
* take points from the start and move them to the end; and then apply a | |
* first-order phase correction). The formula is as such. The number of points | |
* to be cut is floor(GRPDLY), and the 1st order phase correction to be applied | |
* is (GRPDLY - floor(GRPDLY)) * 360. However, in practice, the phase correction | |
* is quite small and so can be left to apk to deal with. | |
* Here we cut a ton of code that was to do with not being able to read GRPDLY | |
* on older spectrometers (I don't know *how* old). It basically involves | |
* reading the parameters DECIM and DSPFVS, then performing a lookup in a list, | |
* which is provided at the bottom of the original pshift script. I don't know | |
* if it's worth (re-)implementing this. | |
*/ | |
double grpdly; | |
FETCHPARS("GRPDLY", &grpdly) | |
int grpdly_points = 2 * ((int) (grpdly + 1)); // *2 because real + imaginary. I don't know why +1. | |
/* Run 2D processing. Error out if not 2D. */ | |
/* TODO: does this work on new TopSpin 4 double-format data? */ | |
int pmod; | |
FETCHPAR("PARMODE", &pmod) | |
if (pmod != 1) STOPMSG("min_pshift: only works on 2D data") | |
// Read in key parameters | |
int td1, td2; | |
double sw1, sw2, sfo1, sfo2; | |
FETCHPAR1S("TD", &td1) | |
FETCHPAR("TD", &td2) | |
FETCHPAR1S("SW", &sw1) | |
FETCHPAR1S("SFO1", &sfo1) | |
FETCHPARS("SW", &sw2) | |
FETCHPARS("SFO1", &sfo2) // this is not actually SFO2 but rather F2 SFO1 | |
sw1 = sw1 * sfo1; // convert both SW's to Hz | |
sw2 = sw2 * sfo2; | |
// ser files are stored in a way such that each individual FID takes up 256n | |
// data points, even if TD is not a multiple of 256. So, we need to round TD2 | |
// up to the nearest multiple of 256. | |
td2 = (td2 + 255) / 256 * 256; | |
// Files to read from and write to. | |
char infile[PATH_MAX], outfile[PATH_MAX]; | |
strcpy(infile, ACQUPATH("ser")); | |
strcpy(outfile, NEWACQUPATH(new_expno, "fid")); | |
// Create the new dataset in the next expno | |
RSER(1, new_expno, 1) | |
// Get the appropriate number of drop points. | |
float cnst4; | |
FETCHPAR("CNST 4", &cnst4) | |
int drop_points = 2 * ROUND(cnst4); // *2 because complex points | |
// Get the number of complex points per chunk. | |
// Note that if sw2/sw1 is not an exact integer, it gets rounded: this means that in effect | |
int points_per_chunk = 2 * ROUND(sw2/sw1); | |
int ser_array[MAX_TD]; | |
FILE *fp_ser, *fp_fid; | |
if ((fp_ser = fopen(infile, "r")) == NULL) STOPMSG("min_pshift: failed to open ser file for reading") | |
if ((fp_fid = fopen(outfile, "w")) == NULL) STOPMSG("min_pshift: failed to open fid file for writing") | |
// Data size. // TODO: change this to be compatible with new double-style data | |
unsigned int dsize = 4; | |
// Loop over chunks in the ser file, copying each to the FID as we go along. | |
int j, k, *p; | |
for (j = 0; j < td1; j++) { // loop over chunks | |
// Read in TD2 points. This includes the group delay plus drop points. | |
if (fread(ser_array, dsize, td2, fp_ser) != td2) STOPMSG("min_pshift: error reading data from ser file") | |
// Reset the pointer to the start of ser_array. | |
p = ser_array; | |
// For the first chunk we want to copy over the group delay. | |
if (j == 0) { | |
if (fwrite(p, dsize, grpdly_points, fp_fid) != grpdly_points) STOPMSG("min_pshift: error writing to fid file") | |
} | |
// Advance the pointer beyond the group delay. If j > 0 then this simply | |
// skips over the group delay without copying it. | |
p += grpdly_points; | |
// Skip over the drop points. | |
p += drop_points; | |
// Write points_per_chunk points to the FID file. | |
if (fwrite(p, dsize, points_per_chunk, fp_fid) != points_per_chunk) STOPMSG("min_pshift: error writing to fid file") | |
} | |
// Close both file pointers. | |
fclose(fp_ser); | |
fclose(fp_fid); | |
// Open new dataset. | |
REXPNO(new_expno) | |
SETCURDATA | |
// Set true length of TD. | |
int new_td = (points_per_chunk * td1) + grpdly_points; | |
STOREPARS("TD", new_td) | |
STOREPAR("TD", new_td) | |
// Perform final processing. | |
XCMD("sendgui browse_update_tree") | |
EFP | |
APK | |
XCMD("abs n") | |
QUIT |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment