From 0e12c8f3f4276bfa9c7dee6b8c485dcd4dbab62f Mon Sep 17 00:00:00 2001 From: Joe Taylor Date: Mon, 12 Dec 2022 13:01:32 -0500 Subject: [PATCH] Get rid of more xpol stuff. --- q65w/libm65/symspec.f90 | 52 ++++--------------------- q65w/libm65/timf2.f90 | 84 +++-------------------------------------- q65w/mainwindow.cpp | 5 +-- q65w/mainwindow.h | 4 +- 4 files changed, 16 insertions(+), 129 deletions(-) diff --git a/q65w/libm65/symspec.f90 b/q65w/libm65/symspec.f90 index 1a0dd90a6..3fcf29a1c 100644 --- a/q65w/libm65/symspec.f90 +++ b/q65w/libm65/symspec.f90 @@ -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 diff --git a/q65w/libm65/timf2.f90 b/q65w/libm65/timf2.f90 index 4873d0e95..681e1dc24 100644 --- a/q65w/libm65/timf2.f90 +++ b/q65w/libm65/timf2.f90 @@ -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 diff --git a/q65w/mainwindow.cpp b/q65w/mainwindow.cpp index edb645e55..73a8bc75b 100644 --- a/q65w/mainwindow.cpp +++ b/q65w/mainwindow.cpp @@ -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); diff --git a/q65w/mainwindow.h b/q65w/mainwindow.h index 3cd57bf9c..442a28ae2 100644 --- a/q65w/mainwindow.h +++ b/q65w/mainwindow.h @@ -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[]);