diff --git a/lib/msk144signalquality.f90 b/lib/msk144signalquality.f90 new file mode 100644 index 000000000..0a26ee1eb --- /dev/null +++ b/lib/msk144signalquality.f90 @@ -0,0 +1,227 @@ +subroutine msk144signalquality(cframe,snr,freq,t0,softbits,msg,dxcall,nbiterrors,eyeopening,trained,pcoeffs) + + character*22 msg,msgsent + character*12 dxcall + character*12 training_dxcall + character*12 trained_dxcall + character*6 mygrid + character*40 pcoeff_filename + character*8 date + character*10 time + character*5 zone + + complex cframe(864) + complex cross(864) + complex cross_avg(864) + complex canalytic(1024) + complex cmodel(1024) + + integer i4tone(144) + integer hardbits(144) + integer msgbits(144) + integer values(8) + + logical*1 bcontest + logical*1 first + logical*1 currently_training + logical*1 trained + logical*1 msg_has_dxcall + logical*1 is_training_frame + + real abscross_avg(864) + real softbits(144) + real softav(144) + real waveform(0:863) + real d(1024) + real phase(864) + real twopi,freq,phi,dphi0,dphi1,dphi + real*8 x(145),y(145),pp(145),sigmay(145),a(5),chisqr,pcoeffs(5) + + data first/.true./ + save cross_avg,abscross_avg,wt_avg,first,currently_training, & + navg,softav,tlast,training_dxcall,trained_dxcall + + if (first) then + navg=0 + cross=cmplx(0.0,0.0) + cross_avg=cmplx(0.0,0.0) + wt_avg=0.0 + tlast=0.0 + trained_dxcall(1:12)=' ' + training_dxcall(1:12)=' ' + trained=.false. + currently_training=.false. + pcoeffs(1:5)=0.0 + first=.false. + endif + + if( (trained .and. (dxcall .ne. trained_dxcall)) .or. & + (currently_training .and. (dxcall .ne. training_dxcall)) .or. & + (navg .gt. 10 )) then !reset and retrain + navg=0 + cross=cmplx(0.0,0.0) + cross_avg=cmplx(0.0,0.0) + wt_avg=0.0 + tlast=0.0 + trained_dxcall(1:12)=' ' + trained=.false. + currently_training=.false. + training_dxcall(1:12)=' ' + trained_dxcall(1:12)=' ' + pcoeffs(1:5)=0.0 +write(*,*) 'reset to untrained state ' + endif + + indx_dxcall=index(msg,trim(dxcall)) + msg_has_dxcall = indx_dxcall .ge. 4 + + if( msg_has_dxcall .and. (.not. trained) .and. .not. currently_training ) then + currently_training=.true. + training_dxcall=trim(dxcall) + trained_dxcall(1:12)=' ' + pcoeffs(1:5)=0.0 +write(*,*) 'start training on call ',training_dxcall + endif + + if( msg_has_dxcall .and. currently_training ) then + trained=.false. ! just to be sure + trained_dxcall(1:12)=' ' + training_dxcall=dxcall + pcoeffs(1:5)=0.0 + endif + +! use decoded message to figure out how many bit errors in the frame + do i=1, 144 + hardbits(i)=0 + if(softbits(i) .gt. 0 ) hardbits(i)=1 + enddo + +! generate tones from decoded message + mygrid="EN50" + ichk=0 + bcontest=.false. + call genmsk144(msg,mygrid,ichk,bcontest,msgsent,i4tone,itype) + +! reconstruct message bits from tones + msgbits(1)=0 + do i=1,143 + if( i4tone(i) .eq. 0 ) then + if( mod(i,2) .eq. 1 ) then + msgbits(i+1)=msgbits(i) + else + msgbits(i+1)=mod(msgbits(i)+1,2) + endif + else + if( mod(i,2) .eq. 1 ) then + msgbits(i+1)=mod(msgbits(i)+1,2) + else + msgbits(i+1)=msgbits(i) + endif + endif + enddo + + nbiterrors=0 + do i=1,144 + if( hardbits(i) .ne. msgbits(i) ) nbiterrors=nbiterrors+1 + enddo + + nplus=0 + nminus=0 + eyetop=1 + eyebot=-1 + do i=1,144 + if( msgbits(i) .eq. 1 ) then + if( softbits(i) .lt. eyetop ) eyetop=softbits(i) + else + if( softbits(i) .gt. eyebot ) eyebot=softbits(i) + endif + enddo + eyeopening=eyetop-eyebot + + is_training_frame = & + (snr.gt.5.0 .and.(nbiterrors.lt.7)) .and. & + (abs(t0-tlast) .gt. 0.072) .and. & + msg_has_dxcall + if( currently_training .and. is_training_frame ) then + twopi=8.0*atan(1.0) + nsym=144 + if( i4tone(41) .lt. 0 ) nsym=40 + dphi0=twopi*(freq-500)/12000.0 + dphi1=twopi*(freq+500)/12000.0 + phi=-twopi/8 + indx=0 + do i=1,nsym + if( i4tone(i) .eq. 0 ) then + dphi=dphi0 + else + dphi=dphi1 + endif + do j=1,6 + waveform(indx)=cos(phi); + indx=indx+1 + phi=mod(phi+dphi,twopi) + enddo + enddo +! convert the passband waveform to complex baseband + npts=864 + nfft=1024 + d=0 + d(1:864)=waveform(0:863) + a=0 + call analytic(d,npts,nfft,canalytic,a,.false.) ! don't equalize the model + call tweak1(canalytic,nfft,-freq,cmodel) + call four2a(cframe(1:864),864,1,-1,1) + call four2a(cmodel(1:864),864,1,-1,1) + +! Cross spectra from different messages can be averaged +! as long as all messages originate from dxcall. + cross=cmodel(1:864)*conjg(cframe)/1000.0 + cross=cshift(cross,864/2) + cross_avg=cross_avg+10**(snr/20.0)*cross + wt_avg=wt_avg+10**(snr/20.0) + navg=navg+1 + tlast=t0 + phase=atan2(imag(cross_avg),real(cross_avg)) + df=12000.0/864.0 + nm=145 + do i=1,145 + x(i)=(i-73)*df/1000.0 + enddo + y=phase((864/2-nm/2):(864/2+nm/2)) + sigmay=wt_avg/abs(cross_avg((864/2-nm/2):(864/2+nm/2))) + mode=1 + npts=145 + nterms=5 + call polyfit(x,y,sigmay,npts,nterms,mode,a,chisqr) + pp=a(1)+x*(a(2)+x*(a(3)+x*(a(4)+x*a(5)))) + rmsdiff=sum( (pp-phase((864/2-nm/2):(864/2+nm/2)))**2 )/145.0 +write(*,*) 'training ',navg,sqrt(chisqr),rmsdiff + if( (sqrt(chisqr).lt.1.5) .and. (rmsdiff.lt.1.1) .and. (navg.ge.5) ) then + open(18,file='phasefit.dat',status='unknown') + do i=1, 145 + write(18,*) x(i),pp(i),y(i),sigmay(i) + enddo + close(18) + open(19,file='cframe.dat',status='unknown') + do i=1,864 + write(19,*) i,real(cframe(i)),imag(cframe(i)) + enddo + close(19) + pcoeffs=a + trained_dxcall=dxcall + training_dxcall(1:12)=' ' + currently_training=.false. + trained=.true. + call date_and_time(date,time,zone,values) + write(pcoeff_filename,'(i2.2,i2.2,i2.2,"_",i2.2,i2.2,i2.2)') & + values(1)-2000,values(2),values(3),values(5),values(6),values(7) + pcoeff_filename=trim(trained_dxcall)//"_"//pcoeff_filename//".pcoeff" +write(*,*) pcoeff_filename + open(17,file=pcoeff_filename,status='new') + write(17,'(i4,2f10.2,5f10.4)') navg,sqrt(chisqr),rmsdiff,a(1),a(2),a(3),a(4),a(5) + close(17) + endif + endif + + return + end subroutine msk144signalquality