func remove_emi(x, win)
 ;win = -1 is a flag that this is an internal call without a window
 if win eq -1 then {
 nx = dimen(x,1)
 ny = dimen(x,0)
 } else {
 nx = dimen(x,0)
 ny = dimen(x,1)
 }
 sc2d,x,a,b,c,d
 ;we want to examine power spectrum so don't need a displayable version,
 ;also must keep the transforms for inverse
 xq=(a*a+b*b+c*c+d*d)
 aq=2.*(c*d-a*b)
 dq=xq+aq
 aq=xq-aq
 nf = dimen(aq,1)*2 -1
 ;not concerned about the edges
 ;look for peaks along the 12.5 and 42 KHz lines
 
 i12k = fix(26*ny/1024)
 i42k = fix(86*ny/1024)
 ty,'power lines =', i12k, i42k
 in12 = peakseek(aq,dq,i12k,win,0.02,8,35*nx/1024,0.05,5)
 inall = in12
 ;also try closer to DC but with higher limit
 in12 = peakseek(aq,dq,i12k,win,0.05,8,10*nx/1024,0.2,5)
 if isarray(inall) then {
   if isarray(in12) then inall = [inall, in12]
   } else inall = in12
 ;and try at 27
 i12k = fix(27*ny/1024)
 ty,'power lines =', i12k 
 in12 = peakseek(aq,dq,i12k,win,0.05,8,10*nx/1024,0.2,10)
 if isarray(inall) then {
   if isarray(in12) then inall = [inall, in12]
   } else inall = in12
 in42 = peakseek(aq,dq,i42k,win,0.01,8,10*nx/1024,0.1,5)
 if isarray(inall) then {
   if isarray(in42) then inall = [inall, in42]
   } else inall = in42
 ;85 sometimes
 i42k = fix(85*ny/1024)
 ty,'power lines =', i42k 
 in42 = peakseek(aq,dq,i42k,win,0.01,8,10*nx/1024,0.1,5)
 if isarray(inall) then {
   if isarray(in42) then inall = [inall, in42]
   } else inall = in42
 ;also try the harmonic of 12.5 KHz
 ty,'power lines =', i12k*2 
 in25 = peakseek(aq,dq,i12k*2,win,0.01,8,35*nx/1024,0.1,5)
 if isarray(inall) then {
   if isarray(in25) then inall = [inall, in25]
   } else inall = in25

 ;optionally show the plots and a map of the peak removed
 if win ge 0 and defined($viewemidebug) then {
 if $viewemidebug(win) ge 1 then {
 ;show these as a bit array that matches the top half of the 2-D power spectrum
 nf2 = dimen(a,0)
 p2d = zero(bytarr(nf, nf2))
 ;the inall indices are already oriented for top half of p2d
 p2d(inall) = 1
 p2d=concat(reverse(p2d),p2d(*,1:*))

 $data4 = p2d
 view, 4
 }}
 
 ;now zero these bad guys
 
 if isscalar(inall) then {
   ;nothing to do
   ty,'no EMI discovered'
   return, x  }
 
 nx = dimen(a,0)	ny = dimen(a,1)
 ;note redef of nx and ny above

 if win eq -1 then {
 ixp = inall%nf - nx + 1
 iyp = inall/nf			;this is always positive
 iq = ixp lt 0
 } else {
 iyp = inall%nf - ny + 1
 ixp = inall/nf			;this is always positive
 iq = iyp lt 0
 }
 iyp = abs(iyp)
 ixp = abs(ixp)
 ind = ixp + nx*iyp
 inp = abs(sieve(ind, iq))
 if isarray(inp) then {
 ;this is the "lower rhs" zone
 xq = .5*( a(inp) - b(inp))
 a(inp) = xq
 b(inp) = - xq
 xq = .5*( c(inp) + d(inp))
 c(inp) = xq
 d(inp) = xq
 }

 ;then the other half
 iq = not(iq)
 inp = abs(sieve(ind, iq))
 if isarray(inp) then {
 ;this is the "upper rhs" zone
 xq = .5*( a(inp) + b(inp))
 a(inp) = xq
 b(inp) = xq
 xq = .5*( c(inp) - d(inp))
 c(inp) = xq
 d(inp) = -xq
 }


 ;finally, need to check for edges, note how x and y are switched
 iq = (ixp eq 0) or (ixp eq (ny-1)) or (iyp eq 0) or (iyp eq (nx-1))
 inp = abs(sieve(ind, iq))
 if isarray(inp) then {
 a(inp) = 0.0
 c(inp) = 0.0
 d(inp) = 0.0
 } 

 sc2db, y, a, b, c, d
 
 if win ge 0 and defined($viewemidebug) then {
   if $viewemidebug(win) ge 1 then { $data8 = x-y  view, 8 }
 }

 return, y
 
 endfunc
 ;=============================================================================
 ;=============================================================================