subroutine apply_smooth(sv,sn,se,nfft,df,nf1,fre,smv,smn, # sme,nf,n_win,lw,jmin,jmax,ier) c c spectrum smoothing c a.tento 07/02 v.1.0 c a.tento 05/03 v.1.1 c a.tento 05/03 v.1.2 c a.tento 05/03 v.2.0 c implicit none include 'main.h' include 'parameters.h' real sv,sn,se, df, fre, smv, smn, sme integer nf1, nf, ier, jmin, jmax, n_win, lw, nfft dimension fre(nf), sv(nfft), sn(nfft), se(nfft), # smv(nf,n_win),smn(nf,n_win),sme(nf,n_win) C real band, sweight, weight, p, c1, c2, bk, c3, bb real fmax, bkmin, bkmax parameter ( bkmin = 3, bkmax = 100 ) integer nw, j, i1, i2, i character sep2*1 parameter (sep2 = ':') character*200 argout(20) c c fmax = float(nf1 - 1) * df +df/10.0 if (fre(nf) .gt. fmax) goto 97 if (fre(1) .lt. df ) goto 97 jmin = 1 jmax = nf c call split3(smooth,argout,sep2,nw) c ier = 1 if ( argout(1) .eq. 'none' )then call split3(freq_spacing,argout,sep2,nw) if( argout(1) .eq. 'fft' ) then do 7 j = 1, nf smv(j,lw) = sv(j+1) smn(j,lw) = sn(j+1) sme(j,lw) = se(j+1) 7 continue ier = 0 c else call interv(fre,nf,df,sv,nf1,smv,n_win,lw) call interv(fre,nf,df,sn,nf1,smn,n_win,lw) call interv(fre,nf,df,se,nf1,sme,n_win,lw) ier = 0 endif endif c if ( argout(1) .eq. 'linear' )then if ( argout(3) .eq. 'box' ) then ier = 0 read(argout(2),*) band if(band .lt. df) then call interv(fre,nf,df,sv,nf1,smv,n_win,lw) call interv(fre,nf,df,sn,nf1,smn,n_win,lw) call interv(fre,nf,df,se,nf1,sme,n_win,lw) else band = band/2. do 2 j = 1, nf i1 = int((fre(j) - band)/df + 0.5) + 1 i2 = int((fre(j) + band)/df + 0.5) + 1 if ( i1 .lt. 1) then jmin = j + 1 ier = 22 elseif(i2 .gt. nf1 ) then jmax = j - 1 ier = 22 goto 52 else smv(j,lw) = 0.0 smn(j,lw) = 0.0 sme(j,lw) = 0.0 do 3 i = i1, i2 smv(j,lw) = sv(i) + smv(j,lw) smn(j,lw) = sn(i) + smn(j,lw) sme(j,lw) = se(i) + sme(j,lw) 3 continue smv(j,lw) = smv(j,lw)/float(i2-i1 +1) smn(j,lw) = smn(j,lw)/float(i2-i1 +1) sme(j,lw) = sme(j,lw)/float(i2-i1 +1) endif 2 continue 52 continue endif endif if ( argout(3) .eq. 'tri' ) then ier = 0 read(argout(2),*) band band = band/2. if(band .lt. df) then call interv(fre,nf,df,sv,nf1,smv,n_win,lw) call interv(fre,nf,df,sn,nf1,smn,n_win,lw) call interv(fre,nf,df,se,nf1,sme,n_win,lw) else do 4 j = 1, nf i1 = int( (fre(j) - band)/df ) + 2 i2 = int( (fre(j) + band)/df ) + 1 if ( i1 .lt. 1) then jmin = j + 1 ier = 22 elseif(i2 .gt. nf1 ) then jmax = j - 1 ier = 22 goto 54 else smv(j,lw) = 0.0 smn(j,lw) = 0.0 sme(j,lw) = 0.0 sweight = 0.0 do 5 i = i1, i2 weight = 1.0 - abs(fre(j) - float(i-1)*df)/band smv(j,lw) = weight*sv(i) + smv(j,lw) smn(j,lw) = weight*sn(i) + smn(j,lw) sme(j,lw) = weight*se(i) + sme(j,lw) sweight = sweight + weight 5 continue smv(j,lw) = smv(j,lw)/sweight smn(j,lw) = smn(j,lw)/sweight sme(j,lw) = sme(j,lw)/sweight endif 4 continue 54 continue endif endif endif c if ( argout(1) .eq. 'log' )then c df <= fre(j) <= fny by deffre if ( argout(3) .eq. 'box' ) then ier = 0 read(argout(2),*) band p= 1 + band/100. do 12 j = 1, nf bb = fre(j) * (p - 1.0/p) if ( bb .lt. df ) then call interv1(j,fre,nf,df,sv,nf1,smv,n_win,lw) call interv1(j,fre,nf,df,sn,nf1,smn,n_win,lw) call interv1(j,fre,nf,df,se,nf1,sme,n_win,lw) else i1 = int((fre(j)/p)/df + 0.5) + 1 i2 = int((fre(j)*p)/df + 0.5) + 1 if(i2 .gt. nf1 ) then jmax = j - 1 ier = 22 goto 62 else smv(j,lw) = 0.0 smn(j,lw) = 0.0 sme(j,lw) = 0.0 do 13 i = i1, i2 smv(j,lw) = sv(i) + smv(j,lw) smn(j,lw) = sn(i) + smn(j,lw) sme(j,lw) = se(i) + sme(j,lw) 13 continue smv(j,lw) = smv(j,lw)/float(i2-i1 +1) smn(j,lw) = smn(j,lw)/float(i2-i1 +1) sme(j,lw) = sme(j,lw)/float(i2-i1 +1) endif endif 12 continue 62 continue endif if ( argout(3) .eq. 'tri' ) then ier = 0 read(argout(2),*) band p= 1. + band/100. c2 = alog(p) do 14 j = 1, nf bb = fre(j) * (p - 1.0/p) if ( bb .lt. 2*df ) then call interv1(j,fre,nf,df,sv,nf1,smv,n_win,lw) call interv1(j,fre,nf,df,sn,nf1,smn,n_win,lw) call interv1(j,fre,nf,df,se,nf1,sme,n_win,lw) else i1 = int( (fre(j)/p)/df ) + 2 i2 = int( (fre(j)*p)/df ) + 1 if(i2 .gt. nf1 ) then jmax = j - 1 ier = 22 goto 64 else c1 = alog(fre(j)/df) smv(j,lw) = 0.0 smn(j,lw) = 0.0 sme(j,lw) = 0.0 sweight = 0.0 do 15 i = i1, i2 weight = 1.0 - abs(c1 - alog(float(i-1)))/c2 smv(j,lw) = weight*sv(i) + smv(j,lw) smn(j,lw) = weight*sn(i) + smn(j,lw) sme(j,lw) = weight*se(i) + sme(j,lw) sweight = sweight + weight 15 continue smv(j,lw) = smv(j,lw)/sweight smn(j,lw) = smn(j,lw)/sweight sme(j,lw) = sme(j,lw)/sweight endif endif 14 continue 64 continue endif endif c if ( argout(1) .eq. 'konno-ohmachi' )then ier = 0 read(argout(2),*) bk if ( bk .lt. bkmin ) goto 99 if ( bk .gt. bkmax ) goto 99 c1 = 10.**(3.14159/bk) do 24 j = 1, nf bb = fre(j) * (c1 - 1.0/c1) if ( bb .lt. 2*df ) then call interv1(j,fre,nf,df,sv,nf1,smv,n_win,lw) call interv1(j,fre,nf,df,sn,nf1,smn,n_win,lw) call interv1(j,fre,nf,df,se,nf1,sme,n_win,lw) else i1 = int( (fre(j)/c1)/df ) + 2 i2 = int( (fre(j)*c1)/df ) + 1 if(i2 .gt. nf1 ) then jmax = j - 1 ier = 22 goto 74 else smv(j,lw) = 0.0 smn(j,lw) = 0.0 sme(j,lw) = 0.0 sweight = 0.0 c2 = alog10(df/fre(j)) do 25 i = i1, i2 c3 = bk * ( c2 + alog10(float(i-1))) if ( abs(c3) .lt. 0.001 ) then weight = 1.0 else weight = (sin(c3)/c3)**4 endif smv(j,lw) = weight*sv(i) + smv(j,lw) smn(j,lw) = weight*sn(i) + smn(j,lw) sme(j,lw) = weight*se(i) + sme(j,lw) sweight = sweight + weight 25 continue smv(j,lw) = smv(j,lw)/sweight smn(j,lw) = smn(j,lw)/sweight sme(j,lw) = sme(j,lw)/sweight endif endif 24 continue 74 continue endif c if ( ier .eq. 1 ) write(*,100)smooth if ( ier .eq. 22 )then if(jmax - jmin .lt. 2 )then ier = 1 write(*,400) else ier = 2 if (lw .eq. 1 ) write(*,300) endif endif c return c 99 write(*,200)smooth,bkmin,bkmax ier = 1 return c 97 write(*,*) ' This error should be impossible! : Fmax > Fny', # ' or Fmin < df (apply_smooth)' ier = 1 return c c 100 format(' ERROR : smoothing parameter not recognized: ',a80) 200 format(' ERROR : smoothing parameter: ',a80,/, # ' allowed Konno_Ohmachi bandwidth is between ',f5.1,' - ',f5.1) 300 format(' WARNING : output frequency limits have been changed', # ' by the smoothing procedure') 400 format(' ERROR : smoothing window too large:',/, # ' number of significant smoothed spectrum values ', # 'less than 2 !!') c c end