WSJT-X/flatten.f
Joe Taylor 06291ab964 First good compile with map65a running inside map65.
git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/map65@334 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
2007-01-11 21:25:52 +00:00

106 lines
3.1 KiB
Fortran

subroutine flatten(s2,nbins,jz,psa,ref,birdie,variance)
C Examines the 2-d spectrum s2(nbins,jz) and makes a reference spectrum
C from the jz/2 spectra below the 50th percentile in total power. Uses
C reference spectrum (with birdies removed) to flatten the passband.
real s2(nbins,jz) !2d spectrum
real psa(nbins) !Grand average spectrum
real ref(nbins) !Ref spect: smoothed ave of lower half
real birdie(nbins) !Spec (with birdies) for plot, in dB
real variance(nbins)
real ref2(750) !Work array
real power(300)
C Find power in each time block, then get median
do j=1,jz
s=0.
do i=1,nbins
s=s+s2(i,j)
enddo
power(j)=s
enddo
call pctile(power,ref2,jz,50,xmedian)
if(jz.lt.5) go to 900
C Get variance in each freq channel, using only those spectra with
C power below the median.
do i=1,nbins
s=0.
nsum=0
do j=1,jz
if(power(j).le.xmedian) then
s=s+s2(i,j)
nsum=nsum+1
endif
enddo
s=s/nsum
sq=0.
do j=1,jz
if(power(j).le.xmedian) sq=sq + (s2(i,j)/s-1.0)**2
enddo
variance(i)=sq/nsum
enddo
C Get grand average, and average of spectra with power below median.
call zero(psa,nbins)
call zero(ref,nbins)
nsum=0
do j=1,jz
call add(psa,s2(1,j),psa,nbins)
if(power(j).le.xmedian) then
call add(ref,s2(1,j),ref,nbins)
nsum=nsum+1
endif
enddo
do i=1,nbins !Normalize the averages
psa(i)=psa(i)/jz
ref(i)=ref(i)/nsum
birdie(i)=ref(i) !Copy ref into birdie
enddo
C Compute smoothed reference spectrum with narrow lines (birdies) removed
do i=4,nbins-3
rmax=-1.e10
do k=i-3,i+3 !Get highest point within +/- 3 bins
if(ref(k).gt.rmax) then
rmax=ref(k)
kpk=k
endif
enddo
sum=0.
nsum=0
do k=i-3,i+3
if(abs(k-kpk).gt.1) then
sum=sum+ref(k)
nsum=nsum+1
endif
enddo
ref2(i)=sum/nsum
enddo
call move(ref2(4),ref(4),nbins-6) !Copy smoothed ref back into ref
call pctile(ref(4),ref2,nbins-6,50,xmedian) !Get median in-band level
C Fix ends of reference spectrum
do i=1,3
ref(i)=ref(4)
ref(nbins+1-i)=ref(nbins-3)
enddo
facmax=30.0/xmedian
do i=1,nbins !Flatten the 2d spectrum
fac=xmedian/ref(i)
fac=min(fac,facmax)
do j=1,jz
s2(i,j)=fac*s2(i,j)
enddo
psa(i)=dB(psa(i)) + 25.
ref(i)=dB(ref(i)) + 25.
birdie(i)=db(birdie(i)) + 25.
enddo
900 continue
return
end