How to generalise the Fourier transform The Fourier transform approximates a signal using a bunch of sine and cosine waves. The inverse Fourier transform then reconstructs the original signal from this information.
I am told that it's possible to decompose a signal using some other set of functions, rather than the usual sine and cosine. My question is, how do you do this? 
For a start, I'm assuming that for your set of functions to be able to approximate every possible signal, you need to have "enough" of these functions, and ideally you want them to be "different" such that each one measures an unrelated aspect of the signal.
To be completely clear, I'm mostly interested in the case of digital sampled data. But the continuous case might be interesting too...

Edit:
I'm not sure why my question isn't producing any answers. Maybe it's because nobody actually knows the answer, or maybe it's because the answer is "too obvious" to somebody who actually possesses formal mathematical training. I'm not sure. But I've wanted to know the answer to this question for years, so let's try one more time...
The discrete Fourier transform works by computing the correlation of the input signal with several different sine waves. The inverse transform then adds together the specified amplitudes of those waves, recovering the original signal. That much seems clear.
It looks like I could just invent a family of functions to use instead of the sine and cosine functions, and do exactly the same process... except that when I do this, it doesn't work in any way, shape or form. If I transform and then inverse-transform, I get gibberish. And I don't know why... but it seems like the key phrase is "complete set of orthonormal functions", whatever that means.

Update:
I had assumed that if I could just find a system of basis functions such that none of them are correlated, and their number equals the number of points in the input, the transform would work. Apparently, it does not.
Consider the following set of functions:
$$f_1 = [1,1,1,1]$$
$$f_2 = [0,1,0,-1]$$
$$f_3 = [1,0,-1,0]$$
$$f_4 = [1,-1,1,-1]$$
Clearly, there are 4 functions. As far as I can tell, none of them are correlated. For example,
$$f_1 * f_4 = (1 * 1) + (1 * -1) + (1 * 1) + (1 * -1) = 1 - 1 + 1 - 1 = 0$$
If we take, say, $x = [1,2,3,4]$ and compute the correlations, we get
$$f_1 * x = 1 + 2 + 3 + 4 = 10$$
$$f_2 * x = 0 + 2 + 0 - 4 = -2$$
$$f_3 * x = 1 + 0 - 3 + 0 = -2$$
$$f_4 * x = 1 - 2 + 3 - 4 = -2$$
Now, computing $10 f_1 - 2 f_2 - 2 f_3 - 2 f_4$, we get
$$10 + 0 - 2 - 2 = 6$$
$$10 - 2 + 0 + 2 = 10$$
$$10 + 0 + 2 - 2 = 10$$
$$10 + 2 + 0 + 2 = 14$$
Clearly $[6, 10, 10, 14]$ is nothing like $[1,2,3,4]$, even with scaling. So... what am I missing?
 A: Your question is at the core of not only signal processing but differential equations and orthogonal special functions, fields of study that have a long history and are still active and evolving, so it's a daunting task to point out where you could start your studies.
