Going backwards

First a “Thank you!” to all that have pointed me towards notch filters. I expect that I'll end up using some sort of notch filter.

Part of the reason I'm doing this the hard way is to understand what I really want to do. I know the trail has already been blazed by Fourier, Laplace, Butterworth and Chebyshev, but I want to find out why this direction is the easiest and best way to get where I want to be. So back to the code.

Having found that enormous spike at 64 Hz, I want to separate my data into two components: the 64 Hz component and the remainder. I had a clever idea. What if I simply trim that spike out of the Fourier transform and then do an inverse Fourier transform to get the original data, but without the spike. There are a couple of problems with this. The first is that the inverse Fourier transform is harder to compute. Here is the simple code:

(define (inverse-dtft x)
  (lambda (n domega)
    (define (f omega)
      (* (x omega)
         (make-polar 1 (* omega (- n)))))

    (/ (integrate f (- *pi*) *pi* domega)
       (* *pi* 2))))

In this code, x is the Fourier transform. To find the value of element n in the data vector, you have to integrate. So n is the index and domega is the step size for the integration. It turns out that you need really small steps. The integral does not converge quickly.

It's a bit harder to understand why this works. The Fourier transform was easy to explain with the rotating arrows and such. This integral has a rotating arrow, but it rotates in the opposite direction. (I got it wrong at first. It doesn't work that way.) Think of the Fourier transform as a water sprinkler spraying your data across the lawn. The rotational speed of the sprinkler determines what data goes where, and if it matches a frequency component in your data, you get an uneven watering pattern. Now the inverse Fourier transform is you with a bucket running like mad trying to catch the water. That's a bit harder to do.

So with the integral converging so slowly, I did a few tricks. First, I memoized the computation of the transform. This takes a pile of memory, but the speedup is worth it. Second, I used Romberg integration. Nonetheless, I found that it still takes quite a while to invert the transform.

I mentioned a couple of problems. Besides it being slow, it doesn't work. After trimming out the spike at 64 Hz I found that the curve still has its hair. Apparently the harmonics carry a fair amount of weight. I tried trimming a couple. Very little effect on the hair, but I was starting to munge up the other data. I tried a different trick.

(define (y1 x)
  ;; Squash all frequencies above .1 (radians per second) 
  (if (> x .1)
      0
      (y x)))

Since the hair is at about .4, I figured this ought to do it.

It doesn't. I have to admit that I'm a bit surprised at this. I figured it would be trivial to trim the hair if I knew what frequency it was at and simply clipped it out. Sure, a notch filter would be far more efficient to compute, but isn't the concept right? I guess it isn't.

Well, then, back to basics.