pro c_spec, xx,yy,tim_secs,coh,phase, thres=thres if n_elements(thres) eq 1 then thres = 0.9 mask_x = xx mask_y = yy mask_x(*) = 1 mask_y(*) = 1 miss = where(xx ge 998.) if miss(0) ne -1 then mask_x(miss) = 0 miss = where(yy ge 998.) if miss(0) ne -1 then mask_y(miss) = 0 ntim = n_elements(xx(*,0)) nh = n_elements(xx(0,*)) rat_x = total(mask_x,1)/ntim rat_y = total(mask_y,1)/ntim fx = complexarr(ntim/2+1,nh) fy = complexarr(ntim/2+1,nh) ;print,ntim,nh,rat_x,fx xx1=fft(xx) yy1=fft(yy) for i=0,nh/2 do begin append_array,xx2,xx1(i) append_array,yy2,yy1(i) endfor for i=0,nh-1 do begin fx(*,i) = ffft2(xx(*,i),tim_secs) fy(*,i) = ffft2(yy(*,i),tim_secs) endfor x_p = abs(fx)^2; *per y_p = abs(fy)^2; *per nomiss = where(rat_x ge thres and rat_y ge thres,nnomiss) if nomiss(0) ne -1 then begin x_s = total(x_p(*,nomiss),1) /nnomiss y_s = total(y_p(*,nomiss),1) /nnomiss endif xy_r = conj(fx)*fy ;xy1=complexarr(n_elements(xy_r)) ;xy1(0)=0.5*xy_r(0)+0.5*xy_r(1) ;for i=1,n_elements(xy_r)-2 do begin ; xy1(i)=0.25*xy_r(i-1)+0.5*xy_r(i)+0.25*xy_r(i+1) ;endfor ;xy1(n_elements(xy_r)-1)=0.5*xy_r(n_elements(xy_r)-2)+0.5*xy_r(n_elements(xy_r)-1) xy1=smooth(xy_r,3) nomiss = where(rat_x ge thres and rat_y ge thres,nnomiss) if nomiss(0) ne -1 then begin xy_e = total(xy_r(*,nomiss),1) /nnomiss endif xy_s = abs(xy1) kxy = float(xy1) qxy = -imaginary(xy1) x_p=smooth(x_p,3) y_p=smooth(y_p,3) coh = xy_s/sqrt(x_p*y_p) phase = atan(qxy,kxy) end