Get rid of more xpol stuff.

This commit is contained in:
Joe Taylor 2022-12-12 13:01:32 -05:00
parent fc2273dc67
commit 0e12c8f3f4
4 changed files with 16 additions and 129 deletions

View File

@ -1,16 +1,13 @@
subroutine symspec(k,nxpol,ndiskdat,nb,nbslider,idphi,nfsample, &
fgreen,iqadjust,iqapply,gainx,gainy,phasex,phasey,rejectx,rejecty, &
subroutine symspec(k,ndiskdat,nb,nbslider,idphi,nfsample, &
fgreen,gainx,gainy,phasex,phasey,rejectx,rejecty, &
pxdb,pydb,ssz5a,nkhz,ihsym,nzap,slimit,lstrong)
! k pointer to the most recent new data
! nxpol 0/1 to indicate single- or dual-polarization
! ndiskdat 0/1 to indicate if data from disk
! nb 0/1 status of noise blanker
! idphi Phase correction for Y channel, degrees
! nfsample sample rate (Hz)
! fgreen Frequency of green marker in I/Q calibrate mode (-48.0 to +48.0 kHz)
! iqadjust 0/1 to indicate whether IQ adjustment is active
! iqapply 0/1 to indicate whether to apply I/Q calibration
! pxdb power in x channel (0-60 dB)
! pydb power in y channel (0-60 dB)
! ssz5a polarized spectrum, for waterfall display
@ -28,10 +25,9 @@ subroutine symspec(k,nxpol,ndiskdat,nb,nbslider,idphi,nfsample, &
real*4 ssz5a(NFFT),w(NFFT),w2a(NFFT),w2b(NFFT)
complex z
complex zsumx,zsumy
complex cx(NFFT),cy(NFFT)
complex cx00(NFFT),cy00(NFFT)
complex cx(NFFT)
complex cx00(NFFT)
complex cx0(0:1023),cx1(0:1023)
complex cy0(0:1023),cy1(0:1023)
logical*1 lstrong(0:1023)
data rms/999.0/,k0/99999999/,nadjx/0/,nadjy/0/
save
@ -75,11 +71,7 @@ subroutine symspec(k,nxpol,ndiskdat,nb,nbslider,idphi,nfsample, &
px=0.
py=0.
iqapply0=0
iqadjust0=0
if(iqadjust.ne.0) iqapply0=0
nwindow=2
! nwindow=0 !### No windowing ###
nfft2=1024
kstep=nfft2
if(nwindow.ne.0) kstep=nfft2/2
@ -88,19 +80,14 @@ subroutine symspec(k,nxpol,ndiskdat,nb,nbslider,idphi,nfsample, &
j=k1+1
do i=0,nfft2-1
cx0(i)=cmplx(dd(1,j+i),dd(2,j+i))
if(nxpol.ne.0) cy0(i)=cmplx(dd(3,j+i),dd(4,j+i))
enddo
call timf2(k,nxpol,nfft2,nwindow,nb,peaklimit,iqadjust0,iqapply0, &
faclim,cx0,cy0,gainx,gainy,phasex,phasey,cx1,cy1,slimit,lstrong, &
call timf2(k,nfft2,nwindow,nb,peaklimit, &
faclim,cx0,gainx,gainy,phasex,phasey,cx1,slimit,lstrong, &
px,py,nzap)
do i=0,kstep-1
dd(1,j+i)=real(cx1(i))
dd(2,j+i)=aimag(cx1(i))
if(nxpol.ne.0) then
dd(3,j+i)=real(cy1(i))
dd(4,j+i)=aimag(cy1(i))
endif
enddo
k1=k1+kstep
enddo
@ -113,19 +100,11 @@ subroutine symspec(k,nxpol,ndiskdat,nb,nbslider,idphi,nfsample, &
i=0
fac=0.0002
do j=ja,jb !Copy data into cx, cy
do j=ja,jb !Copy data into cx
x1=dd(1,j)
x2=dd(2,j)
if(nxpol.ne.0) then
x3=dd(3,j)
x4=dd(4,j)
else
x3=0.
x4=0.
endif
i=i+1
cx(i)=fac*cmplx(x1,x2)
cy(i)=cmplx(x3,x4) !NB: cy includes dphi correction
enddo
if(nzap/178.lt.50 .and. (ndiskdat.eq.0 .or. ihsym.lt.280)) then
@ -134,7 +113,6 @@ subroutine symspec(k,nxpol,ndiskdat,nb,nbslider,idphi,nfsample, &
rmsx=sqrt(0.5*px/nsum)
rmsy=sqrt(0.5*py/nsum)
rms=rmsx
if(nxpol.ne.0) rms=sqrt((px+py)/(4.0*nsum))
endif
pxdb=0.
pydb=0.
@ -144,36 +122,20 @@ subroutine symspec(k,nxpol,ndiskdat,nb,nbslider,idphi,nfsample, &
if(pydb.gt.60.0) pydb=60.0
cx00=cx
if(nxpol.ne.0) cy00=cy
do mm=1,nfast
ihsym=ihsym+1
if(nfast.eq.1) then
cx=w*cx00 !Apply window for 2nd forward FFT
if(nxpol.ne.0) cy=w*cy00
else
if(mm.eq.1) then
cx=w2a*cx00
if(nxpol.ne.0) cy=w2a*cy00
else
cx=w2b*cx00
if(nxpol.ne.0) cy=w2b*cy00
endif
endif
call four2a(cx,NFFT,1,1,1) !Second forward FFT (X)
if(iqadjust.eq.0) nadjx=0
if(iqadjust.ne.0 .and. nadjx.lt.50) call iqcal(nadjx,cx,NFFT, &
gainx,phasex,zsumx,ipkx,rejectx0)
if(iqapply.ne.0) call iqfix(cx,NFFT,gainx,phasex)
if(nxpol.ne.0) then
call four2a(cy,NFFT,1,1,1) !Second forward FFT (Y)
if(iqadjust.eq.0) nadjy=0
if(iqadjust.ne.0 .and. nadjy.lt.50) call iqcal(nadjy,cy,NFFT, &
gainy,phasey,zsumy,ipky,rejecty)
if(iqapply.ne.0) call iqfix(cy,NFFT,gainy,phasey)
endif
n=min(322,ihsym)
do i=1,NFFT

View File

@ -1,16 +1,15 @@
subroutine timf2(k,nxpol,nfft,nwindow,nb,peaklimit,iqadjust,iqapply,faclim, &
cx0,cy0,gainx,gainy,phasex,phasey,cx1,cy1,slimit,lstrong,px,py,nzap)
subroutine timf2(k,nfft,nwindow,nb,peaklimit,faclim, &
cx0,gainx,gainy,phasex,phasey,cx1,slimit,lstrong,px,py,nzap)
! Sequential processing of time-domain I/Q data, using Linrad-like
! "first FFT" and "first backward FFT".
! cx0,cy0 - complex input data
! cx0 - complex input data
! nfft - length of FFTs
! nwindow - 0 for no window, 2 for sin^2 window
! iqapply - 0/1 determines if I/Q phase and amplitude corrections applied
! gainx,y - gain error in Q channel, relative to I
! phasex,y - phase error
! cx1,cy1 - output data
! cx1 - output data
! Non-windowed processing means no overlap, so kstep=nfft.
! Sin^2 window has 50% overlap, kstep=nfft/2.
@ -23,13 +22,9 @@ subroutine timf2(k,nxpol,nfft,nwindow,nb,peaklimit,iqadjust,iqapply,faclim, &
parameter (MAXFFT=1024,MAXNH=MAXFFT/2)
parameter (MAXSIGS=100)
complex cx0(0:nfft-1),cx1(0:nfft-1)
complex cy0(0:nfft-1),cy1(0:nfft-1)
complex cx(0:MAXFFT-1),cxt(0:MAXFFT-1)
complex cy(0:MAXFFT-1),cyt(0:MAXFFT-1)
complex cxs(0:MAXFFT-1),covxs(0:MAXNH-1) !Strong X signals
complex cys(0:MAXFFT-1),covys(0:MAXNH-1) !Strong Y signals
complex cxw(0:MAXFFT-1),covxw(0:MAXNH-1) !Weak X signals
complex cyw(0:MAXFFT-1),covyw(0:MAXNH-1) !Weak Y signals
real*4 w(0:MAXFFT-1)
real*4 s(0:MAXFFT-1)
logical*1 lstrong(0:MAXFFT-1),lprev
@ -40,7 +35,7 @@ subroutine timf2(k,nxpol,nfft,nwindow,nb,peaklimit,iqadjust,iqapply,faclim, &
data k0/99999999/
save
if(faclim+iqadjust.eq.-9999.0) iqadjust=0 !Silence compiler warning.
if(faclim.eq.-9999.0) stop !Silence compiler warning.
if(first) then
pi=4.0*atan(1.0)
do i=0,nfft-1
@ -66,53 +61,13 @@ subroutine timf2(k,nxpol,nfft,nwindow,nb,peaklimit,iqadjust,iqapply,faclim, &
cx(0:nfft-1)=cx0
if(nwindow.eq.2) cx(0:nfft-1)=w(0:nfft-1)*cx(0:nfft-1)
call four2a(cx,nfft,1,1,1) !First forward FFT (X)
if(nxpol.ne.0) then
cy(0:nfft-1)=cy0
if(nwindow.eq.2) cy(0:nfft-1)=w(0:nfft-1)*cy(0:nfft-1)
call four2a(cy,nfft,1,1,1) !First forward FFT (Y)
endif
if(iqapply.ne.0) then !Apply I/Q corrections (X)
h=gainx*cmplx(cos(phasex),sin(phasex))
v=0.
do i=0,nfft-1
u=cx(i)
if(i.gt.0) v=cx(nfft-i)
x=real(u) + real(v) - (aimag(u) + aimag(v))*aimag(h) + &
(real(u) - real(v))*real(h)
y=aimag(u) - aimag(v) + (aimag(u) + aimag(v))*real(h) + &
(real(u) - real(v))*aimag(h)
cxt(i)=0.5*cmplx(x,y)
enddo
else
cxt(0:nfft-1)=cx(0:nfft-1)
endif
if(nxpol.ne.0) then
if(iqapply.ne.0) then !Apply I/Q corrections (Y)
h=gainy*cmplx(cos(phasey),sin(phasey))
v=0.
do i=0,nfft-1
u=cy(i)
if(i.gt.0) v=cy(nfft-i)
x=real(u) + real(v) - (aimag(u) + aimag(v))*aimag(h) + &
(real(u) - real(v))*real(h)
y=aimag(u) - aimag(v) + (aimag(u) + aimag(v))*real(h) + &
(real(u) - real(v))*aimag(h)
cyt(i)=0.5*cmplx(x,y)
enddo
else
cyt(0:nfft-1)=cy(0:nfft-1)
endif
endif
cxt(0:nfft-1)=cx(0:nfft-1)
! Identify frequencies with strong signals, copy frequency-domain
! data into array cs (strong) or cw (weak).
do i=0,nfft-1
p=real(cxt(i))**2 + aimag(cxt(i))**2
if(nxpol.ne.0) p=p + real(cyt(i))**2 + aimag(cyt(i))**2
s(i)=p
enddo
ave=sum(s(0:nfft-1))/nfft
@ -151,50 +106,28 @@ subroutine timf2(k,nxpol,nfft,nwindow,nb,peaklimit,iqadjust,iqapply,faclim, &
if(lstrong(i)) then
cxs(i)=fac*cxt(i)
cxw(i)=0.
if(nxpol.ne.0) then
cys(i)=fac*cyt(i)
cyw(i)=0.
endif
else
cxw(i)=fac*cxt(i)
cxs(i)=0.
if(nxpol.ne.0) then
cyw(i)=fac*cyt(i)
cys(i)=0.
endif
endif
enddo
call four2a(cxw,nfft,1,-1,1) !Transform weak and strong X
call four2a(cxs,nfft,1,-1,1) !back to time domain, separately
if(nxpol.ne.0) then
call four2a(cyw,nfft,1,-1,1) !Transform weak and strong Y
call four2a(cys,nfft,1,-1,1) !back to time domain, separately
endif
if(nwindow.eq.2) then
cxw(0:nh-1)=cxw(0:nh-1)+covxw(0:nh-1) !Add previous segment's 2nd half
covxw(0:nh-1)=cxw(nh:nfft-1) !Save 2nd half
cxs(0:nh-1)=cxs(0:nh-1)+covxs(0:nh-1) !Ditto for strong signals
covxs(0:nh-1)=cxs(nh:nfft-1)
if(nxpol.ne.0) then
cyw(0:nh-1)=cyw(0:nh-1)+covyw(0:nh-1) !Add previous segment's 2nd half
covyw(0:nh-1)=cyw(nh:nfft-1) !Save 2nd half
cys(0:nh-1)=cys(0:nh-1)+covys(0:nh-1) !Ditto for strong signals
covys(0:nh-1)=cys(nh:nfft-1)
endif
endif
! Apply noise blanking to weak data
if(nb.ne.0) then
do i=0,kstep-1
peak=abs(cxw(i))
if(nxpol.ne.0) peak=max(peak,abs(cyw(i)))
if(peak.gt.peaklimit) then
cxw(i)=0.
if(nxpol.ne.0) cyw(i)=0.
nzap=nzap+1
endif
enddo
@ -203,13 +136,8 @@ subroutine timf2(k,nxpol,nfft,nwindow,nb,peaklimit,iqadjust,iqapply,faclim, &
! Compute power levels from weak data only
do i=0,kstep-1
px=px + real(cxw(i))**2 + aimag(cxw(i))**2
if(nxpol.ne.0) py=py + real(cyw(i))**2 + aimag(cyw(i))**2
enddo
cx1(0:kstep-1)=cxw(0:kstep-1) + cxs(0:kstep-1) !Weak + strong (X)
if(nxpol.ne.0) then
cy1(0:kstep-1)=cyw(0:kstep-1) + cys(0:kstep-1) !Weak + strong (Y)
endif
return
end subroutine timf2

View File

@ -328,7 +328,6 @@ void MainWindow::dataSink(int k)
static float fgreen;
static int ndiskdat;
static int nb;
static int nxpol=0;
static float px=0.0,py=0.0;
static uchar lstrong[1024];
static float rejectx;
@ -349,9 +348,7 @@ void MainWindow::dataSink(int k)
nfsample=96000;
if(!m_fs96000) nfsample=95238;
fgreen=m_wide_graph_window->fGreen();
int zero=0;
symspec_(&k, &nxpol, &ndiskdat, &nb, &m_NBslider, &m_dPhi,
&nfsample, &fgreen, &zero, &zero,
symspec_(&k, &ndiskdat, &nb, &m_NBslider, &m_dPhi, &nfsample, &fgreen,
&m_gainx, &m_gainy, &m_phasex, &m_phasey, &rejectx, &rejecty,
&px, &py, s, &nkhz, &ihsym, &nzap, &slimit, lstrong);

View File

@ -205,9 +205,9 @@ extern int killbyname(const char* progName);
extern "C" {
//----------------------------------------------------- C and Fortran routines
void symspec_(int* k, int* nxpol, int* ndiskdat, int* nb,
void symspec_(int* k, int* ndiskdat, int* nb,
int* m_NBslider, int* idphi, int* nfsample, float* fgreen,
int* iqadjust, int* iqapply, float* gainx, float* gainy,
float* gainx, float* gainy,
float* phasex, float* phasey, float* rejectx, float* rejecty,
float* px, float* py, float s[], int* nkhz, int* nhsym,
int* nzap, float* slimit, uchar lstrong[]);