The Wiki leonbloy pointed out, Generalized Fourier Series, and also the Wiki Green's Function with the section on eigenvalue expansions introduce the jargon that you should be thoroughly familiar with.
The basic algorithm is to find dual sets of eigenvectors/eigenfunctions parametrized by a continuous (e.g., $\omega$ below) or discrete index (e.g., $n$ below), that satisfy completeness and orthogonality relations encapsulated in Dirac delta function resolutions such as that for the Fourier transform
$$\delta(x-y)= \int_{-\infty}^{\infty}\exp(i2\pi \omega x)\exp(i2\pi \omega y)d\omega$$ 
giving
$$\int_{-\infty}^{\infty}f(y)\delta(x-y)dy=f(x)=\int_{-\infty}^{\infty}\exp(i2\pi \omega x)\int_{-\infty}^{\infty}f(y)\exp(i2\pi \omega y) dy d\omega$$
$$=\int_{-\infty}^{\infty}\exp(i2\pi \omega x)\hat{f}(\omega) d\omega,$$
or that for the eigenvectors of Sturm-Liouville differential operators over finite domains
$$\delta(x-y)=\sum_{n=0}^{\infty }\Psi_n(x)\Psi_n^*(y)$$
giving
$$f(x)=\sum_{n=0}^{\infty }\Psi_n(x)\int_{-\infty}^{\infty}f(y)\Psi_n^*(y) dy,$$ 
or Kronecker delta resolutions such as that for the associated Laguerre functions 
$$\frac{(n+\alpha)!}{n!}\delta_{mn}=\int_{0}^{\infty}x^{\alpha}e^{-x}L_{n}^{\alpha}(x)L_{m}^{\alpha}(x)dx$$ 
giving
$$f(x)=\sum_{n=0}^{\infty }\frac{n!L_{n}^{\alpha}(x)}{(n+\alpha)!}\hat{f}_n$$
with
$$\hat{f}_n=\int_{0}^{\infty}x^{\alpha}e^{-x}L_{n}^{\alpha}(x)f(x)dx.$$
The Fourier Transform and Its Applications by R. Bracewell is a really good book for grasping the fundamentals of the FT and DFT, as well as G. Strang's Introduction to Applied Mathematics.
Methods of Applied Mathematics by F. Hildebrand and Principles and Techniques of Applied Mathematics by B. Friedman give good intros to Fredholm theory and Green's functions.
More advanced books on harmonic analysis, such as J. Partington's Interpolation, Identification, and Sampling might be the next leap if you are comfortable with complex analysis (e.g., fractional linear transformations) and other integral transforms such as the Laplace transform. 
A: Since I saw your last answer to yourself and that I do not entirely agree, I will add my little stone to the edifice, hoping it might bring some clarity. It relies heavily on Mallat's book (a Wavelet Tour of Signal Processing, see e.g. chap 5) but I think most of it was already present in other's comments.
Let $f\in \mathcal H$ with a Hilbert space $\mathcal H$. Let $\{\phi_n\}_{n\in \Gamma}$ a family of vectors in $\mathcal H$ (anything really, this is where I somehow disagree with your point 1. and 2.) then you can define an operator based on the family of vectors as follows:
$$ \forall n\in \Gamma, \quad Uf[n] = \langle f,\phi_n\rangle $$
So basically this operator does just one thing, it associates to a function $f$ a set of numbers, the $n$th number corresponding to the inner product (or correlation as you put it) of $f$ with $\phi_n$.
Now the big thing is under which general conditions can one recover $f$ from its $Uf[n]$? Which is nothing but a question about the invertibility of $U$. In a practical setting, a typical example of a sequence is a family of sines and cosines up to a certain frequency. Already here, and rather intuitively, you can guess that it is not possible to exactly recover a general $f$ from its $Uf[n]$ but there are some $f$ for which it will be possible etc..
So from this, it should be clear that the characterization of the operator $U$ given the family of $\phi_n$ is essential to see whether or not we'll be able to recover the function (or signal) from its decomposition in the sequence of $\phi_n$.
Some more stuff from Mallat's book: the sequence is a frame of $\mathcal H$ if there exists two positive constants $A,B$ s.t for any $f\in\mathcal H$ one has:
$$ A\|f\|^2 \le \sum_{n\in\Gamma} |\langle f,\phi_n\rangle|^2 \le B\|f\|^2$$
and when $A=B$, the frame is said to be tight.
If you have that, then $U$ is called (very originally) a frame operator and one can prove that iff $U$ is a frame operator, it is invertible on its image with bounded inverse.
Some comments: 


*

*you can have redundancy in the family of $\phi_n$ (sometimes interesting, see e.g. curvelets)

*you can normalize your vectors to have $\|\phi_n\|=1$

*if the $\phi_n$ are linearly independent and form a frame then $A\le 1\le B$. 

*the frame is an orthonormal basis iff $A=B=1$. (e.g. Fourier basis wrt L2)


So basically, you can decompose a function/signal with respect to any kind of family but you might not be able to recover it (fully) from its decomposition. To give some insight about possible applications here are two (and there are many many more):


*

