Now computing (but not yet plotting) spectra s1 and s2.

git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/jtms3@2508 ab8295b8-cf94-4d9e-aec4-7959e3be5d79

libm65/specjtms.f90

@@ -1,185 +1,88 @@
-subroutine specjtms(k)
+subroutine specjtms(k,px,pxsmo,spk0,f0)
 
 ! Starting code for a JTMS3 decoder.
 
   parameter (NSMAX=30*48000)
-  parameter (NFFT=16384,NH=NFFT/2)
+  parameter (MAXFFT=8192)
-  character*22 decoded
-  character*72 c72
   integer*2 id
-  real x(NFFT),w(NFFT)
+  real x(MAXFFT)
-  real p(24)
+  complex cx(MAXFFT),cx2(MAXFFT)
-  real rsent(258),softsym(683),sym2(258)
+  logical first
-  integer nsum(24)
+  common/mscom/id(1440000),s1(215),s2(215)
-  complex cx(NFFT),cx2(NFFT),cx0(NFFT)
-  complex covx(NH)
-  real s1a(NH),s2a(NH)
-  integer mettab(0:255,0:1)             !Metric table
-  integer data4a(9)                     !Decoded data (8-bit byte values)
-  integer data4(12)                     !Decoded data (6-bit byte values)
-  integer*1 data1(13)
-  integer*1 isync(43)
-  integer*1 chansym1(258),datsym2(215)
-  logical first,window
-  integer*1 i1
-  equivalence (i1,i4)
-  common/mscom/id(1440000),s1(215,703),s2(215,703)
   data first/.true./
   save
  
   if(first) then
      pi=4.0*atan(1.0)
      twopi=2.0*pi
-     do i=1,nfft
+     kstep=2048
-        w(i)=(sin(i*pi/nfft))**2
-     enddo
-     df=48000.0/nfft
-     ja=nint(2600.0)/df
-     jb=nint(3400.0)/df
-     iz=3000.0/df
-     covx=0.
-     kstep=4096
-     read(10,3001) rsent
-3001 format(50f1.0)
-     do i=1,258,6
-        rsent(i)=0.
-     enddo
-     rsent=2.0*rsent - 1.0
-
-     open(11,file='bpskmetrics.dat',status='old')
-     bias=0.5
-     scale=20.0
-     do i=0,255
-        read(11,*) xjunk,x0,x1
-        mettab(i,0)=nint(scale*(x0-bias))
-        mettab(i,1)=nint(scale*(x1-bias))
-     enddo
-     close(11)
-     window=.false.
      first=.false.
+     sqsmo=0.
   endif
 
+  t=k/48000.0
+  nfft=4096
+  df=48000.0/nfft
+  nh=nfft/2
   ib=k
   ia=k-kstep+1
   i0=k-nfft+1
   sq=0.
   do i=ia,ib
-     sq=sq + (0.001*id(i))**2
+     d=id(i)
+     sq=sq + d*d
   enddo
-  write(13,1010) t,sq,db(sq)
+  sq=sq/2048.0
-1010 format(3f12.3)
+  sqsmo=0.95*sqsmo + 0.05*sq
+  rms=sqrt(sq)
+  px=db(sq) - 23.0
+  pxsmo=db(sqsmo) - 23.0
+!  write(13,1010) t,rms,sq,px,pxsmo
+!1010 format(5f12.3)
   if(k.lt.nfft) return
 
-  x(1:nfft)=0.001*id(i0:ib)
 
-  fac=2.0/nfft
+  x(1:nfft)=id(i0:ib)
+  fac=0.002/nfft
   cx=fac*x
-  if(window) cx=w*cx
   call four2a(cx,nfft,1,-1,1)                    !Forward c2c FFT
 
+  iz=nint(2500.0/df)
+
   do i=1,iz                                      !Save spectrum for plot
-     s1a(i)=real(cx(i+1))**2 + aimag(cx(i+1))**2
+     s1(i)=real(cx(i+1))**2 + aimag(cx(i+1))**2
   enddo
 
   cx(1)=0.5*cx(1)
   cx(nh+2:nfft)=0.
   call four2a(cx,nfft,1,1,1)                     !Inverse c2c FFT
-  if(window) then
-     cx(1:nh)=cx(1:nh)+covx(1:nh)               !Add previous segment's 2nd half
-     covx(1:nh)=cx(nh+1:nfft)                    !Save 2nd half
-  endif
 
-  t=k/48000.0
+  cx2(1:nfft)=cx(1:nfft)*cx(1:nfft)
-  cx2=cx*cx
+  cx2(nfft+1:)=0.0
+
+  nfft=8192
+  df=48000.0/nfft
+  
   call four2a(cx2,nfft,1,-1,1)                   !Forward c2c FFT of cx2
 
   spk0=0.
+  j0=nint(2.0*1428.57/df)
+  ja=j0-107
+  jb=j0+107
   do j=ja,jb
-     sq=1.e-6*(real(cx2(j))**2 + aimag(cx2(j))**2)
+     sq=1.e-4*(real(cx2(j))**2 + aimag(cx2(j))**2)
-     s2a(j)=sq
+     s2(j-ja+1)=sq
-     f=(j-1)*df
      if(sq.gt.spk0) then
         spk0=sq
