2007-01-11 16:25:52 -05:00
|
|
|
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
|
