From a92e7508b9f8a38c0b96f2288f9d76cd160f2e04 Mon Sep 17 00:00:00 2001 From: Steven Franke Date: Fri, 2 Jun 2017 13:58:32 +0000 Subject: [PATCH] wspr5d_exp now tries to detect sequences of 3, 6, and 9 bits. git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@7700 ab8295b8-cf94-4d9e-aec4-7959e3be5d79 --- lib/fsk4hf/bpdecode300.f90 | 7 +- lib/fsk4hf/getfc2w.f90 | 9 - lib/fsk4hf/osd300.f90 | 10 +- lib/fsk4hf/wspr5d.f90 | 5 +- lib/fsk4hf/wspr5d_exp.f90 | 437 ++++++++++++++++++++++++++--------- lib/fsk4hf/wspr5sim.f90 | 18 +- lib/fsk4hf/wspr_fsk8_sim.f90 | 14 +- 7 files changed, 363 insertions(+), 137 deletions(-) diff --git a/lib/fsk4hf/bpdecode300.f90 b/lib/fsk4hf/bpdecode300.f90 index f4374354e..22f7f885a 100644 --- a/lib/fsk4hf/bpdecode300.f90 +++ b/lib/fsk4hf/bpdecode300.f90 @@ -647,9 +647,14 @@ do iter=0,maxiterations enddo !write(*,*) 'number of unsatisfied parity checks ',ncheck if( ncheck .eq. 0 ) then ! we have a codeword - reorder the columns and return it - niterations=iter +! niterations=iter codeword=cw(colorder+1) decoded=codeword(M+1:N) + nerr=0 + do i=1,N + if( (2*cw(i)-1)*llr(i) .lt. 0.0 ) nerr=nerr+1 + enddo + niterations=nerr return endif diff --git a/lib/fsk4hf/getfc2w.f90 b/lib/fsk4hf/getfc2w.f90 index fb7a199bb..0d100ee56 100644 --- a/lib/fsk4hf/getfc2w.f90 +++ b/lib/fsk4hf/getfc2w.f90 @@ -20,16 +20,8 @@ subroutine getfc2w(c,csync,npeaks,fs,fc1,fpks) ia=nint(0.75*baud/df) cs(ia:NZ-1-ia)=0. !Save only freqs around fc1 -! do i=1,NZ/2 -! filt=1/(1+((i*df)**2/(0.50*baud)**2)**8) -! cs(i)=cs(i)*filt -! cs(NZ+1-i)=cs(NZ+1-i)*filt -! enddo call four2a(cs,NZ,1,1,1) !Back to time domain cs=cs/NZ -!do i=0,NZ-1 -!write(51,*) i,real(cs(i)),imag(cs(i)) -!enddo cs=cs*cs !Square the data call four2a(cs,NZ,1,-1,1) !Compute squared spectrum @@ -37,7 +29,6 @@ subroutine getfc2w(c,csync,npeaks,fs,fc1,fpks) pmax=0. fc2=0. ja=nint(0.3*baud/df) -! ja=nint(0.5*baud/df) k=1 do j=-ja,ja f2=j*df diff --git a/lib/fsk4hf/osd300.f90 b/lib/fsk4hf/osd300.f90 index 9ef2a8397..c4ec9a025 100644 --- a/lib/fsk4hf/osd300.f90 +++ b/lib/fsk4hf/osd300.f90 @@ -80,9 +80,6 @@ do id=1,K ! diagonal element indices enddo g2=transpose(genmrb) -!do i=1,N -! g2(i,1:K)=genmrb(1:K,i) -!enddo ! The hard decisions for the K MRB bits define the order 0 message, m0. ! Encode m0 using the modified generator matrix to find the "order 0" codeword. @@ -140,8 +137,13 @@ enddo ! re-order the codeword to place message bits at the end cw(indices)=cw +hdec(indices)=hdec decoded=cw(M+1:N) -niterations=1 +nerr=0 +do i=1,N + if( hdec(i) .ne. cw(i) ) nerr=nerr+1 +enddo +niterations=nerr return end subroutine osd300 diff --git a/lib/fsk4hf/wspr5d.f90 b/lib/fsk4hf/wspr5d.f90 index 3d18dd629..7b50fff2f 100644 --- a/lib/fsk4hf/wspr5d.f90 +++ b/lib/fsk4hf/wspr5d.f90 @@ -122,10 +122,10 @@ program wspr5d go to 999 endif close(10) - fa=102.0 + fa=100.0 fb=150.0 call getfc1w(c,fs,fa,fb,fc1,xsnr) !First approx for freq -npeaks=20 + npeaks=20 call getfc2w(c,csync,npeaks,fs,fc1,fpks) !Refined freq a(1)=-fc1 @@ -152,6 +152,7 @@ npeaks=20 ibb=NZ-1-j endif z=sum(c(ia:ib)*conjg(csync(iaa:ibb))) +write(51,*) j/fs,real(z),imag(z) if(abs(z).gt.amax) then amax=abs(z) jpk=j diff --git a/lib/fsk4hf/wspr5d_exp.f90 b/lib/fsk4hf/wspr5d_exp.f90 index 781ee17ec..7ce456f48 100644 --- a/lib/fsk4hf/wspr5d_exp.f90 +++ b/lib/fsk4hf/wspr5d_exp.f90 @@ -7,18 +7,33 @@ program wspr5d ! ! Still to do: find and decode more than one signal in the specified passband. - include 'wsprlf_params.f90' +! include 'wsprlf_params.f90' + + parameter (NDOWN=30) + parameter (KK=60) + parameter (ND=300) + parameter (NS=109) + parameter (NR=3) + parameter (NN=NR+NS+ND) + parameter (NSPS0=8640) + parameter (NSPS=16) + parameter (N2=2*NSPS) + parameter (NZ=NSPS*NN) + parameter (NZ400=288*NN) parameter (NMAX=300*12000) + character arg*8,message*22,cbits*50,infile*80,fname*16,datetime*11 character*120 data_dir complex csync(0:NZ-1) !Sync symbols only, from cbb + complex c400(0:NZ400-1) !Complex waveform complex c(0:NZ-1) !Complex waveform - complex cd(0:412*16-1) !Complex waveform - complex ca(0:412*16-1) !Complex waveform + complex cd(0:NZ-1) !Complex waveform + complex ca(0:NZ-1) !Complex waveform + complex zz real*8 fMHz real rxdata(ND),llr(ND) !Soft symbols real pp(32) !Shaped pulse for OQPSK - real ps(0:7),sbits(412) + real sbits(412),softbits(9) real fpks(20) integer id(NS+ND) !NRZ values (+/-1) for Sync and Data integer isync(48) !Long sync vector @@ -28,7 +43,7 @@ program wspr5d integer*2 iwave(NMAX) !Generated full-length waveform integer*1 idat(7) integer*1 decoded(KK),apmask(ND),cw(ND) - integer*1 hbits(412),ebits(411),bits(5) + integer*1 hbits(412),bits(13) data ib13/1,1,1,1,1,-1,-1,1,1,-1,1,-1,1/ nargs=iargc() @@ -98,117 +113,151 @@ program wspr5d j1=index(infile,'.c5') j2=index(infile,'.wav') if(j1.gt.0) then - read(10,end=999) fname,ntrmin,fMHz,c + read(10,end=999) fname,ntrmin,fMHz,c400 read(fname(8:11),*) nutc write(datetime,'(i11)') nutc else if(j2.gt.0) then read(10,end=999) ihdr,iwave read(infile(j2-4:j2-1),*) nutc datetime=infile(j2-11:j2-1) - call wspr5_downsample(iwave,c) + call wspr5_downsample(iwave,c400) else print*,'Wrong file format?' go to 999 endif - close(10) + + fa=100.0 - fb=170.0 - call getfc1w(c,fs,fa,fb,fc1,xsnr) !First approx for freq -! call getfc2w(c,csync,fs,fc1,fc2,fc3) !Refined freq + fb=150.0 + fs400=400.0 + call getfc1(c400,fs400,fa,fb,fc1,xsnr) !First approx for freq +!write(*,*) datetime,'initial guess ',fc1 npeaks=5 - call getfc2w(c,csync,npeaks,fs,fc1,fpks) !Refined freq + call getfc2(c400,npeaks,fs400,fc1,fpks) !Refined freq + do idf=1,npeaks ! consider the top npeak peaks + fc2=fpks(idf) +!write(*,*) 'peak ',idf,fc1+fc2,fc2 + call downsample(c400,fc1+fc2,cd) + s2=sum(cd*conjg(cd))/(16*412) + cd=cd/sqrt(s2) + + do is=0,11 ! search over plus/minus 0.25 seconds for now + idt=is/2 + if( mod(is,2).eq. 1 ) idt=-(is+1)/2 + xdt=real(22+idt)/22.222 - 1.0 + ca=cshift(cd,22+idt) + do iseq=1,3 ! try sequence estimation lengths of 3, 6, and 9 bits. + k=1-2*iseq + nseq=iseq*3 + do i=1,408,iseq*4 + k=k+iseq*2 + j=(i+1)*16 + call mskseqdet(nseq,ca(j),pp,id(k),softbits,1,phase) + hbits(i:i+iseq*4)=bits + sbits(i:i+iseq*4)=bits -!write(*,*) fc1+fc2 -do ipks=1,npeaks - call downsample(c,fc1+fpks(ipks),cd) + sbits(i+1)=softbits(1) + sbits(i+2)=softbits(2) + if( id(k+1) .ne. 0 ) sbits(i+2)=id(k+1)*25 + sbits(i+3)=softbits(3) - do ncoh=1,1,-1 - do is=0,9 - idt=is/2 - if( mod(is,2).eq. 1 ) idt=-is/2 - xdt=idt/22.222 - k=-1 - ca=cshift(cd,22+idt) - do i=1,408,4 - k=k+2 - j=(i+1)*16 - call mskseqdet(ca(j),pp,id(k),bits,ps,ncoh) - r1=max(ps(1),ps(3),ps(5),ps(7))-max(ps(0),ps(2),ps(4),ps(6)) - r2=max(ps(2),ps(3),ps(6),ps(7))-max(ps(0),ps(1),ps(4),ps(5)) - r4=max(ps(4),ps(5),ps(6),ps(7))-max(ps(0),ps(1),ps(2),ps(3)) - hbits(i:i+4)=bits - sbits(i:i+4)=bits - sbits(i+1)=r4 - sbits(i+2)=r2 - if( id(k+1) .ne. 0 ) sbits(i+2)=id(k+1)*25 - sbits(i+3)=r1 + if( iseq .ge. 2 ) then + sbits(i+5)=softbits(4) + sbits(i+6)=softbits(5) + if( id(k+3) .ne. 0 ) sbits(i+6)=id(k+3)*25 + sbits(i+7)=softbits(6) + if( iseq .eq. 3 ) then + sbits(i+9)=softbits(7) + sbits(i+10)=softbits(8) + if( id(k+5) .ne. 0 ) sbits(i+10)=id(k+5)*25 + sbits(i+11)=softbits(9) + endif + endif + enddo + j=1 + do i=1,205 + if( abs(id(i)) .ne. 2 ) then + rxdata(j)=sbits(2*i-1) + j=j+1 + endif + enddo + do i=1,204 + rxdata(j)=sbits(2*i) + j=j+1 + enddo + rxav=sum(rxdata)/ND + rx2av=sum(rxdata*rxdata)/ND + rxsig=sqrt(rx2av-rxav*rxav) + rxdata=rxdata/rxsig +! sigma=0.84 + sigma=1.20 + llr=2*rxdata/(sigma*sigma) + apmask=0 + max_iterations=40 + ifer=0 + nbadcrc=0 + call bpdecode300(llr,apmask,max_iterations,decoded,niterations,cw) +! niterations will be equal to the Hamming distance between hard received word and the codeword + if(niterations.lt.0) call osd300(llr,3,decoded,niterations,cw) + if(niterations.ge.0) call chkcrc10(decoded,nbadcrc) + if(niterations.lt.0 .or. nbadcrc.ne.0) ifer=1 + if( ifer.eq.0 ) then + write(cbits,1200) decoded(1:50) +1200 format(50i1) + read(cbits,1202) idat +1202 format(6b8,b2) + idat(7)=ishft(idat(7),6) + call wqdecode(idat,message,itype) + nsnr=nint(xsnr) + freq=fMHz + 1.d-6*(fc1+fc2) + nfdot=0 + write(13,1210) datetime,0,nsnr,xdt,freq,message,nfdot +1210 format(a11,2i4,f6.2,f12.7,2x,a22,i3) + write(*,1212) datetime(8:11),nsnr,xdt,freq,nfdot,message,'*',idf,nseq,is,niterations +1212 format(a4,i4,f5.1,f11.6,i3,2x,a22,a1,i3,i3,i3,i4) + goto 888 + endif + enddo !iseq enddo - - j=1 - do i=1,205 - if( abs(id(i)) .ne. 2 ) then - rxdata(j)=sbits(2*i-1) - j=j+1 - endif - enddo - do i=1,204 - rxdata(j)=sbits(2*i) - j=j+1 - enddo - rxav=sum(rxdata)/ND - rx2av=sum(rxdata*rxdata)/ND - rxsig=sqrt(rx2av-rxav*rxav) - rxdata=rxdata/rxsig - sigma=0.84 - llr=2*rxdata/(sigma*sigma) - apmask=0 - max_iterations=40 - ifer=0 - nbadcrc=0 - call bpdecode300(llr,apmask,max_iterations,decoded,niterations,cw) - if(niterations.lt.0) call osd300(llr,3,decoded,niterations,cw) - if(niterations.ge.0) call chkcrc10(decoded,nbadcrc) - if(niterations.lt.0 .or. nbadcrc.ne.0) ifer=1 - if( ifer.eq.0 ) then - write(cbits,1200) decoded(1:50) -1200 format(50i1) - read(cbits,1202) idat -1202 format(6b8,b2) - idat(7)=ishft(idat(7),6) - call wqdecode(idat,message,itype) - nsnr=nint(xsnr) - freq=fMHz + 1.d-6*(fc1+fc2) - nfdot=0 - write(13,1210) datetime,0,nsnr,xdt,freq,message,nfdot -1210 format(a11,2i4,f6.2,f12.7,2x,a22,i3) - write(*,1212) datetime(8:11),nsnr,xdt,freq,nfdot,message,'*',ipks -1212 format(a4,i4,f5.1,f11.6,i3,2x,a22,a1,i4) - goto 888 - endif enddo -enddo -enddo ! fpeaks loop -888 enddo +888 continue + enddo write(*,1120) 1120 format("") 999 end program wspr5d -subroutine mskseqdet(cdat,pp,bsync,bestbits,cmbest,ncoh) -complex cdat(16*4),cbest(16*4),cideal(16*4) -complex cdf(16*4),cfac -real cm(0:7),cmbest(0:7) -real pp(32) -integer*1 bits(5),bestbits(5),sgn(5) -integer bsync(3) +subroutine getmetric(ib,ps,xmet) + real ps(0:511) + xm1=0 + xm0=0 + do i=0,511 + if( iand(i/ib,1) .eq. 1 .and. ps(i) .gt. xm1 ) xm1=ps(i) + if( iand(i/ib,1) .eq. 0 .and. ps(i) .gt. xm0 ) xm0=ps(i) + enddo + xmet=xm1-xm0 + return +end subroutine getmetric + +subroutine mskseqdet(ns,cdat,pp,bsync,softbits,ncoh,phase) +! +! Detect sequences of 3, 6, or 9 bits (ns). +! Sync bits are assumed to be known. +! +complex cdat(16*12),cbest(16*12),cideal(16*12) +complex cdf(16*12),cfac,zz +real cm(0:511),cmbest(0:511) +real pp(32),softbits(9) +integer bit(13),bestbits(13),sgn(13) +integer bsync(7) twopi=8.0*atan(1.0) dt=30.0*18.0/12000.0 cmax=0; fbest=0.0; - +np=2**ns-1 idfmax=40 if( ncoh .eq. 1 ) idfmax=0 do idf=0,idfmax @@ -217,40 +266,88 @@ do idf=0,idfmax dphi=twopi*deltaf*dt cfac=cmplx(cos(dphi),sin(dphi)) cdf=1.0 - do i=2,16*4 + do i=2,16*(ns-1) cdf(i)=cdf(i-1)*cfac enddo cm=0 ibflag=0 - do i=0,7 - bits(1)=(bsync(1)+2)/4 - bits(2)=iand(i/4,1) - bits(3)=iand(i/2,1) + do i=0,np + bit(1)=(bsync(1)+2)/4 + bit(2)=iand(i/(2**(ns-1)),1) + bit(3)=iand(i/(2**(ns-2)),1) if( bsync(2).ne.0 ) then ! force the barker bits - bits(3)=(bsync(2)+2)/4 + bit(3)=(bsync(2)+2)/4 + endif + bit(4)=iand(i/(2**(ns-3)),1) + bit(5)=(bsync(3)+2)/4 + + if( ns .ge. 6 ) then + bit(6)=iand(i/(2**(ns-4)),1) + bit(7)=iand(i/(2**(ns-5)),1) + if( bsync(4).ne.0 ) then ! force the barker bits + bit(7)=(bsync(4)+2)/4 + endif + bit(8)=iand(i/(2**(ns-6)),1) + bit(9)=(bsync(5)+2)/4 + if( ns .eq. 9 ) then + bit(10)=iand(i/4,1) + bit(11)=iand(i/2,1) + if( bsync(6).ne.0 ) then ! force the barker bits + bit(11)=(bsync(6)+2)/4 + endif + bit(12)=iand(i/1,1) + bit(13)=(bsync(7)+2)/4 + endif + endif + + sgn=2*bit-1 + cideal(1:16) =cmplx(sgn(1)*pp(17:32),sgn(2)*pp(1:16)) + cideal(17:32) =cmplx(sgn(3)*pp(1:16),sgn(2)*pp(17:32)) + cideal(33:48) =cmplx(sgn(3)*pp(17:32),sgn(4)*pp(1:16)) + cideal(49:64) =cmplx(sgn(5)*pp(1:16),sgn(4)*pp(17:32)) + if( ns .ge. 6 ) then + cideal(65:80) =cmplx(sgn(5)*pp(17:32),sgn(6)*pp(1:16)) + cideal(81:96) =cmplx(sgn(7)*pp(1:16),sgn(6)*pp(17:32)) + cideal(97:112) =cmplx(sgn(7)*pp(17:32),sgn(8)*pp(1:16)) + cideal(113:128)=cmplx(sgn(9)*pp(1:16),sgn(8)*pp(17:32)) + if( ns .eq. 9 ) then + cideal(129:144) =cmplx(sgn(9)*pp(17:32),sgn(10)*pp(1:16)) + cideal(145:160) =cmplx(sgn(11)*pp(1:16),sgn(10)*pp(17:32)) + cideal(161:176) =cmplx(sgn(11)*pp(17:32),sgn(12)*pp(1:16)) + cideal(177:192)=cmplx(sgn(13)*pp(1:16),sgn(12)*pp(17:32)) + endif endif - bits(4)=iand(i/1,1) - bits(5)=(bsync(3)+2)/4 - sgn=2*bits-1 - cideal(1:16)=cmplx(sgn(1)*pp(17:32),sgn(2)*pp(1:16)) - cideal(17:32)=cmplx(sgn(3)*pp(1:16),sgn(2)*pp(17:32)) - cideal(33:48)=cmplx(sgn(3)*pp(17:32),sgn(4)*pp(1:16)) - cideal(49:64)=cmplx(sgn(5)*pp(1:16),sgn(4)*pp(17:32)) cideal=cideal*cdf - cm(i)=abs(sum(cdat*conjg(cideal)))/1.e3 + cm(i)=abs(sum(cdat(1:64*ns/3)*conjg(cideal(1:64*ns/3))))/1.e3 if( cm(i) .gt. cmax ) then ibflag=1 cmax=cm(i) - bestbits=bits + bestbits=bit cbest=cideal fbest=deltaf + zz=sum(cdat*conjg(cbest))/1.e3 + phase=atan2(imag(zz),real(zz)) endif enddo if( ibflag .eq. 1 ) then ! new best found cmbest=cm endif enddo +softbits=0.0 +call getmetric(1,cmbest,softbits(ns)) +call getmetric(2,cmbest,softbits(ns-1)) +call getmetric(4,cmbest,softbits(ns-2)) +if( ns .ge. 6 ) then + call getmetric(8,cmbest,softbits(ns-3)) + call getmetric(16,cmbest,softbits(ns-4)) + call getmetric(32,cmbest,softbits(ns-5)) + if( ns .eq. 9 ) then + call getmetric(64,cmbest,softbits(3)) + call getmetric(128,cmbest,softbits(2)) + call getmetric(256,cmbest,softbits(1)) + endif +endif end subroutine mskseqdet subroutine downsample(ci,f0,co) @@ -260,11 +357,11 @@ subroutine downsample(ci,f0,co) df=400.0/NI ct=ci - call four2a(ct,NI,1,-1,1) !r2c FFT to freq domain + call four2a(ct,NI,1,-1,1) !c2c FFT to freq domain i0=nint(f0/df) co=0.0 co(0)=ct(i0) - b=4.0 + b=6.0 do i=1,NO/2 arg=(i*df/b)**2 filt=exp(-arg) @@ -275,3 +372,133 @@ subroutine downsample(ci,f0,co) call four2a(co,NO,1,1,1) !c2c FFT back to time domain return end subroutine downsample + +subroutine getfc1(c,fs,fa,fb,fc1,xsnr) + +! include 'wsprlf_params.f90' + parameter (NZ=288*412) + parameter (NSPS=288) + parameter (N2=2*NSPS) + parameter (NFFT1=16*NSPS) + parameter (NH1=NFFT1/2) + + complex c(0:NZ-1) !Complex waveform + complex c2(0:NFFT1-1) !Short spectra + real s(-NH1+1:NH1) !Coarse spectrum + nspec=NZ/N2 + df1=fs/NFFT1 + s=0. + do k=1,nspec + ia=(k-1)*N2 + ib=ia+N2-1 + c2(0:N2-1)=c(ia:ib) + c2(N2:)=0. + call four2a(c2,NFFT1,1,-1,1) + do i=0,NFFT1-1 + j=i + if(j.gt.NH1) j=j-NFFT1 + s(j)=s(j) + real(c2(i))**2 + aimag(c2(i))**2 + enddo + enddo +! call smo121(s,NFFT1) + smax=0. + ipk=0 + fc1=0. + ia=nint(fa/df1) + ib=nint(fb/df1) + do i=ia,ib + f=i*df1 + if(s(i).gt.smax) then + smax=s(i) + ipk=i + fc1=f + endif +! write(51,3001) f,s(i),db(s(i)) +! 3001 format(f10.3,e12.3,f10.3) + enddo + +! The following is for testing SNR calibration: + sp3n=(s(ipk-1)+s(ipk)+s(ipk+1)) !Sig + 3*noise + base=(sum(s)-sp3n)/(NFFT1-3.0) !Noise per bin + psig=sp3n-3*base !Sig only + pnoise=(2500.0/df1)*base !Noise in 2500 Hz + xsnr=db(psig/pnoise) + xsnr=xsnr+5.0 + return +end subroutine getfc1 + +subroutine getfc2(c,npeaks,fs,fc1,fpks) + +! include 'wsprlf_params.f90' + parameter (NZ=288*412) + parameter (NSPS=288) + parameter (N2=2*NSPS) + parameter (NFFT1=16*NSPS) + parameter (NH1=NFFT1/2) + + complex c(0:NZ-1) !Complex waveform + complex cs(0:NZ-1) !For computing spectrum + real a(5) + real freqs(413),sp2(413),fpks(npeaks) + integer pkloc(1) + + df=fs/NZ + baud=fs/NSPS + a(1)=-fc1 + a(2:5)=0. + call twkfreq1(c,NZ,fs,a,cs) !Mix down by fc1 + +! Filter, square, then FFT to get refined carrier frequency fc2. + call four2a(cs,NZ,1,-1,1) !To freq domain + + ia=nint(0.75*baud/df) + cs(ia:NZ-1-ia)=0. !Save only freqs around fc1 +! do i=1,NZ/2 +! filt=1/(1+((i*df)**2/(0.50*baud)**2)**8) +! cs(i)=cs(i)*filt +! cs(NZ+1-i)=cs(NZ+1-i)*filt +! enddo + call four2a(cs,NZ,1,1,1) !Back to time domain + cs=cs/NZ +!do i=0,NZ-1 +!write(51,*) i,real(cs(i)),imag(cs(i)) +!enddo + cs=cs*cs !Square the data + call four2a(cs,NZ,1,-1,1) !Compute squared spectrum +! Find two peaks separated by baud + pmax=0. + fc2=0. +! ja=nint(0.3*baud/df) + ja=nint(0.5*baud/df) + k=1 + sp2=0.0 + do j=-ja,ja + f2=j*df + ia=nint((f2-0.5*baud)/df) + if(ia.lt.0) ia=ia+NZ + ib=nint((f2+0.5*baud)/df) + p=real(cs(ia))**2 + aimag(cs(ia))**2 + & + real(cs(ib))**2 + aimag(cs(ib))**2 + if(p.gt.pmax) then + pmax=p + fc2=0.5*f2 + endif + freqs(k)=0.5*f2 + sp2(k)=p + k=k+1 +! write(52,1200) f2,p,db(p) +!1200 format(f10.3,2f15.3) + enddo + + do i=1,npeaks + pkloc=maxloc(sp2) + ipk=pkloc(1) + fpks(i)=freqs(ipk) + ipk0=max(1,ipk-2) + ipk1=min(413,ipk+2) +! ipk0=ipk +! ipk1=ipk + sp2(ipk0:ipk1)=0.0 + enddo + return +end subroutine getfc2 diff --git a/lib/fsk4hf/wspr5sim.f90 b/lib/fsk4hf/wspr5sim.f90 index 8741eb322..5c2147b96 100644 --- a/lib/fsk4hf/wspr5sim.f90 +++ b/lib/fsk4hf/wspr5sim.f90 @@ -51,9 +51,9 @@ program wspr5sim txt=NN*NSPS0/12000.0 call genwspr5(msg,msgsent,itone) !Encode the message, get itone - write(*,1000) f0,xdt,txt,snrdb,nfiles,msgsent + write(*,1000) f0,xdt,txt,snrdb,fspread,delay,nfiles,msgsent 1000 format('f0:',f9.3,' DT:',f6.2,' txt:',f6.1,' SNR:',f6.1, & - ' nfiles:',i3,2x,a22) + ' fspread:',f6.1,' delay:',f6.1,' nfiles:',i3,2x,a22) dphi0=twopi*(f0-0.25*baud)*dt dphi1=twopi*(f0+0.25*baud)*dt @@ -73,19 +73,19 @@ program wspr5sim enddo enddo - if( fspread .ne. 0.0 .or. delay .ne. 0.0 ) then - call watterson(c0,NZ,fs,delay,fspread) - endif - - c0=sig*c0 !Scale to requested sig level - + call sgran() do ifile=1,nfiles + c=c0 if(nwav.eq.0) then + if( fspread .ne. 0.0 .or. delay .ne. 0.0 ) then + call watterson(c,NZ,fs,delay,fspread) + endif + c=c*sig if(snrdb.lt.90) then do i=0,NZ-1 !Add gaussian noise at specified SNR xnoise=gran() ynoise=gran() - c(i)=c0(i) + cmplx(xnoise,ynoise) + c(i)=c(i) + cmplx(xnoise,ynoise) enddo endif write(fname,1100) ifile diff --git a/lib/fsk4hf/wspr_fsk8_sim.f90 b/lib/fsk4hf/wspr_fsk8_sim.f90 index b4003f8da..40d808a0b 100644 --- a/lib/fsk4hf/wspr_fsk8_sim.f90 +++ b/lib/fsk4hf/wspr_fsk8_sim.f90 @@ -70,19 +70,19 @@ program wspr_fsk8_sim enddo enddo - if( fspread .ne. 0.0 .or. delay .ne. 0.0 ) then - call watterson(c0,NZ,fs,delay,fspread) - endif - - c0=sig*c0 !Scale to requested sig level - + call sgran() do ifile=1,nfiles if(nwav.eq.0) then + c=c0 + if( fspread .ne. 0.0 .or. delay .ne. 0.0 ) then + call watterson(c,NZ,fs,delay,fspread) + endif + c=c*sig if( snrdb.lt.90) then do i=0,NZ-1 xnoise=gran() ynoise=gran() - c(i)=c0(i)+cmplx(xnoise,ynoise) + c(i)=c(i)+cmplx(xnoise,ynoise) enddo endif write(fname,1100) ifile