+        f=(j-1)*df
         f0=0.5*(f-3000.0)
         phi0=0.5*atan2(aimag(cx2(j)),real(cx2(j)))
      endif
-     write(15,1020) f,sq
-1020 format(f10.3,f12.3)
   enddo
+  spk0=0.5*db(spk0) - 7.0
+!  write(14,3001) k/2048,spk0,f0,phi0
+!3001 format(i8,3f12.3)
+  call flush(14)
 
-  slimit=2.5
+  return
-!  if(spk0.gt.slimit) then
-  if(abs(spk0-43.5).lt.0.1) then
-     write(*,1030) t,f0,phi0,spk0
-1030 format('t:',f6.2,'   f0:',f7.1,'   phi0:',f6.2,'   spk0:',f8.1)
-
-     phi=phi0
-     phi=3.9                                   !### test ###
-     dphi=twopi*(f0+1500.0 -1.1)/48000.0
-     p=0.
-     nsum=0
-     do i=1,nfft
-        phi=phi+dphi
-        if(phi.gt.twopi) phi=phi-twopi        
-        cx0(i)=cx(i)*cmplx(cos(phi),-sin(phi))
-        pha=atan2(aimag(cx0(i)),real(cx0(i)))
-!        write(18,1060) i,cx0(i),pha
-!1060    format(i6,5f12.3)
-        j=mod(i-1,24) + 1
-!        p(j)=p(j)+abs(cx0(i))
-        p(j)=p(j) + real(cx0(i))**2 + aimag(cx0(i))**2
-        nsum(j)=nsum(j)+1
-     enddo
-
-     do i=1,24
-        p(i)=p(i)/nsum(i)
-        write(20,1070) i,p(i)
-1070    format(i6,f12.3)
-     enddo
-
-     do i=19,nfft,24
-        amp=abs(cx0(i))
-        pha=atan2(aimag(cx0(i)),real(cx0(i)))
-        j=(i+23)/24
-        write(21,1060) j,cx0(i),pha,pha+twopi,amp
-1060    format(i6,5f12.3)
-        softsym(j)=real(cx0(i))
-     enddo
-
-!     do iter=1,5
-     rsent=cshift(rsent,86)
-     lagmax=nfft/24 - 258
-     do lag=0,lagmax
-        sum=dot_product(rsent,softsym(lag+1:lag+258))
-        if(abs(sum).gt.smax) then
-           smax=abs(sum)
-           lagpk=lag
-        endif
-        write(22,1080) lag,sum
-1080    format(i3,f12.3)
-     enddo
-!     rsent=cshift(rsent,43)
-!     enddo
-
-     do i=1,258
-        j=mod(i-1+2580,258) + 1
-        prod=rsent(j)*softsym(lagpk+i)
-        nchsym=nint(0.5*(rsent(j)+1.0))
-        write(23,1090) i,prod,rsent(j),softsym(lagpk+i),j,nchsym,lagpk+i
-1090    format(i5,3f10.3,3i5)
-        
-     enddo
-     
-     sym2=softsym(lagpk+1:lagpk+258)
-     sym2=cshift(sym2,-86)
-     do i=1,258
-        i4=128 + nint(6.0*sym2(i))
-        if(i4.lt.0) i4=0
-        if(i4.gt.255) i4=255
-        chansym1(i)=i1
-        write(24,2001) i,sym2(i),i4,chansym1(i)
-2001    format(i6,f8.3,2i6)
-     enddo
-  endif
-
 end subroutine specjtms

mainwindow.cpp

@@ -392,6 +392,7 @@ void MainWindow::dataSink(int k)
   static int nwrite=0;
   static int k0=99999999;
   static float px=0.0;
+  static float pxsmo,spk0,f0;
   static float sqave=0.0;
   static float green[704];
   static int ig=0;
@@ -418,7 +419,8 @@ void MainWindow::dataSink(int k)
   sqave=0.95*sqave + 0.05*sq;
   float pxave=10.0*log10(sqave/2048.0) - 23.0;
 
-//  specjtms_(&k,&px);
+  specjtms_(&k,&px,&pxsmo,&spk0,&f0);
+  if(spk0 > 3.0) qDebug() << (k-2048.0)/48000.0 << spk0 << f0;
   QString t;
   t.sprintf(" Rx noise: %5.1f ",pxave);
   lab2->setText(t);

mainwindow.h

@@ -235,10 +235,8 @@ extern void getDev(int* numDevices,char hostAPI_DeviceName[][50],
 
 extern "C" {
 //----------------------------------------------------- C and Fortran routines
-  void specjtms_(int* k, float* px);
+  void specjtms_(int* k, float* px, float* pxsmo, float* spk0, float* f0);
-
   void genjtms3_(char* message, short iwave[], int* nwave, int len1);
-
   void makepings_(short iwave[], int* nwave);
 
   void gen65_(char* msg, int* mode65, double* samfac, int* nsendingsh,

plotter.cpp

@@ -134,7 +134,7 @@ void CPlotter::draw(float green[], int ig)                       //draw()
   for(i=0; i<ig; i++) {
     y = green[i];
     painter1.drawPoint(i,0);
-    int y2 = 7*(y-m_plotZero);
+    int y2 = 4*(y-m_plotZero);
     if (y2<0) y2=0;
     if (y2>254) y2=254;
     LineBuf[j].setX(i);