History and Derivation of the Fast Fourier Transform

The Claerbout book also gives a 27-line Fortran program to perform an FFT.

Especially in math, a distinction is often drawn between

    FT   f(real) -> F(real)
    DFS  f(real) -> F(int)
    DTFT f(int)  -> F(real)
    DFT  f(int)  -> F(int)

Only one of these, the DFT, is easily handled by a computer, but they all have an important role in both theory and practice. For example, why does JPEG use a Discrete Cosine Transform rather than a Discrete Fourier Transform? What are the pitfalls of approximating a continuous domain with discrete samples? Why are engineers obsessed with window functions? A continuous domain will lead you to the right answer. You can't easily understand the DFT (or FFT) without the FT and DFT. DFS has its own uses.

To many (applied) mathematicians, the fact that a computer uses the FFT algorithm to compute a DFT is largely immaterial and not particularly worth mentioning, because they are more interested in the questions in the preceding paragraph than questions like "how many points can I take on X hardware while maintaining Y performance?"

My impression is Winograd FFT is still O(N log N), i.e. exactly the same complexity as Cooley-Tukey. That said, they work best on different N, which is possibly what you're referring to?

You are correct on the sine and cosine thing, the minus is switched. Will fix, thanks.

it will depend on the implementation and processor if it's faster to calculate simultaneous sin and cos of a given angle, or sin(w) and sin(w+pi/2), but you are right you could do it all with sin or cos only.

But there is already a different thing called a discrete cosine transform :)

One reason is that in floating-point processing, in some cases adding a constant to an argument can truncate and lose the original argument, so using a single function will fail in a class of pathological edge cases:

    import numpy as np
    a = np.pi/2
    print(a)
    > 1.5707963267948966
    print(a == a + 1e-15)
    > False
    print(a == a + 1e-16)
    > True
    print(np.sin(1e-16))
    > 1e-16
    print(np.cos(np.pi/2+1e-16))
    > 6.12323399574e-17
    print(np.sin(1e-16) == np.cos(np.pi/2+1e-16))
    > False

I believe the reason is due to cosine being even, i.e. cos(-x) = cos(x), and sine being odd, sin(-x) = -sin(x). The sum of even functions is still even, and the sum of odd functions is still odd. If you have a function that is symmetric about the origin, say f(x) = x^2, then you can approximate it with only cosines, but in general functions are not odd or even so you need to use a mix of sines and cosines.

Well, you may not like this answer, but:

If cosine is just sine rotated by 90°, then we really just use sines. Except we write cosine instead of a sine with a 90° rotation because it looks less confusing :)

The real goal is to arrive at the exponential representation, where the sines and cosines are absorbed into the exponential functions, which is later needed to write down the fourier integrals.

Your question boils down to preference. The two things you name are, in fact, equal. Math people just choose to write cos(x) instead of sin(x+pi/2)

The "true" Fourier transform is to write the function as an infinite sum of e^(2 pi i n x) with complex coefficients and n an integer. (This is the "true" Fourier transform because of a connection between periodic functions and circles, and because the representation theory of the circle group says that these are the so-called irreducible representations. This is extremely abridged -- the point is just to mention relevant words.)

Due to Euler's formula, there is a correspondence between a linear combination of e^(2 pi i n x) and e^(-2 pi i n x) and a linear combination of cos(2 pi n x) and sin(2 pi n x). In other words, both pairs of functions span the same subspace (with complex coefficients).

We need all of the exponentials to represent a periodic function, with positive and negative n, and so we need both cos and sin. If you leave out cos, you can only represent odd functions, and leaving out sin, even functions.

But: if your question is just why not just use sin(x) and sin(x-pi/2) instead of sin(x) and cos(x), well, it's the same thing. Just because you write cos(x) like sin(x-pi/2) doesn't mean it's not still cos(x)!

(Fun fact: cosine means "complementary sine," the sine of the complementary angle in the triangle.)

Because with a Fourier series/transform you only get to play with the amplitude of the sine and cosine functions, not the offset. If you only had cosine functions in your sum you'd only be able to represent symmetric functions (and antisymmetric if only sine). Except actually I lied; if you sum both sine and cosine functions you do get to specify offsets. Think of this trig identity: sin(x + offset) = [cos(offset)] sin(x) + [sin(offset)] cos(x). Think about the terms in brackets as just being amplitudes. By having both sine and cosine in the sum you can represent any arbitrary offset by just changing their relative amplitude. This is one of those things which is nicer when you think about it in terms of complex numbers.

The series \sum a_n\cos(n\theta) + \sum b_n\sin(n\theta) and \sum a_n\sin(n\theta + \pi/2) + \sum b_n\sin(n\theta) are literally, exactly the same; it is purely a matter of convention which one you choose to write. In more programming-type language, 'cosine' is just syntactic sugar (or 'sine' is, if you prefer; or, I would rather say, both are, for the real and complex parts of the exponential); but sometimes syntactic sugar is handy!

I am genuinely interested in the source for that if you have it.