function wers_thresh, data, thresh, freq_range, neg=neg

common path, home_path, single_path, avg_path, ps_path, latadj_path, dv_path, $
             pca_path, bhv_path, image_path, ers_path, ica_path
common params, n_ids, n_points, n_chans, period, epoch_range, filt_width, base_range
common wers_params, n_wers_pts, wers_range, wers_base_range, n_scales, scales

data2 = data
data2[n_ids:*,*,*] = 0.

chan_ix = [indgen(61),64]  ; monopolar channels, excluding Nas (not used) (note A1/A2 are mastoids)
n_chan_ix = n_elements(chan_ix)

;scales = dindgen(n_scales) * 0.125
;scales = 2d0^(scales)*(2*period)
scale_range = 1000./freq_range
scale_ix = where((scales GE scale_range[1]) AND (scales LE scale_range[0]), n_scale_ix)

pts_ix = n_ids + indgen(n_wers_pts) + ms_to_pts(wers_range[0],/wers)

for i=0,n_chan_ix-1 do $
  for j=0,n_scale_ix-1 do begin
    cycle = round(scales[scale_ix[j]]/period)
    ;cycle = round(0.5 * scales[scale_ix[j]]/period)
    ;cycle=1

    if (n_elements(thresh) EQ 1) then $
      if (keyword_set(neg)) then $
        for k=0,n_wers_pts-cycle do begin
          sel = data[pts_ix[k]:pts_ix[k]+cycle-1,scale_ix[j],chan_ix[i]] LE -(thresh)
          if (total(sel) EQ cycle) then $
            data2[pts_ix[k]:pts_ix[k]+cycle-1,scale_ix[j],chan_ix[i]] = data[pts_ix[k]:pts_ix[k]+cycle-1,scale_ix[j],chan_ix[i]]
        endfor $
      else $
        for k=0,n_wers_pts-cycle do begin
          sel = data[pts_ix[k]:pts_ix[k]+cycle-1,scale_ix[j],chan_ix[i]] GE thresh
          if (total(sel) EQ cycle) then $
            data2[pts_ix[k]:pts_ix[k]+cycle-1,scale_ix[j],chan_ix[i]] = data[pts_ix[k]:pts_ix[k]+cycle-1,scale_ix[j],chan_ix[i]]
        endfor $
    else $
      for k=0,n_wers_pts-cycle do begin
        sel_lo = data[pts_ix[k]:pts_ix[k]+cycle-1,scale_ix[j],chan_ix[i]] LE thresh[0]
        sel_hi = data[pts_ix[k]:pts_ix[k]+cycle-1,scale_ix[j],chan_ix[i]] GE thresh[1]
        if ((total(sel_lo) EQ cycle) OR (total(sel_hi) EQ cycle)) then $
          data2[pts_ix[k]:pts_ix[k]+cycle-1,scale_ix[j],chan_ix[i]] = data[pts_ix[k]:pts_ix[k]+cycle-1,scale_ix[j],chan_ix[i]]
      endfor
  endfor ;j

print, 'wers_thresh: done'
return, data2
end