From 945b6ef048264d0f59fc60ebe0b44d77971f78b6 Mon Sep 17 00:00:00 2001 From: Joe Taylor Date: Fri, 13 Jul 2012 00:45:43 +0000 Subject: [PATCH] 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 | 185 +++++++++++--------------------------------- mainwindow.cpp | 4 +- mainwindow.h | 4 +- plotter.cpp | 2 +- 4 files changed, 49 insertions(+), 146 deletions(-) diff --git a/libm65/specjtms.f90 b/libm65/specjtms.f90 index d823afe5e..9d2c2ed74 100644 --- a/libm65/specjtms.f90 +++ b/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) - character*22 decoded - character*72 c72 + parameter (MAXFFT=8192) integer*2 id - real x(NFFT),w(NFFT) - real p(24) - real rsent(258),softsym(683),sym2(258) - integer nsum(24) - 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) + real x(MAXFFT) + complex cx(MAXFFT),cx2(MAXFFT) + logical first + common/mscom/id(1440000),s1(215),s2(215) data first/.true./ save - + if(first) then pi=4.0*atan(1.0) twopi=2.0*pi - do i=1,nfft - 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. + kstep=2048 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) -1010 format(3f12.3) + sq=sq/2048.0 + 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=cx*cx + cx2(1:nfft)=cx(1:nfft)*cx(1:nfft) + 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) - s2a(j)=sq - f=(j-1)*df + sq=1.e-4*(real(cx2(j))**2 + aimag(cx2(j))**2) + s2(j-ja+1)=sq 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 -! 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 - + return end subroutine specjtms diff --git a/mainwindow.cpp b/mainwindow.cpp index f9c029033..0825ab8ba 100644 --- a/mainwindow.cpp +++ b/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); diff --git a/mainwindow.h b/mainwindow.h index b4cf634e7..f2707662b 100644 --- a/mainwindow.h +++ b/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, diff --git a/plotter.cpp b/plotter.cpp index b79db09b6..46ceffcfc 100644 --- a/plotter.cpp +++ b/plotter.cpp @@ -134,7 +134,7 @@ void CPlotter::draw(float green[], int ig) //draw() for(i=0; i254) y2=254; LineBuf[j].setX(i);