WSJT-X/symspec.f90
Joe Taylor f419b84d02 Option to accept data from Linrad in floating-point format.
git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/map65@2447 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
2011-09-28 18:25:26 +00:00

174 lines
4.6 KiB
Fortran

subroutine symspec(dd,kbuf,kk,kkdone,nutc,newdat)
! Compute spectra at four polarizations, using half-symbol steps.
parameter (NSMAX=60*96000)
real*4 dd(4,NSMAX,2)
complex z
real*8 ts,hsym
include 'spcom.f90'
include 'gcom2.f90'
complex cx(NFFT),cy(NFFT) ! pad to 32k with zeros
data kbuf0/-999/,n/0/
save
kkk=kk
if(kbuf.eq.2) kkk=kk-5760000
fac=0.0002
hsym=2048.d0*96000.d0/11025.d0 !Samples per half symbol
npts=hsym !Integral samples per half symbol
ntot=322 !Half symbols per transmission
! ntot=279 !Half symbols in 51.8 sec
if(kbuf.ne.kbuf0 .or. ndiskdat.eq.1) then
kkdone=0
kbuf0=kbuf
ts=1.d0 - hsym
n=0
do ip=1,4
do i=1,NFFT
szavg(ip,i)=0.
enddo
enddo
! Get baseline power level for this minute
n1=200 !Block size (somewhat arbitrary)
n2=(kkk-kkdone)/n1 !Number of blocks
k=0 !Starting place
sqq=0.
nsqq=0
do j=1,n2
sq=0.
do i=1,n1 !Find power in each block
k=k+1
x1=dd(1,k,kbuf)
x2=dd(2,k,kbuf)
x3=dd(3,k,kbuf)
x4=dd(4,k,kbuf)
sq=sq + x1*x1 + x2*x2 + x3*x3 + x4*x4
enddo
if(sq.lt.n1*10000.) then !Find power in good blocks
sqq=sqq+sq
nsqq=nsqq+1
endif
enddo
sqave=sqq/nsqq !Average power in good blocks
nclip=0
nz2=0
endif
if(nblank.ne.0) then
! Apply final noise blanking
n2=(kkk-kkdone)/n1
k=kkdone
do j=1,n2
sq=0.
do i=1,n1
k=k+1
x1=dd(1,k,kbuf)
x2=dd(2,k,kbuf)
x3=dd(3,k,kbuf)
x4=dd(4,k,kbuf)
sq=sq + x1*x1 + x2*x2 + x3*x3 + x4*x4
enddo
! If power in this block is excessive, blank it.
if(sq.gt.1.5*sqave) then
do i=k-n1+1,k
dd(1,i,kbuf)=0
dd(2,i,kbuf)=0
dd(3,i,kbuf)=0
dd(4,i,kbuf)=0
enddo
nclip=nclip+1
endif
enddo
nz2=nz2+n2
pctblank=nclip*100.0/nz2
else
pctblank=0.
endif
do nn=1,ntot
i0=ts+hsym !Starting sample pointer
if((i0+npts-1).gt.kkk) go to 998 !See if we have enough points
i1=ts+2*hsym !Next starting sample pointer
ts=ts+hsym !OK, update the exact sample pointer
do i=1,npts !Copy data to FFT arrays
xr=fac*dd(1,i0+i,kbuf)
xi=fac*dd(2,i0+i,kbuf)
cx(i)=cmplx(xr,xi)
yr=fac*dd(3,i0+i,kbuf)
yi=fac*dd(4,i0+i,kbuf)
cy(i)=cmplx(yr,yi)
enddo
do i=npts+1,NFFT !Pad to 32k with zeros
cx(i)=0.
cy(i)=0.
enddo
call four2a(cx,NFFT,1,1,1) !Do the FFTs
call four2a(cy,NFFT,1,1,1)
n=n+1
do i=1,NFFT !Save and accumulate power spectra
sx=real(cx(i))**2 + aimag(cx(i))**2
ssz(1,n,i)=sx ! Pol = 0
szavg(1,i)=szavg(1,i) + sx
z=cx(i) + cy(i)
s45=0.5*(real(z)**2 + aimag(z)**2)
ssz(2,n,i)=s45 ! Pol = 45
szavg(2,i)=szavg(2,i) + s45
sy=real(cy(i))**2 + aimag(cy(i))**2
ssz(3,n,i)=sy ! Pol = 90
szavg(3,i)=szavg(3,i) + sy
z=cx(i) - cy(i)
s135=0.5*(real(z)**2 + aimag(z)**2)
ssz(4,n,i)=s135 ! Pol = 135
szavg(4,i)=szavg(4,i) + s135
z=cx(i)*conjg(cy(i))
! Leif's formula:
! ss5(n,i)=0.5*(sx+sy) + (real(z)**2 + aimag(z)**2 -
! + sx*sy)/(sx+sy)
! Leif's suggestion:
! ss5(n,i)=max(sx,s45,sy,s135)
! Linearly polarized component, from the Stokes parameters:
q=sx - sy
u=2.0*real(z)
! v=2.0*aimag(z)
ssz5(n,i)=0.707*sqrt(q*q + u*u)
enddo
! if(n.eq.ntot) then
if(n.ge.279) then
call move(ssz5,ss5,322*NFFT)
call cs_lock('symspec')
write(utcdata,1002) nutc
1002 format(i4.4)
call cs_unlock
utcdata=utcdata(1:2)//':'//utcdata(3:4)
newspec=1
call move(ssz,ss,4*322*NFFT)
call move(szavg,savg,4*NFFT)
newdat=1
ndecoding=1
go to 999
endif
kkdone=i1-1
nhsym=n
call sleep_msec(0)
enddo
998 kkdone=i1-1
999 continue
return
end subroutine symspec