*Compression: express a signal in a family with a limited number of $\phi_n$ and remove those $Uf[n]$ which are ''too small''. The signal can then be stored with some loss but with (hopefully) very few ''important'' coefficients (hence the compression). For images this can be done very efficiently in a wavelet basis for example (e.g. JPEG2000). A good basis will be a basis in which the decomposition of $f$ has very few ''important coefficients'' $Uf[n]$ and the rest can be ignored. Fourier basis is usually pretty bad in that sense (tends to smear out the data). 

*Denoising: given a noisy signal expressed in a family, can try to recover a less noisy signal by inverting the frame operator on coefficients which are sufficiently big and hence should have a high signal to noise ratio.

A: As the previous answers have stated, your functions need to be an orthonormal basis for your procedure to work.  Your basis is orthogonal, but not normalized.  Try the same thing using 
$ \frac{1}{\sqrt{4}} [1, 1, 1, 1,]$
$ \frac{1}{\sqrt{2}} [0, 1, 0, -1]$
$ \frac{1}{\sqrt{2}} [1, 0, -1, 0]$
$\frac{1}{\sqrt{4}} [1, -1, 1, -1]$
You will have much more luck finding useful information on this by looking up dot products rather than correlation.  
A: You can generalize the Fourier transform by using a family of complete orthogonal functions. For instance if I were interested in a series expansion for functions on the real line $(-\infty,\infty)$ I could use the Hermite polynomial basis. 
The polynomials are defined by,
$$
H_n = (-1)^ne^{x^2}\frac{d^n}{dx^n}e^{-x^2}.
$$ 
These polynomials have two important properties.


*

*Orthogonality,
$$ \int_{-\infty}^{\infty} H_n(x) H_m(x) e^{-x^2}dx  = \begin{cases} \sqrt{\pi}2^n n! \text{ if } n=m \\ 0 \text{ otherwise} \end{cases}$$

*Completeness,
$$ f(x) = \sum_{j=0}^\infty c_j H_j(x)e^{-x^2/2} \qquad \text{(most functions)}$$
If we want the coefficients in the expansion then we just exploit the orthogonality of the basis functions,
$$  c_n =\frac{1}{2^n  n! \sqrt{\pi}} \int_{-\infty}^{\infty} f(x) H_n(x) e^{-x^2/2} dx $$
There are other families of complete orthogonal functions which can do the same thing. There are some limitations to these types of series representations which are shared by fourier series. In the context of the series above you are guaranteed to be able to represent square integrable functions.
For more on orthogonal polynomials see the pdf: "Abromowitz and Stegun Handbook of Mathematical Functions" chapter 22.
A: One aspect of this type of decomposition (i.e. Taylor, Fourier) is that it is an approximation that is made better iteratively. I have always thought of it as $n$-to-$\infty$-dimensional coordinate transformation that is reduced to $n$ dimensions with the folding operation (which is addition in the case of Taylor/Fourier series).
You can rewrite a vector space $(a,b)\in\mathbb{R}^2$ (such as $a+b i\in\mathbb{C}$) in terms of another vector space $(u,v)\in\mathbb{R}^2$ according to the Cauchy-Riemann equations. Maybe this will help you find a way to express the proper conditions for your work (which for CR boils down to "everything must be perpendicular everywhere).
A: Based on the hundreds of comments and the experimental work I've done, it appears the answer to my original question should have looked likt this:
A [discrete] function can be decomposed over a set of basis functions provided that the set of basis functions has the following properties:


*

*The correlation of every function with itself equals exactly one.

*The correlation of every function with every other function equals exactly zero.

*The number of basis functions is at least as large as the number of data points.


Once you have this set, you can "transform" your function simply by computing the correlation between your function and each of the basis functions. To "untransform", simply multiply each basis function by the coefficient you found, and add the results back together.
Somebody said something about a "Gram-Schmidt process" - and looking this up on Wikipedia gives me a nice way of taking a system of functions and making an orthonormal system from it.
The bounty on this question is still open for a few more hours. If anybody can summarise this better than I have, or add any interesting additional information (e.g., the exact conditions under which the set of functions is "complete"), you can still earn yourself a tidy little rep bonus...
