DTFT

Jason Nyberg suggested a Fourier transform. I did a little programming and came up with this:

(define (dtft samples)
  (lambda (omega)
    (sum 0 (vector-length samples)
        (lambda (n)
          (* (vector-ref samples n)
             (make-polar 1.0 (* omega n)))))))

Given a vector of uniformly spaced samples (over time), this returns the Discrete-time Fourier Transform (DTFT) of the sampled function. I vaguely remember doing something with these in my signals and systems course, but it was a long time ago and I've forgotten almost all of it.

The basic idea is fairly simple. Suppose the original data contains a regular signal at some frequency. I an arrow that rotates uniformly at some frequency and multiply it by my samples. If the arrow is rotating at the same frequency as the embedded signal, then it will be pointing in the same direction every time I hit one of the bumps caused by the embedded signal. These will all add up to a really big arrow.

The omega variable is the amount (in radians) that I rotate the arrow for each sample. I represent my arrow by a complex number 1 unit in length. As I walk through the samples, I rotate the arrow (counter-clockwise) by the appropriate angle, multiply it by the sample value, and add that to the running total. I have 1000 samples per second, so I'll need to do unit conversion if I want to keep on the same time scale (and this is confusing enough as it is).

(define (hertz->radians-per-sample samples-per-second)
  (lambda (hz)
    (/ (* *pi* 2 hz) samples-per-second)))

(define (dtft samples-per-second samples)
  (let ((convert (hertz->radians-per-sample samples-per-second)))
    (lambda (hz)
      (let ((omega (convert hz)))
        (sum 0 (vector-length samples)
             (lambda (n)
               (* (vector-ref samples n)
                  (make-polar 1 (* omega n)))))))))

(define Y (dtft 1000 *data*))

(y 10)
;Value: -150.06753930463557+44.941130809509566i

(y 15)
;Value: -213.97389237338032-84.17677317997187i

So far so good, but I'm getting back complex numbers as answers. These encode both the direction and the length of my arrows. Right now, the direction isn't important to me. It tells me the phase of that frequency component, but I just want to know how big it is (the magnitude).

(magnitude (y 15))
;Value: 229.93598187410205

(magnitude (y 20))
;Value: 292.10647664795385

So let's plot the magnitude over a range of frequencies: The hair on my graph is at exactly 64Hz.

I have a couple of comments. First, you might wonder why I'm using complex numbers. I could have done the arrows by using a two-component vector or even by keeping track of X and Y variables. The answer is that complex numbers have the concept of rotation ‘built in’. If I were explicitly modeling the arrows, I'd have a lot of trigonometric functions kicking around. Since I'm modeling frequencies and such, I should have sines and cosines all over the place. These are all tucked in to the complex numbers. There is a 2π visible in the conversion routine, but all the trig is gone.

Second, (and this confused me a bit), my original graph is discrete, but the DTFT is continuous. I can get the magnitude for any frequency, even if it is not a multiple of the sampling frequency. (This confused me because at first I only looked at multiples of the sample frequency and I didn't find the peak. I assumed I got the program wrong and did a lot of

(magnitude (y (+ 40 3/5)))
;Value: 222.50621191745878

In fact, 64 Hz is 15.625 ms, which is not an even multiple of my 1 ms sampling rate.

So what is special about 64 Hz? It seems that this is the rate at which Windows does its time slicing.

Now that I know what the hair is, the next problem is removing it.