function bgndtest,tavg,avg,dsrc,yloc,dmin,dmax

   bl   = [1,2,12,14,15,16,24,25,26]
   chan = [0,1, 5, 7, 9,11,13,15,17]
   indx = (where(bl eq dsrc,ndx))(0)
   if (ndx ne 1) then begin
      print,'Invalid datasource specified.'
      return,-1
   endif

   device,decompos=0
   loadct,39
   xsz = n_elements(tavg)
   if (bl(indx) gt 10) then begin
      datarr = congrid(rotate(reform(sqrt(avg[chan[indx],*,*]^2$
                                      +avg[chan[indx]+1,*,*]^2)),-5),xsz,100)
   endif else begin
      datarr = congrid(rotate(reform(avg[chan[indx],*,*]),-5),xsz,100)
   endelse
   if (n_elements(dmin) eq 0) then dmin = 0
   if (n_elements(dmax) eq 0) then dmax = max(datarr)
   tvscl,datarr<dmax>dmin,0,yloc

   ans = ''
   ; Ask if we want to continue with this scaling.
   read,'Continue?: [Y] ',ans
   ans = strupcase(ans)
   if (ans eq '') then ans ='Y'
   if (ans ne 'Y') then return,-1

   ; Set background points with mouse
   x = 0
   iarr = intarr(100)
   i = 0
   for i = 0, 99 do begin
      cursor,x,y,/device,4
      if (x ge xsz) then goto,done
      plots,/device,[x,x],[0,99]+yloc
      iarr[i] = nint(x)
   endfor

done:
   iarr = iarr[0:i-1]
   iarr = iarr[sort(iarr)]

   nf = n_elements(avg[0,*,0])
   npts = n_elements(avg[0,0,*])

   if (bl(indx) lt 10) then begin
      ; Do Total Power
      dat = (datout = reform(avg[chan[indx],*,*]))
      ; Perform 2-point background subtraction
      for i = 1, n_elements(iarr)-1 do begin
         ; k1 and k2 are the range of indexes in AVG to apply the subtraction
         ; for this pair of background settings
         k1 = iarr[i-1]
         k2 = iarr[i]
         if (i eq 1) then k1 = 0
         if (i eq n_elements(iarr)-1) then k2 = npts-1
         ; Accum background over 3 points around each sub point
         spec1 = (spec2 = fltarr(nf))
         for j = 0, nf-1 do begin
            spec1[j] = (moment(dat[j,iarr[i-1]-1:iarr[i-1]+1],/nan))(0)
            spec2[j] = (moment(dat[j,iarr[i]-1:iarr[i]+1],/nan))(0)
         endfor
         ; Perform subtraction
         for k = k1,k2 do begin
            bgnd = spec1 + (spec2 - spec1)*(k-iarr[i-1])/(iarr[i]-iarr[i-1])
            datout[*,k] = dat[*,k] - bgnd
         endfor
      endfor
   endif else begin
      ; Do correlated data
      y = (yout = reform(avg[chan[indx],*,*]))
      x = (xout = reform(avg[chan[indx]+1,*,*]))
      ; Perform 2-point background subtraction
      for i = 1, n_elements(iarr)-1 do begin
         ; k1 and k2 are the range of indexes in AVG to apply the subtraction
         ; for this pair of background settings
         k1 = iarr[i-1]
         k2 = iarr[i]
         if (i eq 1) then k1 = 0
         if (i eq n_elements(iarr)-1) then k2 = npts-1
         ; Accum background over 3 points around each sub point
         y1 = (y2 = fltarr(nf))
         x1 = (x2 = fltarr(nf))
         for j = 0, nf-1 do begin
            vals = y[j,iarr[i-1]-1:iarr[i-1]+1]
            ifinite = where(finite(vals),nfinite)
            if (nfinite eq 1) then y1[j] = vals[ifinite] $
                              else y1[j] = (moment(vals,/nan))(0)
            vals = y[j,iarr[i]-1:iarr[i]+1]
            ifinite = where(finite(vals),nfinite)
            if (nfinite eq 1) then y2[j] = vals[ifinite] $
                              else y2[j] = (moment(vals,/nan))(0)
            vals = x[j,iarr[i-1]-1:iarr[i-1]+1]
            ifinite = where(finite(vals),nfinite)
            if (nfinite eq 1) then x1[j] = vals[ifinite] $
                              else x1[j] = (moment(vals,/nan))(0)
            vals = x[j,iarr[i]-1:iarr[i]+1]
            ifinite = where(finite(vals),nfinite)
            if (nfinite eq 1) then x2[j] = vals[ifinite] $
                              else x2[j] = (moment(vals,/nan))(0)
         endfor
         ; Perform subtraction
         for k = k1,k2 do begin
            ybgnd = y1 + (y2 - y1)*(k-iarr[i-1])/(iarr[i]-iarr[i-1])
            yout[*,k] = y[*,k] - ybgnd
            xbgnd = x1 + (x2 - x1)*(k-iarr[i-1])/(iarr[i]-iarr[i-1])
            xout[*,k] = x[*,k] - xbgnd
         endfor
      endfor
      datout = sqrt(yout^2+xout^2)
   endelse
   dat = congrid(rotate(reform(datout),-5),xsz,100)
   tvscl,dat>dmin,0,yloc
return,dat
end