From 887f0c99b488901911720413b4b7d80cc15f30ea Mon Sep 17 00:00:00 2001 From: Joe Taylor Date: Tue, 4 Sep 2012 17:43:23 +0000 Subject: [PATCH] OK, I think we're now computing array ss(4,322,32768) correctly for the B2 and C2 modes. git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/map65@2555 ab8295b8-cf94-4d9e-aec4-7959e3be5d79 --- libm65/symspec.f90 | 145 +++++++++++++++++++++++++++++++-------------- mainwindow.cpp | 18 +++--- mainwindow.h | 12 ++-- 3 files changed, 114 insertions(+), 61 deletions(-) diff --git a/libm65/symspec.f90 b/libm65/symspec.f90 index 647977775..8bd26be87 100644 --- a/libm65/symspec.f90 +++ b/libm65/symspec.f90 @@ -1,5 +1,5 @@ -subroutine symspec(k,nxpol,ndiskdat,nb,nbslider,idphi,nfsample,fgreen, & - iqadjust,iqapply,gainx,gainy,phasex,phasey,rejectx,rejecty, & +subroutine symspec(k,nfast, nxpol,ndiskdat,nb,nbslider,idphi,nfsample, & + fgreen,iqadjust,iqapply,gainx,gainy,phasex,phasey,rejectx,rejecty, & pxdb,pydb,ssz5a,nkhz,ihsym,nzap,slimit,lstrong) ! k pointer to the most recent new data @@ -23,10 +23,11 @@ subroutine symspec(k,nxpol,ndiskdat,nb,nbslider,idphi,nfsample,fgreen, & real*8 ts,hsym real*8 fcenter common/datcom/dd(4,5760000),ss(4,322,NFFT),savg(4,NFFT),fcenter,nutc,junk(36) - real*4 ssz5a(NFFT),w(NFFT) + real*4 ssz5a(NFFT),w(NFFT),w2a(NFFT),w2b(NFFT) complex z,zfac complex zsumx,zsumy complex cx(NFFT),cy(NFFT) + complex cx00(NFFT),cy00(NFFT) complex cx0(0:1023),cx1(0:1023) complex cy0(0:1023),cy1(0:1023) logical*1 lstrong(0:1023) @@ -40,10 +41,19 @@ subroutine symspec(k,nxpol,ndiskdat,nb,nbslider,idphi,nfsample,fgreen, & endif if(k0.eq.99999999) then pi=4.0*atan(1.0) + w2a=0. + w2b=0. do i=1,NFFT - w(i)=(sin(i*pi/NFFT))**2 + w(i)=(sin(i*pi/NFFT))**2 !Window for nfast=1 + if(i.lt.17833) w2a(i)=(sin(i*pi/17832.925))**2 !Window a for nfast=2 + j=i-8916 + if(j.lt.17833) w2b(i)=(sin(j*pi/17832.925))**2 !Window b for nfast=2 enddo endif + + hsym=2048.d0*96000.d0/11025.d0 !Samples per JT65 half-symbol + if(nfsample.eq.95238) hsym=2048.d0*95238.1d0/11025.d0 + if(k.lt.k0) then ts=1.d0 - hsym savg=0. @@ -90,15 +100,12 @@ subroutine symspec(k,nxpol,ndiskdat,nb,nbslider,idphi,nfsample,fgreen, & k1=k1+kstep enddo - hsym=2048.d0*96000.d0/11025.d0 !Samples per JT65 half-symbol - if(nfsample.eq.95238) hsym=2048.d0*95238.1d0/11025.d0 npts=NFFT !Samples used in each half-symbol FFT - ihsym=ihsym+1 - ja=ts+hsym !Index of first sample + ts=ts+hsym + ja=ts !Index of first sample jb=ja+npts-1 !Last sample - ts=ts+hsym i=0 fac=0.0002 dphi=idphi/57.2957795 @@ -133,54 +140,102 @@ subroutine symspec(k,nxpol,ndiskdat,nb,nbslider,idphi,nfsample,fgreen, & if(pxdb.gt.60.0) pxdb=60.0 if(pydb.gt.60.0) pydb=60.0 - cx=w*cx !Apply window for 2nd forward FFT - if(nxpol.ne.0) cy=w*cy + cx00=cx + if(nxpol.ne.0) cy00=cy - call four2a(cx,NFFT,1,1,1) !Second forward FFT - 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) + do mm=1,nfast + 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 - if(nxpol.ne.0) then - call four2a(cy,NFFT,1,1,1) - 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 + 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) - n=ihsym - do i=1,NFFT - sx=real(cx(i))**2 + aimag(cx(i))**2 - ss(1,n,i)=sx ! Pol = 0 - savg(1,i)=savg(1,i) + sx - if(nxpol.ne.0) then - z=cx(i) + cy(i) - s45=0.5*(real(z)**2 + aimag(z)**2) - ss(2,n,i)=s45 ! Pol = 45 - savg(2,i)=savg(2,i) + s45 + 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 - sy=real(cy(i))**2 + aimag(cy(i))**2 - ss(3,n,i)=sy ! Pol = 90 - savg(3,i)=savg(3,i) + sy + ihsym=ihsym+1 + n=ihsym + if(mm.eq.1) then + do i=1,NFFT + sx=real(cx(i))**2 + aimag(cx(i))**2 + ss(1,n,i)=sx ! Pol = 0 + savg(1,i)=savg(1,i) + sx + + if(nxpol.ne.0) then + z=cx(i) + cy(i) + s45=0.5*(real(z)**2 + aimag(z)**2) + ss(2,n,i)=s45 ! Pol = 45 + savg(2,i)=savg(2,i) + s45 + + sy=real(cy(i))**2 + aimag(cy(i))**2 + ss(3,n,i)=sy ! Pol = 90 + savg(3,i)=savg(3,i) + sy - z=cx(i) - cy(i) - s135=0.5*(real(z)**2 + aimag(z)**2) - ss(4,n,i)=s135 ! Pol = 135 - savg(4,i)=savg(4,i) + s135 + z=cx(i) - cy(i) + s135=0.5*(real(z)**2 + aimag(z)**2) + ss(4,n,i)=s135 ! Pol = 135 + savg(4,i)=savg(4,i) + s135 - z=cx(i)*conjg(cy(i)) - q=sx - sy - u=2.0*real(z) - ssz5a(i)=0.707*sqrt(q*q + u*u) !Spectrum of linear polarization + z=cx(i)*conjg(cy(i)) + q=sx - sy + u=2.0*real(z) + ssz5a(i)=0.707*sqrt(q*q + u*u) !Spectrum of linear polarization ! Leif's formula: ! ssz5a(i)=0.5*(sx+sy) + (real(z)**2 + aimag(z)**2 - sx*sy)/(sx+sy) + else + ssz5a(i)=sx + endif + enddo else - ssz5a(i)=sx + do i=1,NFFT + sx=real(cx(i))**2 + aimag(cx(i))**2 + ss(1,n,i)=ss(1,n,i) + sx ! Pol = 0 + savg(1,i)=savg(1,i) + sx + + if(nxpol.ne.0) then + z=cx(i) + cy(i) + s45=0.5*(real(z)**2 + aimag(z)**2) + ss(2,n,i)=ss(2,n,i) + s45 ! Pol = 45 + savg(2,i)=savg(2,i) + s45 + + sy=real(cy(i))**2 + aimag(cy(i))**2 + ss(3,n,i)=ss(3,n,i) + sy ! Pol = 90 + savg(3,i)=savg(3,i) + sy + + z=cx(i) - cy(i) + s135=0.5*(real(z)**2 + aimag(z)**2) + ss(4,n,i)=ss(4,n,i) + s135 ! Pol = 135 + savg(4,i)=savg(4,i) + s135 + + z=cx(i)*conjg(cy(i)) + q=sx - sy + u=2.0*real(z) + ssz5a(i)=ssz5a(i) + 0.707*sqrt(q*q + u*u) !Spectrum of lin pol + else + ssz5a(i)=ssz5a(i) + sx + endif + enddo endif enddo + if(ihsym.eq.278) then if(iqadjust.ne.0 .and. ipkx.ne.0 .and. ipky.ne.0) then rejectx=10.0*log10(savg(1,1+nfft-ipkx)/savg(1,1+ipkx)) diff --git a/mainwindow.cpp b/mainwindow.cpp index 61a08837a..acc8b985a 100644 --- a/mainwindow.cpp +++ b/mainwindow.cpp @@ -497,10 +497,10 @@ void MainWindow::dataSink(int k) fgreen=(float)g_pWideGraph->fGreen(); nadj++; if(m_adjustIQ==0) nadj=0; - symspec_(&k, &nxpol, &ndiskdat, &nb, &m_NBslider, &m_dPhi, &nfsample, - &fgreen, &m_adjustIQ, &m_applyIQcal, &m_gainx, &m_gainy, &m_phasex, - &m_phasey, &rejectx, &rejecty, &px, &py, s, &nkhz, - &ihsym, &nzap, &slimit, lstrong); + symspec_(&k, &m_nfast, &nxpol, &ndiskdat, &nb, &m_NBslider, &m_dPhi, + &nfsample, &fgreen, &m_adjustIQ, &m_applyIQcal, + &m_gainx, &m_gainy, &m_phasex, &m_phasey, &rejectx, &rejecty, + &px, &py, s, &nkhz, &ihsym, &nzap, &slimit, lstrong); QString t; m_pctZap=nzap/178.3; if(m_xpol) t.sprintf(" Rx noise: %5.1f %5.1f %5.1f %% ",px,py,m_pctZap); @@ -549,7 +549,7 @@ void MainWindow::dataSink(int k) ntrz=ntr; n=0; } - if(ihsym == 279/m_nfast) { + if(ihsym == 280) { datcom_.newdat=1; datcom_.nagain=0; QDateTime t = QDateTime::currentDateTimeUtc(); @@ -1038,7 +1038,7 @@ void MainWindow::diskDat() //diskDat() if(m_fs96000) hsym=2048.0*96000.0/11025.0; //Samples per JT65 half-symbol if(!m_fs96000) hsym=2048.0*95238.1/11025.0; - for(int i=0; i<282/m_nfast; i++) { // Do the half-symbol FFTs + for(int i=0; i<284/m_nfast; i++) { // Do the half-symbol FFTs int k = i*hsym + 2048.5; dataSink(k); if(i%10 == 0) qApp->processEvents(); //Keep the GUI responsive @@ -1399,10 +1399,9 @@ void MainWindow::guiUpdate() ba2msg(ba,message); int len1=22; int mode65=m_mode65; - int nfast=m_nfast; double samfac=1.0; - gen65_(message,&mode65,&nfast,&samfac,&nsendingsh,msgsent,iwave, + gen65_(message,&mode65,&m_nfast,&samfac,&nsendingsh,msgsent,iwave, &nwave,len1,len1); msgsent[22]=0; @@ -1821,8 +1820,7 @@ void MainWindow::msgtype(QString t, QLineEdit* tx) //msgtype() int i1=t.indexOf(" OOO"); QByteArray s=t.toUpper().toLocal8Bit(); ba2msg(s,message); - int nfast=m_nfast; - gen65_(message,&mode65,&nfast,&samfac,&nsendingsh,msgsent,iwave, + gen65_(message,&mode65,&m_nfast,&samfac,&nsendingsh,msgsent,iwave, &mwave,len1,len1); QPalette p(tx->palette()); diff --git a/mainwindow.h b/mainwindow.h index a68effa3f..968bed887 100644 --- a/mainwindow.h +++ b/mainwindow.h @@ -267,12 +267,12 @@ extern void getDev(int* numDevices,char hostAPI_DeviceName[][50], extern "C" { //----------------------------------------------------- C and Fortran routines - void symspec_(int* k, int* nxpol, int* ndiskdat, int* nb, int* m_NBslider, - int* idphi, int* nfsample, float* fgreen, int* iqadjust, - int* iqapply, 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[]); + void symspec_(int* k, int* nfast, int* nxpol, int* ndiskdat, int* nb, + int* m_NBslider, int* idphi, int* nfsample, float* fgreen, + int* iqadjust, int* iqapply, 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[]); void gen65_(char* msg, int* mode65, int* nfast, double* samfac, int* nsendingsh, char* msgsent, short iwave[], int* nwave,