spacecruft
/
strf
1
0
Fork 0
Code Releases Activity

Updates

pull/10/head
Cees Bassa 2015-04-20 23:48:26 +02:00
parent 7b75ff732e
commit 3a7dd8e6e3
2 changed files with 177 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

31
frequencies.txt
View File

@ -358,3 +358,34 @@
25476 2264.857
40588 2215.993
23802 2264.978
27640 2215.876
27640 2216.900
32289 2216.817
27640 2215.364
40336 2213.126
29522 2262.196
29522 2264.859
27640 2221.508
40013 2220.005
39766 2220.239
38337 2224.003
35937 2237.258
28054 2235.463
24753 2237.498
29522 2239.570
35951 2239.567
24753 2239.567
39630 2239.567
29522 2238.545
35951 2238.545
24753 2238.545
39630 2238.545
29522 2236.497
35951 2236.497
24753 2236.497
39630 2236.497
28414 2239.846
35937 2238.283
35937 2239.307
37180 2239.844
39159 2238.851

146
rfplot.c
View File

@ -21,6 +21,9 @@ void time_axis(double *mjd,int n,float xmin,float xmax,float ymin,float ymax);
void usage(void);
void plot_traces(struct trace *t,int nsat);
struct trace fit_trace(struct spectrogram s,struct select sel,int site_id);
void convolve(float *y,int n,float *w,int m,float *z);
float gauss(float x,float w);
void quadfit(float x[],float y[],int n,float a[]);
// Fit trace
struct trace locate_trace(struct spectrogram s,struct trace t,int site_id)
@ -183,6 +186,79 @@ void filter(struct spectrogram s,int site_id)
return;
}
void peakfind(struct spectrogram s,int site_id,int i0,int i1,int j0,int j1)
{
int i,j,k,l,m=21,n;
float *w,*y,*sy,*a,*b,*c,d[3],dx[3],dw=2.0,x0,c0=-0.005;
double f;
FILE *file;
n=j1-j0;
// Allocate
y=(float *) malloc(sizeof(float)*n);
sy=(float *) malloc(sizeof(float)*n);
a=(float *) malloc(sizeof(float)*n);
b=(float *) malloc(sizeof(float)*n);
c=(float *) malloc(sizeof(float)*n);
// Make gaussian smoothing filter
w=(float *) malloc(sizeof(float)*m);
for (i=0;i<m;i++)
w[i]=gauss((float) (i-m/2),dw);
// Open file
file=fopen("peakfind.dat","w");
// Loop over subints
for (i=i0;i<i1;i++) {
if (s.mjd[i]==0.0)
continue;
// Fill array
for (j=0;j<n;j++)
y[j]=s.z[i+s.nsub*(j0+j)];
// Convolve
convolve(y,n,w,m,sy);
// Fit parabolas
dx[0]=-1.0;
dx[1]=0.0;
dx[2]=1.0;
for (j=1;j<n-1;j++) {
quadfit(dx,&sy[j-1],3,d);
a[j]=d[0];
b[j]=d[1];
c[j]=d[2];
}
// Mark points
for (j=0;j<n-1;j++) {
if (b[j]>0.0 && b[j+1]<0.0 && c[j]<c0) {
x0=(float) (j+j0)+b[j]/(b[j]-b[j+1]);
f=s.freq-0.5*s.samp_rate+(double) x0*s.samp_rate/(double) s.nchan;
if (s.mjd[i]>1.0)
fprintf(file,"%lf %lf %f %d\n",s.mjd[i],f,s.z[i+s.nsub*j],site_id);
cpgpt1((float) i+0.5,x0+0.5,17);
}
}
}
// Close
fclose(file);
// Free
free(y);
free(sy);
free(a);
free(b);
free(c);
free(w);
return;
}
int main(int argc,char *argv[])
{
struct spectrogram s;
@ -413,6 +489,22 @@ int main(int argc,char *argv[])
if (c=='g')
filter(s,site_id);
if (c=='G') {
i0=(int) floor(xmin);
i1=(int) ceil(xmax);
j0=(int) floor(ymin);
j1=(int) ceil(ymax);
if (i0<0)
i0=0;
if (i1>=s.nsub)
i1=s.nsub-1;
if (j0<0)
j0=0;
if (j1>=s.nchan)
j1=s.nchan-1;
peakfind(s,site_id,i0,i1,j0,j1);
}
// Fit
if (c=='f') {
tf=fit_trace(s,sel,site_id);
@ -985,3 +1077,57 @@ struct trace fit_trace(struct spectrogram s,struct select sel,int site_id)
return t;
}
// Discrete convolution
void convolve(float *y,int n,float *w,int m,float *z)
{
int i,j,k,imid;
imid=m/2;
for (i=0;i<n;i++) {
z[i]=0.0;
for (j=0;j<m;j++) {
k=i-imid+j;
if (k<0 || k>n-1)
continue;
z[i]+=w[j]*y[k];
}
}
return;
}
// Gaussian
float gauss(float x,float w)
{
float c;
c=pow(x/w,2);
return exp(-0.5*c)/(sqrt(2.0*M_PI)*w);
}
// Quadratic fit
void quadfit(float x[],float y[],int n,float a[])
{
int i;
float p,q,r,s,t,u,v,d;
p=q=r=s=t=u=v=0.;
for (i=0;i<n;i++) {
p+=x[i];
q+=pow(x[i],2);
r+=pow(x[i],3);
s+=pow(x[i],4);
t+=y[i];
u+=x[i]*y[i];
v+=pow(x[i],2)*y[i];
}
d=n*q*s+2.*p*q*r-q*q*q-p*p*s-(float) n*r*r;
a[0]=(q*s*t+q*r*u+p*r*v-q*q*v-p*s*u-r*r*t)/d;
a[1]=((float) n*s*u+p*q*v+q*r*t-q*q*u-p*s*t-(float) n*r*v)/d;
a[2]=((float) n*q*v+p*r*t+p*q*u-q*q*t-p*p*v-(float) n*r*u)/d;
return;
}