From be724611429037e6e8ef7f53e7030de5bdd29f24 Mon Sep 17 00:00:00 2001 From: Steve Franke Date: Thu, 18 Apr 2019 14:15:24 -0500 Subject: [PATCH] FT4: Some groundwork for subtraction. --- CMakeLists.txt | 1 + lib/ft4/ft4sim_mult.f90 | 4 ++- lib/ft4/gen_ft4wave.f90 | 28 ++++++++++++----- lib/ft4/subtractft4.f90 | 66 +++++++++++++++++++++++++++++++++++++++++ lib/ft8/gen_ft8wave.f90 | 1 + widgets/mainwindow.cpp | 6 ++-- 6 files changed, 96 insertions(+), 10 deletions(-) create mode 100644 lib/ft4/subtractft4.f90 diff --git a/CMakeLists.txt b/CMakeLists.txt index 5d39c437d..3fbb879dc 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -555,6 +555,7 @@ set (wsjt_FSRCS lib/stdmsg.f90 lib/subtract65.f90 lib/ft8/subtractft8.f90 + lib/ft4/subtractft4.f90 lib/sun.f90 lib/symspec.f90 lib/symspec2.f90 diff --git a/lib/ft4/ft4sim_mult.f90 b/lib/ft4/ft4sim_mult.f90 index 065b2a76b..5acfaccef 100644 --- a/lib/ft4/ft4sim_mult.f90 +++ b/lib/ft4/ft4sim_mult.f90 @@ -10,6 +10,7 @@ program ft4sim_mult type(hdr) h !Header for .wav file character arg*12,fname*17,cjunk*4 character msg37*37,msgsent37*37,c77*77 + complex cwave0((NN+2)*NSPS) real wave0((NN+2)*NSPS) real wave(NZZ) real tmp(NZZ) @@ -63,7 +64,8 @@ program ft4sim_mult call pack77(msg37,i3,n3,c77) call genft4(msg37,0,msgsent37,itone) nwave0=(NN+2)*NSPS - call gen_ft4wave(itone,NN,NSPS,12000.0,f0,wave0,nwave0) + icmplx=0 + call gen_ft4wave(itone,NN,NSPS,12000.0,f0,cwave0,wave0,icmplx,nwave0) k0=nint((xdt+0.5)/dt) if(k0.lt.1) k0=1 diff --git a/lib/ft4/gen_ft4wave.f90 b/lib/ft4/gen_ft4wave.f90 index 305310e68..aa4174645 100644 --- a/lib/ft4/gen_ft4wave.f90 +++ b/lib/ft4/gen_ft4wave.f90 @@ -1,6 +1,7 @@ -subroutine gen_ft4wave(itone,nsym,nsps,fsample,f0,wave,nwave) +subroutine gen_ft4wave(itone,nsym,nsps,fsample,f0,cwave,wave,icmplx,nwave) real wave(nwave) + complex cwave(nwave) real pulse(6144) !512*4*3 real dphi(0:240000-1) integer itone(nsym) @@ -34,19 +35,32 @@ subroutine gen_ft4wave(itone,nsym,nsps,fsample,f0,wave,nwave) phi=0.0 dphi = dphi + twopi*f0*dt !Shift frequency up by f0 wave=0. + cwave=0. k=0 do j=0,nwave-1 k=k+1 - wave(k)=sin(phi) + if(icmplx.eq.0) then + wave(k)=sin(phi) + else + cwave(k)=cmplx(cos(phi),sin(phi)) + endif phi=mod(phi+dphi(j),twopi) enddo ! Compute the ramp-up and ramp-down symbols - wave(1:nsps)=wave(1:nsps) * & - (1.0-cos(twopi*(/(i,i=0,nsps-1)/)/(2.0*nsps)))/2.0 - k1=(nsym+1)*nsps+1 - wave(k1:k1+nsps-1)=wave(k1:k1+nsps-1) * & - (1.0+cos(twopi*(/(i,i=0,nsps-1)/)/(2.0*nsps)))/2.0 + if(icmplx.eq.0) then + wave(1:nsps)=wave(1:nsps) * & + (1.0-cos(twopi*(/(i,i=0,nsps-1)/)/(2.0*nsps)))/2.0 + k1=(nsym+1)*nsps+1 + wave(k1:k1+nsps-1)=wave(k1:k1+nsps-1) * & + (1.0+cos(twopi*(/(i,i=0,nsps-1)/)/(2.0*nsps)))/2.0 + else + cwave(1:nsps)=cwave(1:nsps) * & + (1.0-cos(twopi*(/(i,i=0,nsps-1)/)/(2.0*nsps)))/2.0 + k1=(nsym+1)*nsps+1 + cwave(k1:k1+nsps-1)=cwave(k1:k1+nsps-1) * & + (1.0+cos(twopi*(/(i,i=0,nsps-1)/)/(2.0*nsps)))/2.0 + endif return end subroutine gen_ft4wave diff --git a/lib/ft4/subtractft4.f90 b/lib/ft4/subtractft4.f90 new file mode 100644 index 000000000..16e3b9e7f --- /dev/null +++ b/lib/ft4/subtractft4.f90 @@ -0,0 +1,66 @@ +subroutine subtractft4(dd,itone,f0,dt) + +! Subtract an ft4 signal +! +! Measured signal : dd(t) = a(t)cos(2*pi*f0*t+theta(t)) +! Reference signal : cref(t) = exp( j*(2*pi*f0*t+phi(t)) ) +! Complex amp : cfilt(t) = LPF[ dd(t)*CONJG(cref(t)) ] +! Subtract : dd(t) = dd(t) - 2*REAL{cref*cfilt} + + use timer_module, only: timer + + parameter (NMAX=6*12000,NFRAME=(103+2)*79) + parameter (NFFT=NMAX,NFILT=1400) + real*4 dd(NMAX), window(-NFILT/2:NFILT/2), xjunk + complex cref,camp,cfilt,cw + integer itone(79) + logical first + data first/.true./ + common/heap8/cref(NFRAME),camp(NMAX),cfilt(NMAX),cw(NMAX),xjunk(NFRAME) + save first + + nstart=dt*12000+1 + nsym=79 + nsps=1920 + fs=12000.0 + icmplx=1 + bt=1.0 ! Temporary compromise? + call gen_ft4wave(itone,nsym,nsps,fs,f0,cref,xjunk,icmplx,NFRAME) + camp=0. + do i=1,nframe + id=nstart-1+i + if(id.ge.1.and.id.le.NMAX) camp(i)=dd(id)*conjg(cref(i)) + enddo + + if(first) then +! Create and normalize the filter + pi=4.0*atan(1.0) + fac=1.0/float(nfft) + sum=0.0 + do j=-NFILT/2,NFILT/2 + window(j)=cos(pi*j/NFILT)**2 + sum=sum+window(j) + enddo + cw=0. + cw(1:NFILT+1)=window/sum + cw=cshift(cw,NFILT/2+1) + call four2a(cw,nfft,1,-1,1) + cw=cw*fac + first=.false. + endif + + cfilt=0.0 + cfilt(1:nframe)=camp(1:nframe) + call four2a(cfilt,nfft,1,-1,1) + cfilt(1:nfft)=cfilt(1:nfft)*cw(1:nfft) + call four2a(cfilt,nfft,1,1,1) + +! Subtract the reconstructed signal + do i=1,nframe + j=nstart+i-1 + if(j.ge.1 .and. j.le.NMAX) dd(j)=dd(j)-2*REAL(cfilt(i)*cref(i)) + enddo + + return +end subroutine subtractft4 + diff --git a/lib/ft8/gen_ft8wave.f90 b/lib/ft8/gen_ft8wave.f90 index 550011301..37f3b0cbe 100644 --- a/lib/ft8/gen_ft8wave.f90 +++ b/lib/ft8/gen_ft8wave.f90 @@ -42,6 +42,7 @@ subroutine gen_ft8wave(itone,nsym,nsps,bt,fsample,f0,cwave,wave,icmplx,nwave) phi=0.0 dphi = dphi + twopi*f0*dt !Shift frequency up by f0 wave=0. + cwave=0. k=0 do j=nsps,nsps+nwave-1 !Don't include dummy symbols k=k+1 diff --git a/widgets/mainwindow.cpp b/widgets/mainwindow.cpp index 5d59abb6a..a78b12a7f 100644 --- a/widgets/mainwindow.cpp +++ b/widgets/mainwindow.cpp @@ -106,7 +106,7 @@ extern "C" { float xjunk[], float wave[], int* icmplx, int* nwave); void gen_ft4wave_(int itone[], int* nsym, int* nsps, float* fsample, float* f0, - float wave[], int* nwave); + float xjunk[], float wave[], int* icmplx, int* nwave); void gen4_(char* msg, int* ichk, char* msgsent, int itone[], int* itext, fortran_charlen_t, fortran_charlen_t); @@ -3745,7 +3745,9 @@ void MainWindow::guiUpdate() float fsample=48000.0; float f0=ui->TxFreqSpinBox->value() - m_XIT; int nwave=(nsym+2)*nsps; - gen_ft4wave_(const_cast(itone),&nsym,&nsps,&fsample,&f0,foxcom_.wave,&nwave); + int icmplx=0; + gen_ft4wave_(const_cast(itone),&nsym,&nsps,&fsample,&f0,foxcom_.wave, + foxcom_.wave,&icmplx,&nwave); } if(SpecOp::EU_VHF==m_config.special_op_id()) {