subroutine refspectrum(id2,id2b,kk,bclear,brefspec,buseref,fname) ! Input: ! id2 i*2 Raw 16-bit integer data, 12000 Hz sample rate ! brefspec logical True when accumulating a reference spectrum parameter (NFFT=6912,NH=NFFT/2,NPOLYLOW=400,NPOLYHIGH=2600) integer*2 id2(NFFT),id2b(120*12000) logical*1 bclear,brefspec,buseref,blastuse real x0(0:NH-1) !Input samples real x1(0:NH-1) !Output samples (delayed by one block) real x0s(0:NH-1) !Saved upper half of input samples real x1s(0:NH-1) !Saved upper half of output samples real x(0:NFFT-1) !Work array real*4 w(0:NFFT-1) !Window function real*4 s(0:NH) !Average spectrum real*4 fil(0:NH) real*8 xfit(1500),yfit(1500),sigmay(1500),a(5),chisqr !Polyfit arrays logical first complex cx(0:NH) !Complex frequency-domain work array character*(*) fname common/spectra/syellow(6827),ref(0:NH),filter(0:NH) equivalence(x,cx) data first/.true./,blastuse/.false./ save if(first) then pi=4.0*atan(1.0) do i=0,NFFT-1 ww=sin(i*pi/NFFT) w(i)=ww*ww/NFFT enddo nsave=0 s=0.0 filter=1.0 x0s=0. x1s=0. first=.false. endif if(bclear) s=0. if(brefspec) then x(0:NH-1)=0.001*id2(1:NH) x(NH:NFFT-1)=0.0 call four2a(x,NFFT,1,-1,0) !r2c FFT do i=1,NH s(i)=s(i) + real(cx(i))**2 + aimag(cx(i))**2 enddo nsave=nsave+1 fac0=0.9 if(mod(nsave,4).eq.0) then df=12000.0/NFFT ia=nint(1000.0/df) ib=nint(2000.0/df) avemid=sum(s(ia:ib))/(ib-ia+1) do i=0,NH fil(i)=0. if(s(i).gt.0.0) then fil(i)=sqrt(avemid/s(i)) endif enddo ! Default range is 240 - 4000 Hz. For narrower filters, use frequencies ! at which gain is -20 dB relative to 1500 Hz. ia=nint(240.0/df) ib=nint(4000.0/df) i0=nint(1500.0/df) do i=i0,ia,-1 if(s(i)/s(i0).lt.0.01) exit enddo ia=i do i=i0,ib,1 if(s(i)/s(i0).lt.0.01) exit enddo ib=i fac=fac0 do i=ia,1,-1 fac=fac*fac0 fil(i)=fac*fil(i) enddo fac=fac0 do i=ib,NH fac=fac*fac0 fil(i)=fac*fil(i) enddo do iter=1,100 !### ??? ### call smo121(fil,NH) enddo do i=0,NH filter(i)=-60.0 if(s(i).gt.0.0) filter(i)=20.0*log10(fil(i)) enddo il=nint(NPOLYLOW/df) ih=nint(NPOLYHIGH/df) nfit=ih-il+1 mode=0 nterms=5 do i=1,nfit xfit(i)=((i+il-1)*df-1500.0)/1000.0 yfit(i)=fil(i+il-1) sigmay(i)=1.0 enddo call polyfit(xfit,yfit,sigmay,nfit,nterms,mode,a,chisqr) open(16,file=fname,status='unknown') write(16,1003) NPOLYLOW,NPOLYHIGH,nterms,a 1003 format(3i5,5e25.16) do i=1,NH freq=i*df ref(i)=db(s(i)/avemid) write(16,1005) freq,s(i),ref(i),fil(i),filter(i) 1005 format(f10.3,e12.3,f12.6,e12.3,f12.6) enddo close(16) endif return endif if(buseref) then if(blastuse.neqv.buseref) then !just enabled so read filter fil=1.0 open(16,file=fname,status='old',err=110) read(16,1003,err=20,end=100) ndummy,ndummy,nterms,a goto 30 20 rewind(16) !allow for old style refspec.dat with no header 30 do i=1,NH read(16,1005,err=100,end=100) freq,s(i),ref(i),fil(i),filter(i) enddo 100 close(16) 110 continue endif x0=id2(1:NH) x(0:NH-1)=x0s !Previous 2nd half to new 1st half x(NH:NFFT-1)=x0 !New 2nd half x0s=x0 !Save the new 2nd half x=w*x !Apply window call four2a(x,NFFT,1,-1,0) !r2c FFT (to frequency domain) cx=fil*cx call four2a(cx,NFFT,1,1,-1) !c2r FFT (back to time domain) x1=x1s + x(0:NH-1) !Add previous segment's 2nd half ! id2(1:NH)=nint(x1) if(kk.ge.6912) id2b(kk-6192+1:kk-6192+NH)=nint(x1) x1s=x(NH:NFFT-1) !Save the new 2nd half endif blastuse=buseref return end subroutine refspectrum