From df8b038b8f60f26ae5d8797d504a24986f656457 Mon Sep 17 00:00:00 2001 From: "cbassa@gmail.com" Date: Sun, 22 Nov 2020 14:04:04 +0100 Subject: [PATCH 1/8] Wideband point fitting --- makefile | 4 +- rfplot.c | 227 +++++++++++++++++++++++++++++++++++++++++++++++++++++-- 2 files changed, 221 insertions(+), 10 deletions(-) diff --git a/makefile b/makefile index a3d2003..87f874d 100644 --- a/makefile +++ b/makefile @@ -34,8 +34,8 @@ rffind: rffind.o rfio.o rftime.o rftrack: rftrack.o rfio.o rftime.o rftrace.o sgdp4.o satutl.o deep.o ferror.o $(CC) -o rftrack rftrack.o rfio.o rftime.o rftrace.o sgdp4.o satutl.o deep.o ferror.o -lm -rfplot: rfplot.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o - gfortran -o rfplot rfplot.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o $(LFLAGS) +rfplot: rfplot.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o versafit.o dsmin.o simplex.o + gfortran -o rfplot rfplot.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o versafit.o dsmin.o simplex.o $(LFLAGS) rffft: rffft.o rftime.o $(CC) -o rffft rffft.o rftime.o -lfftw3f -lm diff --git a/rfplot.c b/rfplot.c index d4f02f7..f5e9cea 100644 --- a/rfplot.c +++ b/rfplot.c @@ -17,6 +17,10 @@ struct select { float x[NMAX],y[NMAX],w; float sigma; }; +struct data { + int n; + double *x, *y, *ym; +} gd; void dec2sex(double x,char *s,int f,int len); void time_axis(double *mjd,int n,float xmin,float xmax,float ymin,float ymax); @@ -30,6 +34,9 @@ void quadfit(float x[],float y[],int n,float a[]); struct trace locate_trace(struct spectrogram s,struct trace t,int site_id,float sigmamin,float wmin,float wmax,int graves); void filter(struct spectrogram s,int site_id,float sigma,int graves); void peakfind(struct spectrogram s,int site_id,int i0,int i1,int j0,int j1); +void versafit(int m,int n,double *a,double *da,double (*func)(double *),double dchisq,double tol,char *opt); +double chisq_gaussian(double a[]); +float fit_gaussian_point(struct spectrogram s,float x,float y,struct select sel,int site_id,int graves); int main(int argc,char *argv[]) { @@ -54,7 +61,7 @@ int main(int argc,char *argv[]) int ix=0,iy=0,isub=0; int i0,j0,i1,j1,jmax; float width=1500; - float x,y,x0,y0; + float x,y,x0,y0,yfit; char c; char path[128],xlabel[128],ylabel[64],filename[32],tlefile[128]; int sec,lsec,ssec; @@ -83,9 +90,15 @@ int main(int argc,char *argv[]) env=getenv("ST_TLEDIR"); sprintf(tlefile,"%s/bulk.tle",env); + // Set selection + isel=0; + sel.n=0; + sel.w=100.0; + sel.sigma=5.0; + // Read arguments if (argc>1) { - while ((arg=getopt(argc,argv,"p:f:w:s:l:b:z:hc:C:gm:o:"))!=-1) { + while ((arg=getopt(argc,argv,"p:f:w:s:l:b:z:hc:C:gm:o:S:W:"))!=-1) { switch (arg) { case 'p': @@ -111,6 +124,14 @@ int main(int argc,char *argv[]) case 'w': df0=(double) atof(optarg); break; + + case 'W': + sel.w=atof(optarg); + break; + + case 'S': + sel.sigma=atof(optarg); + break; case 'z': zmax=atof(optarg); @@ -180,12 +201,6 @@ int main(int argc,char *argv[]) // Set trace tf.n=0; - // Set selection - isel=0; - sel.n=0; - sel.w=100.0; - sel.sigma=5.0; - // Forever loop for (;;) { if (redraw==1) { @@ -315,7 +330,9 @@ int main(int argc,char *argv[]) printf("t Locate traces\n"); printf("s Select track points\n"); printf("u Undo track point selection\n"); + printf("S Set track selection width and sigma\n"); printf("f Fits selected track points (generates out.dat)\n"); + printf("F Fits selected wideband track points (generates out.dat)\n"); printf("g Filter points above a certain sigma limit (generates filter.dat)\n"); printf("G Find peaks in current window (generates peakfind.dat)\n"); printf("i Identify selected track points\n"); @@ -324,6 +341,7 @@ int main(int argc,char *argv[]) printf("b Increase dynamic range\n"); printf("Z Provide manual dynamic range limits\n"); printf("D/mid Mark point (appends to mark.dat)\n"); + printf("w Fit wideband point (appends to mark.dat)\n"); printf("c Center view\n"); printf("C Toggle through color maps\n"); printf("p/right Toggle overlays\n"); @@ -400,6 +418,13 @@ int main(int argc,char *argv[]) continue; } + // Fit + if (c=='F') { + tf=fit_gaussian_trace(s,sel,site_id,graves); + tf.site=site_id; + continue; + } + // Identify if (c=='i') { if (graves==0) @@ -530,6 +555,26 @@ int main(int argc,char *argv[]) fclose(file); } + // Fit single wideband point + if (c=='w') { + file=fopen("mark.dat","a"); + + // Fit point + yfit=fit_gaussian_point(s,x,y,sel,site_id,graves); + cpgpt1(x,yfit,17); + i=(int) floor(x); + f=s.freq-0.5*s.samp_rate+(double) yfit*s.samp_rate/(double) s.nchan; + if (s.mjd[i]>1.0) { + if (graves==0) + fprintf(file,"%lf %lf %f %d\n",s.mjd[i],f,s.z[i+s.nsub*j],site_id); + else + fprintf(file,"%lf %lf %f %d 9999\n",s.mjd[i],f,s.z[i+s.nsub*j],site_id); + printf("%lf %lf %f %d\n",s.mjd[i],f,s.z[i+s.nsub*j],site_id); + } + fclose(file); + } + + // Mark if (c=='m') { i0=(int) floor(xmin); @@ -732,6 +777,9 @@ int main(int argc,char *argv[]) free(t[i].freq); free(t[i].za); } + free(gd.x); + free(gd.y); + free(gd.ym); return 0; } @@ -849,6 +897,8 @@ void usage(void) printf("-f Frequency to zoom into (Hz)\n"); printf("-w Bandwidth to zoom into (Hz)\n"); printf("-o Frequency offset to apply (Hz) [0]\n"); + printf("-W Track selection width (in pixels) [100]\n"); + printf("-S Track selection significance [5]\n"); printf("-C Site ID\n"); printf("-c TLE catalog\n"); printf("-g GRAVES data\n"); @@ -985,6 +1035,99 @@ struct trace fit_trace(struct spectrogram s,struct select sel,int site_id,int gr return t; } +// Fit gaussian trace +struct trace fit_gaussian_trace(struct spectrogram s,struct select sel,int site_id,int graves) +{ + int i,j,k,l,sn; + int i0,i1,j0,j1,jmax; + double f,chi2; + double a[4],da[4]; + float x,y,s1,s2,z,za,zs,zm,sigma,rms; + struct trace t; + FILE *file; + + // Set up trace + t.satno=99999; + t.n=(int) ceil(sel.x[sel.n-1]-sel.x[0]); + t.mjd=(double *) malloc(sizeof(double)*t.n); + t.freq=(double *) malloc(sizeof(double)*t.n); + t.za=(float *) malloc(sizeof(float)*t.n); + + // Set up data + gd.n=(int) 2 * sel.w; + gd.x=(double *) malloc(sizeof(double) * gd.n); + gd.y=(double *) malloc(sizeof(double) * gd.n); + gd.ym=(double *) malloc(sizeof(double) * gd.n); + + // Open file + file=fopen("out.dat","w"); + + // Loop over selected regions + for (k=0,l=0;k=s.nchan) + j1=s.nchan; + + // Loop over points + for (j=j0,l=0;jsel.sigma && s.mjd[i]>1.0) { + f=s.freq-0.5*s.samp_rate+(double) a[0]*s.samp_rate/(double) s.nchan; + if (graves==0) + fprintf(file,"%lf %lf %f %d\n",s.mjd[i],f,sigma,site_id); + else + fprintf(file,"%lf %lf %f %d 9999\n",s.mjd[i],f,sigma,site_id); + cpgpt1((float) i,(float) a[0], 17); + cpgmove((float) i,(float) a[0]-a[1]); + cpgdraw((float) i,(float) a[0]+a[1]); + t.mjd[l]=s.mjd[i]; + t.freq[l]=f; + t.za[l]=0.0; + l++; + } + } + } + t.n=l; + + // Close file + fclose(file); + + return t; +} + // Discrete convolution void convolve(float *y,int n,float *w,int m,float *z) { @@ -1285,3 +1428,71 @@ void peakfind(struct spectrogram s,int site_id,int i0,int i1,int j0,int j1) return; } + +// Gaussian point +double chisq_gaussian(double a[]) +{ + int i; + double arg,amp,chisq; + + // Compute chisq + for (i=0,chisq=0.0;i=s.nchan) + j1=s.nchan; + + // Loop over points + for (j=j0,l=0;j Date: Sat, 28 Nov 2020 20:33:05 +0100 Subject: [PATCH 2/8] Add user output filename option --- rffft.c | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) diff --git a/rffft.c b/rffft.c index 823e734..77d69c7 100644 --- a/rffft.c +++ b/rffft.c @@ -13,6 +13,7 @@ void usage(void) printf("rffft: FFT RF observations\n\n"); printf("-i Input file (can be fifo) [stdin]\n"); printf("-p Output prefix\n"); + printf("-o Output filename [default: YYYY-MM-DDTHH:MM:SS.sss_XXXXXX.bin]\n"); printf("-f Center frequency (Hz)\n"); printf("-s Sample rate (Hz)\n"); printf("-c Channel size [100Hz]\n"); @@ -31,11 +32,11 @@ void usage(void) int main(int argc,char *argv[]) { - int i,j,k,l,m,nchan,nint=1,arg=0,nbytes,nsub=60,flag,nuse=1,realtime=1,quiet=0,imin,imax,partial=0; + int i,j,k,l,m,nchan,nint=1,arg=0,nbytes,nsub=60,flag,nuse=1,realtime=1,quiet=0,imin,imax,partial=0,useoutput=0; fftwf_complex *c,*d; fftwf_plan fft; FILE *infile,*outfile; - char infname[128]="",outfname[128]="",path[64]=".",prefix[32]=""; + char infname[128]="",outfname[128]="",path[64]=".",prefix[32]="",output[128]=""; char informat='i',outformat='f'; int16_t *ibuf; char *cbuf; @@ -48,7 +49,7 @@ int main(int argc,char *argv[]) // Read arguments if (argc>1) { - while ((arg=getopt(argc,argv,"i:f:s:c:t:p:n:hm:F:T:bqR:"))!=-1) { + while ((arg=getopt(argc,argv,"i:f:s:c:t:p:n:hm:F:T:bqR:o:"))!=-1) { switch(arg) { case 'i': @@ -59,6 +60,11 @@ int main(int argc,char *argv[]) strcpy(path,optarg); break; + case 'o': + strcpy(output,optarg); + useoutput=1; + break; + case 'f': freq=(double) atof(optarg); break; @@ -191,7 +197,11 @@ int main(int argc,char *argv[]) // Forever loop for (m=0;;m++) { // File name - sprintf(outfname,"%s/%s_%06d.bin",path,prefix,m); + if (useoutput==0) { + sprintf(outfname,"%s/%s_%06d.bin",path,prefix,m); + } else { + sprintf(outfname,"%s/%s_%06d.bin",path,output,m); + } outfile=fopen(outfname,"w"); // Loop over subints to dump From 0dc92e2d312bc60cb1ffed42c231fccc84010690 Mon Sep 17 00:00:00 2001 From: Jan - PE0SAT Date: Thu, 8 Apr 2021 15:36:46 +0200 Subject: [PATCH 3/8] Update Makefile.osx Changed this Makefile so it can be used with more recent OSX/Macports solutions. --- Makefile.osx | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/Makefile.osx b/Makefile.osx index bfbb04d..9c38776 100644 --- a/Makefile.osx +++ b/Makefile.osx @@ -19,7 +19,8 @@ LFLAGS = -L$(prefix)/lib -lcpgplot -lpgplot -lX11 -lpng -lm -lgsl -lgslcblas # NOTE: STRF will not compile or link correctly with the system gcc (which is actually clang) # It's best to build with gcc provided by Macports or Homebrew # Under Macports, this is provided as gcc-mp-7, as below, on Homebrew this may be different. -CC = $(prefix)/bin/gcc-mp-7 +# Newer Macports releases use gcc-mp-10, therefor changed the default to gcc-mp-10 +CC = $(prefix)/bin/gcc-mp-10 # Installation INSTALL_PROGRAM = install -m 557 @@ -44,8 +45,8 @@ rffind: rffind.o rfio.o rftime.o rftrack: rftrack.o rfio.o rftime.o rftrace.o sgdp4.o satutl.o deep.o ferror.o $(CC) -o rftrack rftrack.o rfio.o rftime.o rftrace.o sgdp4.o satutl.o deep.o ferror.o -lm -rfplot: rfplot.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o - $(CC) -o rfplot rfplot.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o $(LFLAGS) +rfplot: rfplot.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o versafit.o dsmin.o simplex.o + $(CC) -o rfplot rfplot.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o versafit.o dsmin.o simplex.o $(LFLAGS) rffft: rffft.o rftime.o $(CC) -o rffft rffft.o rftime.o -lfftw3f -lm $(LFLAGS) From 0bd116be5a315f9126b24157291b3820c4e5e4a0 Mon Sep 17 00:00:00 2001 From: Cees Bassa Date: Thu, 8 Apr 2021 17:27:05 +0200 Subject: [PATCH 4/8] Add option to allow band inversion --- rffft.c | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/rffft.c b/rffft.c index 77d69c7..67a0e3a 100644 --- a/rffft.c +++ b/rffft.c @@ -23,6 +23,7 @@ void usage(void) printf("-F Input format char, int, float [int]\n"); printf("-T YYYY-MM-DDTHH:MM:SSS.sss\n"); printf("-R Frequency range to store (Hz)\n"); + printf("-I Invert frequencies\n"); printf("-b Digitize output to bytes [off]\n"); printf("-q Quiet mode, no output [off]\n"); printf("-h This help\n"); @@ -46,10 +47,11 @@ int main(int argc,char *argv[]) double freq,samp_rate,mjd,freqmin=-1,freqmax=-1; struct timeval start,end; char tbuf[30],nfd[32],header[256]=""; + int sign=1; // Read arguments if (argc>1) { - while ((arg=getopt(argc,argv,"i:f:s:c:t:p:n:hm:F:T:bqR:o:"))!=-1) { + while ((arg=getopt(argc,argv,"i:f:s:c:t:p:n:hm:F:T:bqR:o:I"))!=-1) { switch(arg) { case 'i': @@ -115,6 +117,10 @@ int main(int argc,char *argv[]) realtime=0; break; + case 'I': + sign=-1; + break; + case 'h': usage(); return 0; @@ -235,17 +241,17 @@ int main(int argc,char *argv[]) if (informat=='i') { for (i=0;i Date: Thu, 8 Apr 2021 17:27:53 +0200 Subject: [PATCH 5/8] Change coincidence algorithm and output frequency precision --- rffit.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/rffit.c b/rffit.c index be57894..ef77c85 100644 --- a/rffit.c +++ b/rffit.c @@ -257,7 +257,7 @@ int identify_satellite_from_visibility(char *catalog,double altmin) } } frac=(float) nalt/(float) nsel; - if (nsel==nalt) + if (frac>0.95) printf("%5d %d/%d %.4f\n",orb.satno,nalt,nsel,frac); } rewind(fp); @@ -727,7 +727,7 @@ int main(int argc,char *argv[]) } else { rms=fit_curve(orb,ia); //printf("rms: %.3f kHz, cf: %.6f MHz, TCA: %s\n",rms,d.ffit/1000.0,nfd); - printf("%05d %.3f %.3f %s %04d %02d%010.6f\n",orb.satno,d.ffit/1000.0,rms,nfd,d.p[0].site_id,orb.ep_year-2000,orb.ep_day); + printf("%05d %.6f %.3f %s %04d %02d%010.6f\n",orb.satno,d.ffit/1000.0,rms,nfd,d.p[0].site_id,orb.ep_year-2000,orb.ep_day); redraw=1; } } From 6cd277083de8015dfbca7c45c408f00a91ce8f48 Mon Sep 17 00:00:00 2001 From: Cees Bassa Date: Thu, 8 Apr 2021 17:34:41 +0200 Subject: [PATCH 6/8] Updated gitignore --- .gitignore | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/.gitignore b/.gitignore index 9fa3b1b..8ed5565 100644 --- a/.gitignore +++ b/.gitignore @@ -16,3 +16,12 @@ *.exe *.out *.app + +/rfdop +/rfedit +/rffft +/rffind +/rffit +/rfplot +/rfpng +/rfinfo \ No newline at end of file From 3c49fdc8311825310e170bc0dea39a55cc51da73 Mon Sep 17 00:00:00 2001 From: Cees Bassa Date: Thu, 8 Apr 2021 17:36:00 +0200 Subject: [PATCH 7/8] Free memory --- rffft.c | 2 ++ 1 file changed, 2 insertions(+) diff --git a/rffft.c b/rffft.c index 67a0e3a..d2127f4 100644 --- a/rffft.c +++ b/rffft.c @@ -360,6 +360,8 @@ int main(int argc,char *argv[]) // Deallocate free(ibuf); + free(cbuf); + free(fbuf); fftwf_free(c); fftwf_free(d); free(z); From 0177a2d38f3d88408bc2961c53cf4e41002b5238 Mon Sep 17 00:00:00 2001 From: Cees Bassa Date: Thu, 8 Apr 2021 17:39:25 +0200 Subject: [PATCH 8/8] Add site --- rffind.c | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/rffind.c b/rffind.c index db945d1..8b035e4 100644 --- a/rffind.c +++ b/rffind.c @@ -117,7 +117,6 @@ void usage(void) printf("-w Bandwidth to zoom into (Hz)\n"); printf("-o Frequency offset to apply\n"); printf("-C Site ID\n"); - printf("-c TLE catalog\n"); printf("-g GRAVES data\n"); printf("-S Sigma limit [default: 5.0]\n"); printf("-h This help\n"); @@ -148,7 +147,7 @@ int main(int argc,char *argv[]) // Read arguments if (argc>1) { - while ((arg=getopt(argc,argv,"p:f:w:s:l:hc:o:S:g"))!=-1) { + while ((arg=getopt(argc,argv,"p:f:w:s:l:hC:o:S:g"))!=-1) { switch (arg) { case 'p': @@ -158,6 +157,10 @@ int main(int argc,char *argv[]) case 's': isub=atoi(optarg); break; + + case 'C': + site_id=atoi(optarg); + break; case 'l': nsub=atoi(optarg);