w4kek/WSJT-X
1
0
Fork 0
mirror of https://github.com/saitohirga/WSJT-X.git synced 2024-11-27 06:38:44 -05:00
Code Issues Packages Projects Releases Wiki Activity

Bug fixes for wsprd and wsprsim utilities. Implemented fully-coherent signal subtraction in wsprd_exp.

Browse Source 
git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@5621 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
This commit is contained in:
Steven Franke 2015-06-18 02:24:06 +00:00
parent 0864449f21
commit 8e4c20b8a5
5 changed files with 146 additions and 44 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

162
lib/wsprd/wsprd_exp.c
View File

@ -326,8 +326,10 @@ void sync_and_demodulate(double *id, double *qd, long np,
}
return;
}
//***************************************************************************
void subtract_signal(double *id, double *qd, long np,
/***************************************************************************
symbol-by-symbol signal subtraction
****************************************************************************/
void subtract_signal(double *id, double *qd, long np,
float f0, int shift0, float drift0, unsigned char* channel_symbols)
{
float dt=1.0/375.0, df=375.0/256.0;
@ -364,10 +366,12 @@ void subtract_signal(double *id, double *qd, long np,
}
}
// subtract the signal here.
i0=i0/256.0; //will be wrong for partial symbols at the edges...
q0=q0/256.0;
for (j=0; j<256; j++) {
k=shift0+i*256+j;
if( (k>0) & (k<np) ) {
@ -378,6 +382,92 @@ void subtract_signal(double *id, double *qd, long np,
}
return;
}
/******************************************************************************
Fully coherent signal subtraction
*******************************************************************************/
void subtract_signal2(double *id, double *qd, long np,
float f0, int shift0, float drift0, unsigned char* channel_symbols)
{
double dt=1.0/375.0, df=375.0/256.0;
int i, j, k, ii;
double pi=4.*atan(1.0),twopidt;
double refi[45000],refq[45000];
double ci[45000],cq[45000],cfi[45000],cfq[45000];
memset(refi,0,sizeof(double)*45000);
memset(refq,0,sizeof(double)*45000);
memset(ci,0,sizeof(double)*45000);
memset(cq,0,sizeof(double)*45000);
memset(cfi,0,sizeof(double)*45000);
memset(cfq,0,sizeof(double)*45000);
double phi=0, dphi;
// double dphi, cdphi, sdphi;
twopidt=2.0*pi*dt;
// measured signal is: s(t)=a(t)*exp( j*theta(t) )
// reference is: r(t) = exp( j*phi(t) )
// complex amplitude is estimated as: c(t)=LPF[s(t)*conjugate(r(t))]
// so c(t) has phase angle theta-phi
// multiply r(t) by c(t) and subtract from s(t), i.e. s'(t)=s(t)-c(t)r(t)
// create reference wspr signal vector, centered on f0.
for (i=0; i<162; i++) {
dphi=twopidt*
(
f0 + ((float)drift0/2.0)*((float)i-81.0)/81.0
+ ((double)channel_symbols[i]-1.5)*df
);
for ( j=0; j<256; j++ ) {
ii=256*i+j;
refi[ii]=refi[ii]+cos(phi); //cannot precompute sin/cos because dphi is changing
refq[ii]=refq[ii]+sin(phi);
phi=phi+dphi;
}
}
// s(t) * conjugate(r(t))
// place signal 1 impulse response width in so that we don't have to deal
// with partial convolutions at the beginning.
for (i=0; i<41472; i++) {
k=shift0+i;
if( (k>0) & (k<np) ) {
ci[i+500] = id[k]*refi[i] + qd[k]*refq[i];
cq[i+500] = qd[k]*refi[i] - id[k]*refq[i];
}
}
//quick and dirty filter - may want to do better
double w[500], norm=0;
for (i=0; i<500; i++) {
w[i]=sin(pi*i/499.0);
norm=norm+w[i];
}
for (i=0; i<500; i++) {
w[i]=w[i]/norm;
}
// LPF
for (i=500; i<45000-375; i++) {
cfi[i]=0.0; cfq[i]=0.0;
for (j=0; j<500; j++) {
cfi[i]=cfi[i]+w[j]*ci[i-250+j];
cfq[i]=cfq[i]+w[j]*cq[i-250+j];
}
}
// subtract c(t)*ref(i) here
// (ci+j*cq)(refi+j*refq)=(ci*refi-cq*refq)+j(ci*refq)+cq*refi)
for (i=0; i<41472; i++) {
k=shift0+i;
if( (k>0) & (k<np) ) {
id[k]=id[k] - (cfi[i+500]*refi[i]-cfq[i+500]*refq[i]);
qd[k]=qd[k] - (cfi[i+500]*refq[i]+cfq[i+500]*refi[i]);
}
}
return;
}
unsigned long writec2file(char *c2filename, int trmin, double freq
, double *idat, double *qdat)
@ -417,6 +507,7 @@ void usage(void)
printf("\n");
printf("Options:\n");
printf(" -a <path> path to writeable data files, default=\".\"\n");
printf(" -c write .c2 file at the end of the first pass\n");
printf(" -e x (x is transceiver dial frequency error in Hz)\n");
printf(" -f x (x is transceiver dial frequency in MHz)\n");
// blanking is not yet implemented. The options are accepted for compatibility
@ -425,7 +516,6 @@ void usage(void)
// printf(" -b n (n is pct of time that is blanked)\n");
printf(" -H do not use (or update) the hash table\n");
printf(" -m decode wspr-15 .wav file\n");
printf(" -n write noise estimates to file noise.dat\n");
printf(" -q quick mode - doesn't dig deep for weak signals\n");
printf(" -s signal subtraction mode\n");
printf(" -t signal subtraction followed by a second pass\n");
@ -449,7 +539,8 @@ int main(int argc, char *argv[])
char timer_fname[200],hash_fname[200];
char uttime[5],date[7];
int c,delta,maxpts=65536,verbose=0,quickmode=0;
int writenoise=0,usehashtable=1,wspr_type=2, subtraction=0, ipass, npasses=1;
int writenoise=0,usehashtable=1,wspr_type=2, subtraction=0, ipass;
int writec2=0, npasses=1;
int shift1, lagmin, lagmax, lagstep, worth_a_try, not_decoded;
unsigned int nbits=81;
unsigned int npoints, metric, maxcycles, cycles, maxnp;
@ -478,7 +569,7 @@ int main(int argc, char *argv[])
char allcalls[100][13];
memset(allfreqs,0,sizeof(float)*100);
memset(allcalls,0,sizeof(char)*100*13);
int uniques=0, noprint=0;
// Parameters used for performance-tuning:
@ -501,11 +592,14 @@ int main(int argc, char *argv[])
idat=malloc(sizeof(double)*maxpts);
qdat=malloc(sizeof(double)*maxpts);
while ( (c = getopt(argc, argv, "a:b:e:f:Hmnqstwvz:")) !=-1 ) {
while ( (c = getopt(argc, argv, "a:ce:f:Hmqstwvz:")) !=-1 ) {
switch (c) {
case 'a':
data_dir = optarg;
break;
case 'c':
writec2=1;
break;
case 'e':
dialfreq_error = strtof(optarg,NULL); // units of Hz
// dialfreq_error = dial reading - actual, correct frequency
@ -519,17 +613,14 @@ int main(int argc, char *argv[])
case 'm':
wspr_type = 15;
break;
case 'n':
writenoise = 1;
break;
case 'q':
quickmode = 1;
break;
case 's':
subtraction = 1;
subtraction = 1; //subtraction only, no second pass
break;
case 't':
subtraction = 1;
subtraction = 1; //subtraction and a second pass (t for two)
npasses = 2; //npasses defaults to 1
break;
case 'v':
@ -654,6 +745,14 @@ int main(int argc, char *argv[])
//*************** main loop starts here *****************
for (ipass=0; ipass<npasses; ipass++) {
if( ipass == 1 && uniques == 0 ) break;
if( ipass == 1 ) { //otherwise we bog down on the second pass
minsync1=0.18;
minsync2=0.2;
maxcycles=5000;
printf("-------------- 2 ----------------\n");
}
memset(ps,0.0, sizeof(float)*512*nffts);
for (i=0; i<nffts; i++) {
@ -750,8 +849,7 @@ int main(int argc, char *argv[])
}
}
npk=i;
/*
// bubble sort on snr, bringing freq along for the ride
int pass;
float tmp;
@ -767,7 +865,7 @@ int main(int argc, char *argv[])
}
}
}
*/
t0=clock();
/* Make coarse estimates of shift (DT), freq, and drift
@ -853,7 +951,7 @@ int main(int argc, char *argv[])
memset(symbols,0,sizeof(char)*nbits*2);
memset(callsign,0,sizeof(char)*13);
memset(call_loc_pow,0,sizeof(char)*23);
f1=freq0[j];
drift1=drift0[j];
shift1=shift0[j];
@ -906,11 +1004,11 @@ int main(int argc, char *argv[])
sq += y*y;
}
rms=sqrt(sq/162.0);
if((sync1 > minsync2) && (rms > minrms)) {
deinterleave(symbols);
t0 = clock();
not_decoded = fano(&metric,&cycles,&maxnp,decdata,symbols,nbits,
mettab,delta,maxcycles);
tfano += (double)(clock()-t0)/CLOCKS_PER_SEC;
@ -927,9 +1025,9 @@ int main(int argc, char *argv[])
}
if( worth_a_try && !not_decoded ) {
for(i=0; i<11; i++) {
if( decdata[i]>127 ) {
message[i]=decdata[i]-256;
} else {
@ -943,12 +1041,12 @@ int main(int argc, char *argv[])
// call_loc_pow string and also callsign (for de-duping).
noprint=unpk_(message,hashtab,call_loc_pow,callsign);
if( subtraction && !noprint ) {
unsigned char channel_symbols[162];
if( subtraction && (ipass == 0) && !noprint ) {
unsigned char channel_symbols[162];
if( get_wspr_channel_symbols(call_loc_pow, channel_symbols) ) {
subtract_signal(idat, qdat, npoints, f1, shift1, drift1, channel_symbols);
subtract_signal2(idat, qdat, npoints, f1, shift1, drift1, channel_symbols);
} else {
break;
}
@ -965,7 +1063,7 @@ int main(int argc, char *argv[])
strcpy(allcalls[uniques],callsign);
allfreqs[uniques]=f1;
uniques++;
// Add an extra space at the end of each line so that wspr-x doesn't
// truncate the power (TNX to DL8FCL!)
@ -1002,6 +1100,15 @@ int main(int argc, char *argv[])
}
}
}
if( ipass == 0 && writec2 ) {
char c2filename[15];
double carrierfreq=dialfreq;
int wsprtype=2;
strcpy(c2filename,"000000_0001.c2");
printf("Writing %s\n",c2filename);
writec2file(c2filename, wsprtype, carrierfreq, idat, qdat);
}
}
printf("<DecodeFinished>\n");
@ -1048,13 +1155,6 @@ int main(int argc, char *argv[])
fclose(fhash);
}
// char c2filename[15];
// double carrierfreq=dialfreq;
// int wsprtype=2;
// strcpy(c2filename,"000000_0001.c2");
// printf("Writing %s\n",c2filename);
// writec2file(c2filename, wsprtype, carrierfreq, idat, qdat);
if(writenoise == 999) return -1; //Silence compiler warning
return 0;
}

13
lib/wsprd/wsprd_utils.c
View File

@ -66,7 +66,7 @@ void unpack50( signed char *dat, int32_t *n1, int32_t *n2 )
*n2=*n2+((i4>>6)&3);
}
void unpackcall( int32_t ncall, char *call )
int unpackcall( int32_t ncall, char *call )
{
char c[]={'0','1','2','3','4','5','6','7','8','9','A','B','C','D','E',
'F','G','H','I','J','K','L','M','N','O','P','Q','R','S','T',
@ -108,10 +108,13 @@ void unpackcall( int32_t ncall, char *call )
call[i]='\0';
}
}
} else {
return 0;
}
return 1;
}
void unpackgrid( int32_t ngrid, char *grid)
int unpackgrid( int32_t ngrid, char *grid)
{
char c[]={'0','1','2','3','4','5','6','7','8','9','A','B','C','D','E',
'F','G','H','I','J','K','L','M','N','O','P','Q','R','S','T',
@ -139,7 +142,9 @@ void unpackgrid( int32_t ngrid, char *grid)
grid[3]=c[n2];
} else {
strcpy(grid,"XXXX");
return 0;
}
return 1;
}
int unpackpfx( int32_t nprefix, char *call)
@ -233,8 +238,8 @@ int unpk_(signed char *message, char hashtab[32768][13], char *call_loc_pow, cha
char grid[5],grid6[7],cdbm[3];
unpack50(message,&n1,&n2);
unpackcall(n1,callsign);
unpackgrid(n2, grid);
if( !unpackcall(n1,callsign) ) return 1;
if( !unpackgrid(n2, grid) ) return 1;
int ntype = (n2&127) - 64;
callsign[12]=0;
grid[4]=0;

4
lib/wsprd/wsprd_utils.h
View File

@ -12,9 +12,9 @@
void unpack50( signed char *dat, int32_t *n1, int32_t *n2 );
void unpackcall( int32_t ncall, char *call );
int unpackcall( int32_t ncall, char *call );
void unpackgrid( int32_t ngrid, char *grid);
int unpackgrid( int32_t ngrid, char *grid);
int unpackpfx( int32_t nprefix, char *call);

1
lib/wsprd/wsprsim.c
View File

@ -133,7 +133,6 @@ int main(int argc, char *argv[])
message=malloc(sizeof(char)*23);
strcpy(c2filename,"000000_0001.c2");
printf("%s\n",c2filename);
while ( (c = getopt(argc, argv, "cdo:s:")) !=-1 ) {
switch (c) {
case 'c':

10
lib/wsprd/wsprsim_utils.c
View File

@ -84,7 +84,6 @@ void pack_prefix(char *callsign, int32_t *n, int32_t *m, int32_t *nadd ) {
char *call6;
call6=malloc(sizeof(char)*6);
memset(call6,32,sizeof(char)*6);
int i1=strcspn(callsign,"/");
if( callsign[i1+2] == 0 ) {
@ -200,10 +199,9 @@ int get_wspr_channel_symbols(char* rawmessage, unsigned char* symbols) {
if( (i1>3) & (i1<7) & (i2==mlen) & (i3==mlen) ) {
// Type 1 message: K9AN EN50 33
// xxnxxxx xxnn nn
const char s[2]=" ";
callsign = strtok(message,s);
grid = strtok(NULL,s);
powstr = strtok(NULL,s);
callsign = strtok(message," ");
grid = strtok(NULL," ");
powstr = strtok(NULL," ");
int power = atoi(powstr);
n = pack_call(callsign);
@ -239,6 +237,7 @@ int get_wspr_channel_symbols(char* rawmessage, unsigned char* symbols) {
} else if ( i2 < mlen ) { // just looks for a right slash
// Type 2: PJ4/K1ABC 37
callsign=strtok(message," ");
if( strlen(callsign) < i2 ) return 0; //guards against pathological case
powstr=strtok(NULL," ");
int power = atoi(powstr);
if( power < 0 ) power=0;
@ -250,7 +249,6 @@ int get_wspr_channel_symbols(char* rawmessage, unsigned char* symbols) {
m=128*ng+ntype+64;
n=n1;
} else {
// printf("Error: bad message format\n");
return 0;
}