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,
xis the Fourier transform. To find the value of element
nin the data vector, you have to integrate. So
nis the index and
domegais 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.