Updates to remove_tone() routine. Currently disabled pending more tests.

This commit is contained in:
Steven Franke 2024-09-18 13:18:07 -05:00
parent 6b98c85473
commit acbc3095e9
2 changed files with 45 additions and 19 deletions

View File

@ -149,16 +149,17 @@ subroutine smo121a(x,nz,a,b)
return return
end subroutine smo121a end subroutine smo121a
subroutine remove_tone(c0,npts) subroutine remove_tone(c0,fsync)
parameter (NMAX=15*12000)
parameter (NFILT=8000) parameter (NFILT=8000)
complex c0(npts) complex c0(NMAX)
complex cwindow(15*12000) complex cwindow(15*12000)
complex cref(npts) complex cref(NMAX)
complex cfilt(npts) complex cfilt(NMAX)
real window(-NFILT/2:NFILT/2) real window(-NFILT/2:NFILT/2)
! real endcorrection(NFILT/2+1) ! real endcorrection(NFILT/2+1)
real s(npts/4) real s(NMAX/4)
integer ipk(1) integer ipk(1)
logical first logical first
data first/.true./ data first/.true./
@ -166,7 +167,7 @@ subroutine remove_tone(c0,npts)
if(first) then if(first) then
pi=4.0*atan(1.0) pi=4.0*atan(1.0)
fac=1.0/float(npts) fac=1.0/float(NMAX)
sumw=0.0 sumw=0.0
do j=-NFILT/2,NFILT/2 do j=-NFILT/2,NFILT/2
window(j)=cos(pi*j/NFILT)**2 window(j)=cos(pi*j/NFILT)**2
@ -175,37 +176,60 @@ subroutine remove_tone(c0,npts)
cwindow=0. cwindow=0.
cwindow(1:NFILT+1)=window/sumw cwindow(1:NFILT+1)=window/sumw
cwindow=cshift(cwindow,NFILT/2+1) cwindow=cshift(cwindow,NFILT/2+1)
call four2a(cwindow,npts,1,-1,1) call four2a(cwindow,NMAX,1,-1,1)
cwindow=cwindow*fac ! frequency domain smoothing filter cwindow=cwindow*fac ! frequency domain smoothing filter
first=.false. first=.false.
endif endif
fsample=12000.0 fsample=12000.0
df=fsample/npts df=fsample/NMAX
fac=1.0/npts fac=1.0/NMAX
cfilt=fac*c0 cfilt=fac*c0
call four2a(cfilt,npts,1,-1,1) ! fourier transform of input data call four2a(cfilt,NMAX,1,-1,1) ! fourier transform of input data
iz=npts/4 iz=NMAX/4
do i=1,iz do i=1,iz
s(i)=real(cfilt(i))**2 + aimag(cfilt(i))**2 s(i)=real(cfilt(i))**2 + aimag(cfilt(i))**2
enddo enddo
ia=nint(700.0/df) ia=nint((fsync-100.0)/df)
ib=nint(800.0/df) ib=nint((fsync+100.0)/df)
ipk=maxloc(s(ia:ib)) ipk=maxloc(s(ia:ib))
i0=ipk(1) + ia - 1 i0=ipk(1) + ia - 1
f2=df*i0
write(*,*) 'remove_tone - frequency: ',f2 i10=nint(10.0/df)
ia=i0-i10
ib=i0+i10
s0=0.0
s1=0.0
s2=0.0
do i=ia,ib
s0=s0+s(i)
s1=s1+(i-i0)*s(i)
enddo
delta=s1/s0
i0=nint(i0+delta)
f2=i0*df
ia=i0-i10
ib=i0+i10
do i=ia,ib
s2=s2 + s(i)*(i-i0)**2
enddo
sigma=sqrt(s2/s0)*df
! write(61,*) 'frequency, spectral width ',f2,sigma
if(sigma .gt. 2.5) go to 999
! write(61,*) 'remove_tone - frequency: ',f2
dt=1.0/fsample dt=1.0/fsample
do i=1, npts do i=1, NMAX
arg=2*pi*f2*i*dt arg=2*pi*f2*i*dt
cref(i)=cmplx(cos(arg),sin(arg)) cref(i)=cmplx(cos(arg),sin(arg))
enddo enddo
cfilt=c0*conjg(cref) ! baseband to be filtered cfilt=c0*conjg(cref) ! baseband to be filtered
call four2a(cfilt,npts,1,-1,1) call four2a(cfilt,NMAX,1,-1,1)
cfilt=cfilt*cwindow cfilt=cfilt*cwindow
call four2a(cfilt,npts,1,1,1) call four2a(cfilt,NMAX,1,1,1)
nframe=50*3456 nframe=50*3456
do i=1,nframe do i=1,nframe
@ -213,6 +237,6 @@ subroutine remove_tone(c0,npts)
c0(i)=c0(i)-cref(i) c0(i)=c0(i)-cref(i)
enddo enddo
return 999 return
end subroutine remove_tone end subroutine remove_tone

View File

@ -36,6 +36,8 @@ subroutine sfrx_sub(nyymmdd,nutc,nfqso,ntol,iwave)
call sfox_ana(dd,npts,c0,npts) call sfox_ana(dd,npts,c0,npts)
! call remove_tone(c0,fsync) ! Needs testing
ndepth=3 ndepth=3
dth=0.5 dth=0.5
damp=1.0 damp=1.0