diff --git a/CMakeLists.txt b/CMakeLists.txt index caad358be..dd5477e70 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -382,6 +382,7 @@ set (wsjt_FSRCS lib/demod64a.f90 lib/determ.f90 lib/downsam9.f90 + lib/echo_snr.f90 lib/encode232.f90 lib/encode4.f90 lib/encode_msk40.f90 diff --git a/lib/echo_snr.f90 b/lib/echo_snr.f90 new file mode 100644 index 000000000..8a0ba4b84 --- /dev/null +++ b/lib/echo_snr.f90 @@ -0,0 +1,56 @@ +subroutine echo_snr(sa,sb,fspread0,blue,red,snrdb,db_err,fpeak,snr_detect) + + parameter (NZ=4096) + real sa(NZ) + real sb(NZ) + real blue(NZ) + real red(NZ) + integer ipkv(1) + equivalence (ipk,ipkv) + + df=12000.0/32768.0 + wh=0.5*fspread0+10.0 + i1=nint((1500.0 - 2.0*wh)/df) - 2048 + i2=nint((1500.0 - wh)/df) - 2048 + i3=nint((1500.0 + wh)/df) - 2048 + i4=nint((1500.0 + 2.0*wh)/df) - 2048 + + call pctile(sb(i1),i2-i1,50,r0) + call pctile(sb(i3+1),i4-i3,50,r1) + scale=(r1-r0)/(0.5*(i3+i4) - 0.5*(i1+i2)) + do i=1,NZ + x=i + y=r0 + scale*(i-0.5*(i1+i2)) + blue(i)=sa(i)/y + red(i)=sb(i)/y +! write(13,1100) i,i*df+750.0,sb(i),y,red(i) +!1100 format(i6,f10.2,2f10.2,f10.6) +enddo + +! baseline=(sum(sb(i1:i2-1)) + sum(sb(i3+1:i4)))/(i2+i4-i1-i3) +! blue=sa/baseline +! red=sb/baseline + psig=sum(red(i2:i3)-1.0) + pnoise_2500 = 2500.0/df + snrdb=db(psig/pnoise_2500) + + smax=0. + mh=max(1,nint(0.2*fspread0/df)) + do i=i2,i3 + ssum=sum(red(i-mh:i+mh)) + if(ssum.gt.smax) then + smax=ssum + ipk=i + endif + enddo + fpeak=ipk*df - 750.0 + + call averms(red(i1:i2-1),i2-i1,-1,ave1,rms1) + call averms(red(i3+1:i4),i4-i3,-1,ave2,rms2) + perr=0.707*(rms1+rms2)*sqrt(float(i2-i1+i4-i3)) + snr_detect=psig/perr + db_err=99.0 + if(psig.gt.perr) db_err=snrdb - db((psig-perr)/pnoise_2500) + + return +end subroutine echo_snr diff --git a/lib/test_snr.f90 b/lib/test_snr.f90 index d1ea1685f..9c61bdad6 100644 --- a/lib/test_snr.f90 +++ b/lib/test_snr.f90 @@ -35,13 +35,12 @@ program test_snr nfsample=h%nsamrate df=12000.0/NFFT s=0. + fac=1.0/NMAX do iping=1,npings read(10) id2(1:NMAX) - x(1:NMAX)=id2(1:NMAX) - x(NMAX+1:)=0. - rms=sqrt(dot_product(x(1:NMAX),x(1:NMAX))/NMAX) - + x(1:NMAX)=fac*id2(1:NMAX) + x(NMAX+1:)=0. call four2a(x,NFFT,1,-1,0) do i=1,8192 !Accumulate spectrum 0 - 3 kHz s(i)=s(i) + real(c(i))**2 + aimag(c(i))**2 @@ -53,14 +52,7 @@ program test_snr call echo_snr(sa,sb,fspread0,blue,red,snrdb,db_err,fpeak,snr_detect) - nsum=npings - nqual=0 - if(nsum.ge.2 .and. nsum.lt.4) nqual=(snr_detect-4)/5 - if(nsum.ge.4 .and. nsum.lt.8) nqual=(snr_detect-3)/4 - if(nsum.ge.8 .and. nsum.lt.12) nqual=(snr_detect-3)/3 - if(nsum.ge.12) nqual=(snr_detect-2.5)/2.5 - if(nqual.lt.0) nqual=0 - if(nqual.gt.10) nqual=10 + nqual=min(10,int(snr_detect-4.0)) write(*,1010) fspread0,snrdb0,snrdb,snrdb-snrdb0,db_err,fpeak, & snr_detect,nqual @@ -69,52 +61,7 @@ program test_snr do i=1,8192 write(12,1100) i*df,s(i) -1100 format(f10.3,e12.3) +1100 format(f10.3,e15.6) enddo 999 end program test_snr - -subroutine echo_snr(sa,sb,fspread0,blue,red,snrdb,db_err,fpeak,snr_detect) - - parameter (NZ=4096) - real sa(NZ) - real sb(NZ) - real blue(NZ) - real red(NZ) - integer ipkv(1) - equivalence (ipk,ipkv) - - df=12000.0/32768.0 - wh=0.5*fspread0+10.0 - i1=nint((1500.0 - 2.0*wh)/df) - 2048 - i2=nint((1500.0 - wh)/df) - 2048 - i3=nint((1500.0 + wh)/df) - 2048 - i4=nint((1500.0 + 2.0*wh)/df) - 2048 - - baseline=(sum(sb(i1:i2-1)) + sum(sb(i3+1:i4)))/(i2+i4-i1-i3) - blue=sa/baseline - red=sb/baseline - psig=sum(red(i2:i3)-1.0) - pnoise_2500 = 2500.0/df - snrdb=db(psig/pnoise_2500) - - smax=0. - mh=max(1,nint(0.2*fspread0/df)) - do i=i2,i3 - ssum=sum(red(i-mh:i+mh)) - if(ssum.gt.smax) then - smax=ssum - ipk=i - endif - enddo - fpeak=ipk*df - 750.0 - - call averms(red(i1:i2-1),i2-i1,-1,ave1,rms1) - call averms(red(i3+1:i4),i4-i3,-1,ave2,rms2) - perr=0.707*(rms1+rms2)*sqrt(float(i2-i1+i4-i3)) - snr_detect=psig/perr - db_err=99.0 - if(psig.gt.perr) db_err=snrdb - db((psig-perr)/pnoise_2500) - - return -end subroutine echo_snr