From a9fc41b0ee9aa443629106b43faf0199a9ddba9e Mon Sep 17 00:00:00 2001 From: Joe Taylor Date: Thu, 27 Dec 2007 20:16:42 +0000 Subject: [PATCH] Histogram to get median in spec.F90, instead of pctile. git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/map65@606 ab8295b8-cf94-4d9e-aec4-7959e3be5d79 --- map65.py | 2 +- spec.f90 | 42 ++++++++++++++++++++++++++++++++++++++++-- 2 files changed, 41 insertions(+), 3 deletions(-) diff --git a/map65.py b/map65.py index f90742581..631af4923 100644 --- a/map65.py +++ b/map65.py @@ -1,4 +1,4 @@ -#---------------------------------------------------------------------- MAP65 +#--------------------------------------------------------------------- MAP65 # $Date$ $Revision$ # from Tkinter import * diff --git a/spec.f90 b/spec.f90 index af56dfe12..40f399453 100644 --- a/spec.f90 +++ b/spec.f90 @@ -15,9 +15,10 @@ subroutine spec(brightness,contrast,ngain,nspeed,a,a2) integer nstep(5) integer b0,c0 + integer hist(0:1000) ! Could save memory by doing the averaging-by-7 (or 10?) of ss5 in symspec. include 'spcom.f90' - real s(NFFT,NY) + real s(NFFT,NY),savg2(NFFT) include 'gcom1.f90' include 'gcom2.f90' include 'gcom3.f90' @@ -45,6 +46,42 @@ subroutine spec(brightness,contrast,ngain,nspeed,a,a2) enddo enddo enddo + call zero(savg2,NFFT) + do j=1,nlines + do i=1,NFFT + savg2(i)=savg2(i) + s(i,j) + enddo + enddo + + ia=0.08*NFFT + ib=0.92*NFFT + smin=1.e30 + smax=-smin + sum=0. + nsum=0 + do i=ia,ib + smin=min(savg2(i),smin) + smax=max(savg2(i),smax) + if(savg2(i).lt.10000.0) then + sum=sum + savg2(i) + nsum=nsum+1 + endif + enddo + ave=sum/nsum + call zero(hist,1001) + do i=ia,ib + n=savg2(i) * (300.0/ave) + if(n.gt.1000) n=1000 + hist(n)=hist(n)+1 + enddo + + sum=0. + do i=0,1000 + sum=sum + float(hist(i))/(ib-ia+1) + if(sum.gt.0.4) go to 10 + enddo +10 base=i*ave/300.0 + base=base/(nadd*nlines) newpts=NX*nlines do i=newpts+1,NX*NY @@ -56,6 +93,8 @@ subroutine spec(brightness,contrast,ngain,nspeed,a,a2) gamma=1.3 + 0.01*contrast offset=(brightness+64.0)/2 fac=20.0/nadd + fac=fac*0.065/base + ! fac=fac*(0.1537/base) foffset=0.001*(1270+nfcal) nbpp=(nfb-nfa)*NFFT/(96.0*NX) !Bins per pixel in wideband (upper) waterfall fselect=mousefqso + foffset @@ -73,7 +112,6 @@ subroutine spec(brightness,contrast,ngain,nspeed,a,a2) do j=nlines,1,-1 !Reverse order so last will be on top do i=1,NX k=k+1 - n=0 x=0. iia=(i-1)*nbpp + ii0 + 1