Wednesday, July 1, 2009

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.

2 comments:

jpc said...

DFT is it's own inverse (the scalling may be different depending on what scalling you chosen for forward DFT). [1]

Taking inverse DFT of modified spectrum is problematic because of the implied periodicity of the signal (DFT computes as if the signal was wrapped around from end to the beginning). This can be worked around (look for explanations of the overlap-add method).

For a very good introduction (or refreshment) to DSP try the book that can be downloaded at [2]. I have a deadtree copy on my bookshelf and I cannot overstate how great it is.

[1]: http://en.wikipedia.org/wiki/Discrete_Fourier_transform#Definition

[2]: http://dspguide.com/

jpc said...

PS. As it was mentioned earlier you should try and experiment with median filters. They are non-linear (so they don't fit the Fourier aka frequency domain theory) but can be really useful for eliminating impulse noise without blurring the signal.