mikeash.com pyblog/friday-qa-2012-10-26-fourier-transforms-and-ffts.html comments

Stephan Müller - 2014-04-08 13:42:21

Tue, 08 Apr 2014 13:42:21 GMT

[…] and sometimes with inverted components […]

float imag = buf[i] * sin(-angle);

Its a convolution… one of the signals needs to go "backwards".

Stephan Müller - 2014-04-08 13:33:34

Tue, 08 Apr 2014 13:33:34 GMT

What helped me understand the discrete fourier transform was looking at easy examples: A signal consisting of only two samples. Such a signal contains two frequencies, [1,1] (the DC offset) and [1,-1] (the Nyquist component). It is pretty obvious that by combining those two vectors, any two samples can be "created".

Larger signals require phase offsets (hence the complex numbers): [0,1,0,-1] (sampled sine) can be shifted to create [1,0,-1,0] (sampled cosine).

And the final missing link to understanding fourier (at least for me) was realising that different frequencies are orthogonal: if you multiply two different frequencies sample by sample and sum up the results, the sum will always be 0. Hence if you do this with a complex signal and a single frequency, all other frequencies in the signal sum up to 0. (try summing up [0,1,0,-1] and [1,-1,1,-1]).

Stephan Müller - 2014-04-08 13:02:26

Tue, 08 Apr 2014 13:02:26 GMT

[…] for reasons I don't fully understand, the results of this function are off by a factor of two […]

(Note that the whole article is not about the fourier transform, but the discrete fourier transform. But anyway…)

While the scaling of the transformed result is nothing but convention, it can still be helpful to understand where the 2 comes from. Lets take an easy signal, say sin(x) over a whole period (x in [0, 2pi]).

This signal contains exactly one frequency. So the transformation basically becomes summing up sin(angle)*sin(angle). Lets just assume the sampling rate is really high and integrate the function from 0 to 2pi (which gives us pi). Now I want the result to only depend on the signal amplitude and not the signal length, so I divide by the signal length of 2pi, which gives us our famous 1/2 (even though the original amplitude was 1).

In short: the factor of 2 is created by squaring the sine! The transformed signal has more direct meaning if it is multiplied by 2, but any scaling is just convention.

tgree - 2014-02-10 00:56:01

Mon, 10 Feb 2014 00:56:01 GMT

Your spinny description of the FFT was fantastic, and breaking it down into sin/cos instead of complex exponential arithmetic really helped to clarify what is going on for me.

Thanks for your article!

:) - 2013-10-12 20:42:40

Sat, 12 Oct 2013 20:42:40 GMT

*Of course* it'll be twice as large.

It's a sine wave; it has a positive and a negative peak. During the positive peak you're facing in the right direction, during the negative peak you're facing the opposite direction and walking backwards.

So you're walking in the same direction with a speed corresponding to the sampled intensity twice per cycle -- and that works out to a speed corresponding to one times the amplitude.

The inverted sign can be explained just as easily: you're in the *negative* halfwave when you think you are in the *positive* halfwave, and the other way around.

Jarryd - 2013-06-03 20:50:24

Mon, 03 Jun 2013 20:50:24 GMT

Thanks for the great explanation of the FFT.

Should all of these steps be done.

    vDSP_ctoz((COMPLEX *)data, 2, &out, 1, bufferFrames / 2);
    vDSP_fft_zrip(fftSetup, &out, 1, bufferlog2, FFT_FORWARD);
    vDSP_fft_zrip(fftSetup, &out, 1, bufferlog2, FFT_INVERSE);
    float outData[bufferFrames];
    vDSP_ztoc(&out, 1, (COMPLEX *)outData, 2, bufferFrames / 2);
    cblas_sscal(bufferFrames, 1.0 / (bufferFrames * 2), outData, 1);

I dont get where the dominant frequency can be found.

mikeash - 2012-11-09 21:40:26

Fri, 09 Nov 2012 21:40:26 GMT

SSteve: Yes, I can imagine. The fun thing about this bug is that, since OS X's audio pipeline is all floating-point, values in excess of 1.0 will, at least on some hardware, be louder than your set output volume. I had my speakers turned all the way down, but this stuff still came out at full blast.

SSteve - 2012-11-09 21:18:01

Fri, 09 Nov 2012 21:18:01 GMT

"This needs to be removed before treating the result as audio data, unless you like frightening everybody in the restaurant with the sudden blast of sound from your MacBook (ask me how I know this)."

I was doing beta testing for an audio software company during the PPC to Intel transition. The occasional endian bug could be deafening when we started playback.

jamie - 2012-11-01 22:36:23

Thu, 01 Nov 2012 22:36:23 GMT

By the way, your narrative description of the Fourier transform is the best I've read. I've been trying to come up with a slightly more memorable one involving two turtles on a beach arguing over how long the tide lasts :)

nevyn - 2012-10-28 08:55:05

Sun, 28 Oct 2012 08:55:05 GMT

Yay, thank you! Will have to experiment more using this knowledge...

TJ Usiyan - 2012-10-27 00:54:07

Sat, 27 Oct 2012 00:54:07 GMT

The reason that they are scaled and how is explained here

http://developer.apple.com/library/mac/documentation/Performance/Conceptual/vDSP_Programming_Guide/UsingFourierTransforms/UsingFourierTransforms.html#//apple_ref/doc/uid/TP40005147-CH202-15952

mikeash - 2012-10-26 20:37:08

Fri, 26 Oct 2012 20:37:08 GMT

vasi: Sorry, a bit of habitual shorthand there. rate refers to the sample rate of the audio, e.g. 44100. hz is the target frequency for the analysis.

jamie - 2012-10-26 19:54:47

Fri, 26 Oct 2012 19:54:47 GMT

"rate" is the sample rate of the float *.

"The result of this inverse transformation, for reasons I don't fully understand, comes out multiplied by a factor of bufferFrames * 2."

A bare FFT always gives amplitudes as a product of the block size. You might be getting the 2 from the fact that the output is the real part of an analytic signal? I dunno the conventions for veclib, my FFT experience is completely from Pure Data (which is great for learning but a bit tough for making polished products).

Also while you could definitely use something like this for cool visualization, you're just doing a straight rectangular window here, so you're going to end up with a lot of spectral leakage in the forward transform and when you resynthesize it'll sound distorted.

vasi - 2012-10-26 16:33:59

Fri, 26 Oct 2012 16:33:59 GMT

The signature for Rotate() is really confusing:

COMPLEX Rotate(float *buf, int samples, float hz, float rate)

It's not clear to me what the last two parameters mean: 'hz' is presumably a rate, and 'rate' is presumably measured in Hertz. I guess one of them indicates the frequency of our rotation, but which one, and what does the other one mean?

Clark - 2012-10-26 15:31:56

Fri, 26 Oct 2012 15:31:56 GMT

Next up spherical bessel functions for fun and profit…

More seriously, do you know what the two other main classes of functions for analysis aren't in vDSP? Admittedly they aren't used quite as much but still. Spheres and cylinders are very useful for anlayzing symmetries and the like.