Hasta ahora hemos visto cómo aproximar funciones continuas por funciones polinómicas minimizando el error cometido basado en la norma euclídea (caso de aproximación por mínimos cuadrados) o la norma infinito (caso de la aproximación minimax).

En el caso de aproximación por mínimos cuadrados, hemos visto cómo aproximar una nube de puntos \((x_k,y_k),\ k=1,\ldots,m\) (caso discreto) o cómo aproximar una función continua en un determinado intervalo \([a,b]\) (caso continuo).

Sin embargo, el uso de polinomios no siempre es la mejor opción para aproximar una función o una nube de puntos.

Los polinomios tienen una tendencia natura a oscilar mucho, sobre todo cuando aumentamos el grado de los mismos.

Esto hace que las cotas de los errores cometidos sea mucho más grande que los errores reales.

En esta sección, consideraremos otro tipo de funciones para aproximar una nube de puntos o una función continua: las llamadas funciones trigonométricas.

Nos centraremos para simplificar los cálculos y para un mayor entendimiento de los conceptos que vamos a ver al intervalo \([-\pi,\pi]\).

Dado un entero positivo \(n\), las funciones trigonométricas consideradas \(\{\psi_0(x),\psi_1(x),\ldots,\psi_n(x),\psi_{n+1}(x),\ldots,\psi_{2n-1}(x)\}\) tienen la expresión siguiente: \[ \begin{align*} \psi_{0}(x)= & \frac{1}{2},\\ \psi_{k}(x)=&\cos(kx),\quad k=1,\ldots,n\\ \psi_{n+1}(x)=&\sin(x),\\ \psi_{n+k}(x)=&\sin(kx),\quad k=1,\ldots,n-1. \end{align*} \]

El espacio de aproximación \(T_n\) considerado será el espacio vectorial cuya base de funciones es el conjunto anterior \(\{\psi_0(x),\psi_1(x),\ldots,\psi_n(x),\psi_{n+1}(x),\ldots,\psi_{2n-1}(x)\}\): \[ \begin{align*} T_n = \{S_n(x)\ |\ S_n(x)=& \frac{a_0}{2}+\sum_{j=1}^{n-1} (a_j\psi_j(x)+b_j\psi_{n+j}(x))+a_n \psi_n(x)\\ = & \frac{a_0}{2}+\sum_{j=1}^{n-1} (a_j\cos(jx)+b_j\sin(jx))+a_n\cos(nx)\}. \end{align*} \] El problema que queremos resolver es el siguiente:

Sea \(f\in {\cal}[-\pi,\pi]\) una función continua en el intervalo \([-\pi,\pi]\), queremos hallar la función \(S_n(x)\in T_n\) que pertenece al espacio de aproximación definido anteriormente tal que la función error siguiente sea mínima: \[ \begin{align*} & E(a_0,a_1,\ldots,a_n,b_1,b_2,\ldots,b_{n-1})= \int_{-\pi}^{\pi}(f(x)-S_n(x))^2\, dx\\ &\quad = \int_{-\pi}^{\pi}\left(f(x)-\left(\frac{a_0}{2}+\sum_{j=1}^{n-1} (a_j\cos(jx)+b_j\sin(jx))+a_n\cos(nx)\right)\right)^2\, dx. \end{align*} \] La función error depende de \(2n-1\) coeficientes \(a_k\), \(k=0,\ldots,n\) y \(b_k\), \(k=1,\ldots,n-1\) que tenemos que hallar.

El Teorema siguiente nos da la función \(S_n(x)\) que minimiza la función error:

Teorema: función que minimiza la función error.

Sea \(f\) una función definida e integrable en el intervalo \([-\pi,\pi]\). Entonces la función que minimiza la función error es la siguiente: \[ S_n(x)=\frac{a_0}{2}+\sum_{j=1}^{n-1} (a_j\cos(jx)+b_j\sin(jx))+a_n\cos(nx), \] con \[ \begin{align*} a_k = & \frac{1}{\pi}\int_{-\pi}^\pi f(x)\cos(kx)\,dx,\ k=0,\ldots,n\\ b_k = & \frac{1}{\pi}\int_{-\pi}^\pi f(x)\sin(kx)\,dx,\ k=1,\ldots,n-1. \end{align*} \]

Observación.

La razón de que en la expresión \(S_n(x)\) aparezca \(\frac{a_0}{2}\) es que la expresión de \(a_k\) queda compacta, para \(k=0,\ldots,n\). De hecho, \(\displaystyle a_0 =\frac{1}{\pi}\int_{-\pi}^\pi f(x)\,dx.\)

Observación.

La función \(S_n(x)\) cuando \(n\to\infty\) se denomina desarrollo de Fourier de la función \(f(x)\): \[ S(x)=\frac{a_0}{2}+\sum_{k=1}^\infty (a_k\cos(kx)+b_k\sin(kx)). \]

Antes de desarrollar las ecuaciones anteriores veamos el siguiente lema que nos será muy útil:

Lema.

Sean \(j,k\) valores enteros. Entonces se cumplen las igualdades siguientes: \[ \begin{align*} \int_{-\pi}^{\pi}\sin(jx)\, dx= & 0,\ j=0,1,\ldots,\quad \int_{-\pi}^\pi \cos(jx)\, dx= 0,\ j=1,2,\ldots,\\ \int_{-\pi}^\pi \cos(jx)\cos(kx)\, dx = & 0,\mbox{ si }j\neq k,\quad \int_{-\pi}^\pi \cos^2(kx)\, dx = \pi,\ k=1,2,\ldots,\\ \int_{-\pi}^\pi \sin(jx)\cos(kx)\, dx = & 0,\ j,k=1,2,\ldots,\\ \int_{-\pi}^\pi \sin(jx)\sin(kx)\, dx = & 0,\ \mbox{ si }j\neq k,\quad \int_{-\pi}^\pi \sin^2(kx)\, dx = \pi,\ k=1,2,\ldots \end{align*} \]

El término del error para la función trigonométrica \(S_n(x)\) que da el mínimo queda simplificado tal como dice el resultado siguiente:

Proposición. Expresión del término del error.

Sea \(\displaystyle S_n(x)=\frac{a_0}{2}+\sum_{j=1}^{n-1} (a_j\cos(jx)+b_j\sin(jx))+a_n\cos(nx),\) la expresión que minimiza el error \[ \begin{align*} & E(a_0,a_1,\ldots,a_n,b_1,b_2,\ldots,b_{n-1})= \int_{-\pi}^{\pi}(f(x)-S_n(x))^2\, dx\\ &\quad = \int_{-\pi}^{\pi}\left(f(x)-\left(\frac{a_0}{2}+\sum_{j=1}^{n} a_j\cos(jx)+\sum_{j=1}^{n-1}b_j\sin(jx))\right)\right)^2\, dx. \end{align*} \]

Proposición. Expresión del término del error (continuación)

es decir: \[ \begin{align*} a_k=& \frac{1}{\pi}\int_{-\pi}^\pi f(x)\cos(kx)\, dx,\ k=0,1,\ldots,n,\\ b_k= &\frac{1}{\pi}\int_{-\pi}^\pi f(x)\sin(kx)\, dx,\ k=1,\ldots,n-1. \end{align*} \] Entonces la expresión del término del error puede escribirse de la forma siguiente: \[ E(a_i,b_j)=\int_{-\pi}^\pi f(x)^2\, dx-\pi\left(\frac{a_0^2}{2}+\sum_{k=1}^n a_k^2+\sum_{k=1}^{n-1}b_k^2\right). \]

En la sección anterior, vimos cómo aproximar una función \(f\) integrable en \([-\pi,\pi]\) por una función trigonométrica \[ S_n(x)=\frac{a_0}{2}+\sum_{j=1}^n a_j\cos(jx)+\sum_{j=1}^{n-1} b_j\sin(jx), \] hallando los coeficientes \(a_j\), \(j=0,1,\ldots,n\) y \(b_j\), \(j=1,\ldots,n-1\) que minimizan la función error: \[ \begin{align*} & E(a_0,a_1,\ldots,a_n,b_1,b_2,\ldots,b_{n-1})= \int_{-\pi}^{\pi}(f(x)-S_n(x))^2\, dx\\ &\quad = \int_{-\pi}^{\pi}\left(f(x)-\left(\frac{a_0}{2}+\sum_{j=1}^{n} a_j\cos(jx)+\sum_{j=1}^{n-1}b_j\sin(jx))\right)\right)^2\, dx. \end{align*} \] En esta sección vamos a ver la versión discreta de dicha aproximación.

Concretamente, sea \(m\) un entero positivo. Consideramos los \(2m\) puntos siguientes equiespaciados en el intervalo \([-\pi,\pi]\): \[ x_j =-\pi +\frac{j\pi}{m},\ j=0,\ldots,2m-1. \] Supongamos que nos dan \(y_j\), \(j=0,1,\ldots,2m-1\) las \(2m\) ordenadas correspondientes a dichos \(2m\) puntos.

Es decir, los puntos en el plano \((x_j,y_j)\), \(j=0,1,\ldots,2m-1\) forman nuestra nube de puntos.

Nos planteamos el problema siguiente:

Problema de minimización trigonométrica discreta.

Queremos hallar una función \(S_n(x)\in T_n\), con \(n<m\), \[ S_n(x)=\frac{a_0}{2}+\sum_{j=1}^n a_j\cos(jx)+\sum_{j=1}^{n-1} b_j\sin(jx), \] que minimice la función error siguiente: \[ \begin{align*} & E(a_0,a_1,\ldots,a_n,b_1,b_2,\ldots,b_{n-1})= \sum_{i=0}^{2m-1}(y_i-S_n(x_i))^2\\ &\quad = \sum_{i=0}^{2m-1}\left(y_i-\left(\frac{a_0}{2}+\sum_{j=1}^{n} a_j\cos(jx_i)+\sum_{j=1}^{n-1}b_j\sin(jx_i))\right)\right)^2\,. \end{align*} \]

Observación.

La condición \(n<m\) se impone ya que si \(n=m\), basta hallar la función trigonométrica \(S_n(x)\) que interpole la nube de puntos \((x_k,y_k),\ k=0,1,\ldots,2m-1\), es decir \(S_n(x_k)=y_k,\ k=0,1,\ldots,2m-1\) y la función error sería mínima al ser nula: \(E(a_i,b_j)=0\).

En el caso \(n=m\), tenemos \(2m\) ecuaciones \(S_n(x_k)=y_k,\ k=0,1,\ldots,2m-1\) y \(2n=2m\) incógnitas \(a_0,a_1,\ldots,a_n,b_1,\ldots,b_{n-1}\), (\(n+1+n-1=2n\)).

El Teorema siguiente nos da la función \(S_n(x)\) buscada:

Teorema: función que minimiza la función error.

La función \(S_n(x)\) que minimiza la función error es la siguiente: \[ S_n(x)=\frac{a_0}{2}+\sum_{j=1}^n a_j\cos(jx)+\sum_{j=1}^{n-1} b_j\sin(jx), \] con: \[ \begin{align*} a_k = & \frac{1}{m}\sum_{i=0}^{2m-1}y_i\cos(kx_i),\ k=0,1,\ldots,n,\\ b_k = & \frac{1}{m}\sum_{i=0}^{2m-1}y_i\sin(kx_i),\ k=1,2,\ldots,n-1. \end{align*} \]

De forma parecida al caso continuo, necesitamos el lema siguiente para la demostración del Teorema:

Lema.

Sean \(k\), \(l\) y \(m\) tres enteros. Entonces se cumplen las igualdades siguientes:

• si \(k\) no es múltiplo de \(2m\), \[ \sum_{j=0}^{2m-1}\cos(kx_j)= 0,\quad \sum_{j=0}^{2m-1}\sin(kx_j)=0. \]
• si \(k\) no es múltiplo de \(m\), \[ \sum_{j=0}^{2m-1}\cos(kx_j)^2 = m,\quad \sum_{j=0}^{2m-1}\sin^2(kx_j)=m. \]

Lema (continuación).

• si \(k+l\) y \(l-k\) no son múltiplos de \(2m\), \[ \sum_{j=0}^{2m-1}\cos(kx_j)\sin(lx_j) = 0. \]

• si \(k+l\) y \(l-k\) no son múltiplos de \(2m\) y \(k\neq l\), \[ \sum_{j=0}^{2m-1}\cos(kx_j)\cos(lx_j) = 0,\quad \sum_{j=0}^{2m-1}\sin(kx_j)\sin(lx_j) = 0. \]

Caso 1: Supondremos que el lenguaje de programación en el que se va a programar el código empieza a inicializar en \(1\):

• INPUT número de puntos de la nube \(m\), puntos de la nube \((x_k,y_k)\), \(k=1,\ldots,2\cdot m\), valor \(n\).
• Set \(\mathbf{a}=\mathbf{0}\). (Inicializamos los vectores \(\mathbf{a}\) y \(\mathbf{b}\) donde guardaremos \(a_k\), \(k=1,\ldots,n+1\) y \(b_k,\ k=1,\ldots,n-1\).)
• Set \(\mathbf{b}=\mathbf{0}\).
• For i=1,...,2*m (primero calculamos \(a_1\))
  • Set \(a_1=a_1+y_i\)
• Set \(a_1=\frac{a_1}{m}\)

• For k=2,...,n (Calculamos \(a_{k}\), las componentes del vector \(\mathbf{a}\) y las componentes del vector \(\mathbf{b}\), para \(k=2,\ldots,n\))
  • For i=1,...,2*m
    • Set \(a_k=a_k+y_i\cdot\cos((k-1)\cdot x_i)\)
    • Set \(b_{k-1}=b_{k-1}+y_i\cdot\sin((k-1)\cdot x_i)\)
  • Set \(a_k =\frac{a_k}{m}\)
  • Set \(b_k =\frac{b_k}{m}\)
• For i=1,...,2*m (Calculamos \(a_{n+1}\))
  • Set \(a_{n+1}=a_{n+1}+y_i\cdot\cos(n\cdot x_i)\)
• Set \(a_{n+1}=\frac{a_{n+1}}{m}\).
• Print \(a_1,a_2,\ldots,a_{n+1},b_1,\ldots,b_{n-1}\).

Caso 2: Supondremos que el lenguaje de programación en el que se va a programar el código empieza a inicializar en \(0\) como Python:

• INPUT número de puntos de la nube \(m\), puntos de la nube \((x_k,y_k)\), \(k=0,\ldots,2\cdot m-1\), valor \(n\).
• Set \(\mathbf{a}=\mathbf{0}\). (Inicializamos los vectores \(\mathbf{a}\) y \(\mathbf{b}\) donde guardaremos \(a_k\), \(k=0,\ldots,n\) y \(b_k,\ k=1,\ldots,n-1\).)
• Set \(\mathbf{b}=\mathbf{0}\).
• For i=0,...,2*m-1 (primero calculamos \(a_0\))
  • Set \(a_0=a_0+y_i\)
• Set \(a_0=\frac{a_0}{m}\)

• For k=1,...,n-1 (Calculamos \(a_{k}\), las componentes del vector \(\mathbf{a}\) y las componentes del vector \(\mathbf{b}\), para \(k=1,\ldots,n-1\))
  • For i=0,...,2*m-1
    • Set \(a_k=a_k+y_i\cdot\cos(k\cdot x_i)\)
    • Set \(b_k=b_k+y_i\cdot\sin(k\cdot x_i)\)
  • Set \(a_k =\frac{a_k}{m}\)
  • Set \(b_k =\frac{b_k}{m}\)
• For i=0,...,2*m-1 (Calculamos \(a_{n}\))
  • Set \(a_{n}=a_{n}+y_i\cdot\cos(n\cdot x_i)\)
• Set \(a_{n}=\frac{a_{n}}{m}\).
• Print \(a_0,a_1,\ldots,a_{n},b_1,\ldots,b_{n-1}\).

Cuando en la aproximación trigonométrica discreta, el valor de \(m\) y \(n\) coinciden, \(n=m\), el valor de la aproximación \(S_n(x)\) se denomina serie de Fourier discreta correspondiente a la nube de puntos \((x_k,y_k),\ k=0,\ldots,2m-1\): \[ S_n(x)=\frac{a_0}{2}+\sum_{j=1}^n a_j\cos(jx)+\sum_{j=1}^{n-1} b_j\sin(jx). \] En este caso, la función \(S_n(x)\) interpola a los puntos \((x_k,y_k)\), \(k=0,1,\ldots,2m-1\) ya que tenemos \(n+1+n-1=2n\) coeficientes \(a_k,\ k=0,\ldots,n\) y \(b_k,\ k=1,\ldots,n-1\) y \(2m=2n\) condiciones: \[ S_n(x_k)=y_k,\ k=0,1,\ldots,2m-1=2n-1. \]

Ahora bien, las expresiones de las \(a_k\), \(k=0,1,\ldots,n\), \(b_k\), \(k=1,\ldots,n-1\) son las que hemos deducido anteriormente: \[ \begin{align*} a_k = & \frac{1}{m}\sum_{i=0}^{2m-1}y_i\cos(kx_i),\ k=0,1,\ldots,n-1,\\ b_k = & \frac{1}{m}\sum_{i=0}^{2m-1}y_i\sin(kx_i),\ k=0,1,\ldots,n-1, \end{align*} \] excepto la expresión para \(a_n\).

Para deducir la expresión para \(a_n\), hacíamos lo siguiente: \[ \begin{align*} \frac{\partial E}{\partial a_n}= &-2\sum_{i=0}^{2m-1}\left(\left(y_i-\left(\frac{a_0}{2}+ \sum_{j=1}^{n}a_j\cos(jx_i)+\sum_{j=1}^{n-1}b_j\sin(jx_i)\right)\right)\cos(nx_i)\right) \\ = & -2\left(\sum_{i=0}^{2m-1} y_i\cos(nx_i)-\frac{a_0}{2}\sum_{i=0}^{2m-1}\cos(nx_i)-\sum_{j=1}^{n} a_j\sum_{i=0}^{2m-1}\cos(jx_i)\cos(nx_i)-\sum_{j=1}^{n-1} b_j\sum_{i=0}^{2m-1}\sin(jx_i)\cos(nx_i)\right)\\ = & -2\left(\sum_{i=0}^{2m-1} y_i\cos(nx_i)-a_n\sum_{i=0}^{2m-1}\cos^2(nx_i)\right). \end{align*} \] Nos falta calcular \(\displaystyle\sum_{i=0}^{2m-1}\cos^2(nx_i)=\sum_{i=0}^{2m-1}\cos^2(mx_i)\) pero no podemos aplicar el lema ya que recordemos que el lema nos da la suma \(\displaystyle\sum_{i=0}^{2m-1}\cos^2(k x_i)=m\), si \(k\) no es múltiplo de \(m\) pero en nuestro caso \(k=n\) y como \(n=m\), \(k=m\) es múltiplo de \(m\).

El lema siguiente nos calcula la suma solicitada:

Lema. \[ \sum_{i=0}^{2m-1}\cos^2(mx_i) =2m. \]

Entonces: \[ \begin{align*} \frac{\partial E}{\partial a_n}= & -2\left(\sum_{i=0}^{2m-1} y_i\cos(nx_i)-a_n\sum_{i=0}^{2m-1}\cos^2(nx_i)\right)=0,\ \Rightarrow \\ \sum_{i=0}^{2m-1} y_i\cos(nx_i)-2ma_n= & 0. \end{align*} \] El valor de \(a_n\) será, pues, · \[ a_n=\frac{1}{2m}\sum_{i=0}^{2m-1}y_i\cos(n x_i). \]

La función trigonométrica interpoladora o la expresión de la serie de Fourier discreta correspondiente a la nube de puntos \((x_k,y_k)\), \(k=0,1,\ldots,2m-1\) será: \[ S_m(x)=\frac{a_0}{2}+\sum_{j=1}^m a_j\cos(jx)+\sum_{j=1}^{m-1}b_j\sin(jx), \] con \[ \begin{align*} a_k = & \frac{1}{m}\sum_{i=0}^{2m-1}y_i\cos(kx_i),\ k=0,1,\ldots,m-1,\ a_m = \frac{1}{2m}\sum_{i=0}^{2m-1}y_i\cos(m x_i),\\ b_k = & \frac{1}{m}\sum_{i=0}^{2m-1}y_i\sin(kx_i),\ k=0,1,\ldots,m-1,\\ \end{align*} \] donde sólo hemos escrito el valor \(m\) ya que suponemos que \(n=m\).

De cara a tener una expresión más compacta en la expresión de los coeficientes \(a_k\), redefinimos \(a_m\) como \(\displaystyle a_m=\frac{1}{m}\sum_{i=0}^{2m-1}y_i\cos(m x_i)\) y escribimos la serie de Fourier discreta como: \[ S_m(x)=\frac{a_0+a_m\cos(mx)}{2}+\sum_{j=1}^{m-1} a_j\cos(jx)+\sum_{j=1}^{m-1}b_j\sin(jx), \] con \[ \begin{align*} a_k = & \frac{1}{m}\sum_{i=0}^{2m-1}y_i\cos(kx_i),\ k=0,1,\ldots,m,\\ b_k = & \frac{1}{m}\sum_{i=0}^{2m-1}y_i\sin(kx_i),\ k=0,1,\ldots,m-1.\\ \end{align*} \]

Calcular los valores \(a_k\), \(k=0,\ldots,m-1\) y \(b_k\), \(k=1,\ldots,m-1\) requiere en total: \((2m)^2=4m^2\) multiplicaciones/divisiones y \((2m)^2=4m^2\) sumas/restas ya que:

El número total de operaciones es de orden \(m^2\): \(O(m^2)\).

Vamos a diseñar un algoritmo que reduce el número de operaciones a \(O(m\ln_2(m))\). Dicho algoritmo se denomina transformada rápida de Fourier.

En primer lugar, consideramos la función siguiente: \[ \frac{1}{m}\sum_{k=0}^{2m-1}c_k \mathrm{e}^{\mathrm{i}kx},\mbox{ con }c_k=\sum_{j=0}^{2m-1}y_j\mathrm{e}^{\frac{\mathrm{i}k\pi j}{m}},\ k=0,1,\ldots,2m-1. \]

Veamos a continuación que los coeficientes \(a_k\) y \(b_k\) se pueden calcular a partir de los coeficientes \(c_k\):

Proposición.

La relación entre los coeficientes \(c_k\), \(a_k\) y \(b_k\) es la siguiente: \[ a_k+\mathrm{i}b_k=\frac{(-1)^k}{m}c_k,\ k=0,1,\ldots,m, \] donde suponemos que \(b_0=b_m=0\).

Por tanto, si \(k=0\), la relación anterior nos dice que \(a_0 =\frac{(-1)^0}{m}c_0=\frac{c_0}{m}\) y si \(k=m\), \(a_m=\frac{(-1)^m}{m}c_m\).

A partir de ahora, supongamos que \(m\) es un potencia de \(2\), es decir, que existe un entero \(p\) tal que \(m=2^p\).

Para empezar, calculemos para \(k=0,1,\ldots,m-1\), la suma siguiente: \[ c_k+c_{m+k}=\sum_{j=0}^{2m-1}y_j\mathrm{e}^{\frac{\mathrm{i}k\pi j}{m}}+\sum_{j=0}^{2m-1}y_j\mathrm{e}^{\frac{\mathrm{i}(m+k)\pi j}{m}}=\sum_{j=0}^{2m-1}y_j\mathrm{e}^{\frac{\mathrm{i}k\pi j}{m}}\left(1+\mathrm{e}^{\mathrm{i}\pi j}\right). \] El valor de \(1+\mathrm{e}^{\mathrm{i}\pi j}\) vale: \[ 1+\mathrm{e}^{\mathrm{i}\pi j} =\begin{cases}2,& \mbox{ si } j \mbox{ es par}\\ 0, & \mbox{ si } j\mbox{ es impar.}\end{cases} \]

Entonces en la suma \(\displaystyle \sum_{j=0}^{2m-1}y_j\mathrm{e}^{\frac{\mathrm{i}k\pi j}{m}}\left(1+\mathrm{e}^{\mathrm{i}\pi j}\right)\) que aparece en la expresión de \(c_k+c_{m+k}\) sólo se suman los términos en los que \(j\) es par, es decir: \[ c_k + c_{m+k}=2\sum_{l=0}^{m-1}y_{2l}\mathrm{e}^{\frac{\mathrm{i}\pi k(2l)}{m}}=2\sum_{l=0}^{m-1}y_{2l}\mathrm{e}^{\frac{\mathrm{i}\pi kl}{\frac{m}{2}}}, \] donde hemos hecho el cambio de índice \(j=2l\) y como \(j=0,1,\ldots,2m-1\), \(l=0,1,\ldots,m-1\).

A continuación, calculamos la suma siguiente: \[ c_k-c_{m+k}=\sum_{j=0}^{2m-1}y_j\mathrm{e}^{\frac{\mathrm{i}k\pi j}{m}}-\sum_{j=0}^{2m-1}y_j\mathrm{e}^{\frac{\mathrm{i}(m+k)\pi j}{m}}=\sum_{j=0}^{2m-1}y_j\mathrm{e}^{\frac{\mathrm{i}k\pi j}{m}}\left(1-\mathrm{e}^{\mathrm{i}\pi j}\right). \] El valor de \(1-\mathrm{e}^{\mathrm{i}\pi j}\) vale: \[ 1-\mathrm{e}^{\mathrm{i}\pi j} =\begin{cases}0,& \mbox{ si } j \mbox{ es par}\\ 2, & \mbox{ si } j\mbox{ es impar.}\end{cases} \]

Entonces en la suma \(\displaystyle \sum_{j=0}^{2m-1}y_j\mathrm{e}^{\frac{\mathrm{i}k\pi j}{m}}\left(1-\mathrm{e}^{\mathrm{i}\pi j}\right)\) que aparece en la expresión de \(c_k-c_{m+k}\) sólo se suman los términos en los que \(j\) es impar, es decir: \[ c_k - c_{m+k}=2\sum_{l=0}^{m-1}y_{2l+1}\mathrm{e}^{\frac{\mathrm{i}\pi k(2l+1)}{m}}=2\mathrm{e}^{\frac{\mathrm{i}k\pi}{m}}\sum_{l=0}^{m-1}y_{2l+1}\mathrm{e}^{\frac{\mathrm{i}\pi kl}{\frac{m}{2}}}, \] donde hemos hecho el cambio de índice \(j=2l+1\) y como \(j=0,1,\ldots,2m-1\), \(l=0,1,\ldots,m-1\).

En resumen, hemos deducido lo siguiente: \[ \begin{align*} c_k +c_{m+k}=& 2\sum_{l=0}^{m-1}y_{2l}\mathrm{e}^{\frac{\mathrm{i}\pi kl}{\frac{m}{2}}},\\ c_k -c_{m+k}= & 2\mathrm{e}^{\frac{\mathrm{i}k\pi}{m}}\sum_{l=0}^{m-1}y_{2l+1}\mathrm{e}^{\frac{\mathrm{i}\pi kl}{\frac{m}{2}}}. \end{align*} \]

Ahora viene la idea clave: las sumas \(\displaystyle\sum_{l=0}^{m-1}y_{2l}\mathrm{e}^{\frac{\mathrm{i}\pi kl}{\frac{m}{2}}}\) y \(\displaystyle\sum_{l=0}^{m-1}y_{2l+1}\mathrm{e}^{\frac{\mathrm{i}\pi kl}{\frac{m}{2}}}\) corresponden a la definición de los coeficientes \(c_k\) para las muestras \(y_0,y_2,\ldots,y_{2m-2}\) (términos pares) y \(y_1,y_3,\ldots,y_{2m-1}\) (términos impares), respectivamente.

Es decir, si sustituimos \(m\) por \(\frac{m}{2}\) en la definición de los \(c_k\) obtenemos los nuevos \(c_k\) que aparecen en los términos de la derecha en las expresiones de \(c_k+c_{m+k}\) y \(c_{k}-c_{m+k}\).

Introduzcamos un poco de notación:

• llamamos \(c_k^{(2m)}\) a los coeficientes \(c_k\) para la muestra \(y_0,y_1,\ldots,y_{2m-1}\), y
• llamamos \(c_{k,p}^{(m)}\) y \(c_{k,i}^{(m)}\) a los coeficientes \(c_k\) para las muestras \(y_0,y_2,\ldots,y_{2m-2}\) (términos pares) y \(y_1,y_3,\ldots,y_{2m-1}\) (términos impares) respectivamente.

La relación anterior sería: \[ \begin{align*} c_k^{(2m)} +c_{m+k}^{(2m)}=& 2c_{k,p}^{(m)},\\ c_k^{(2m)} -c_{m+k}^{(2m)}= & 2\mathrm{e}^{\frac{\mathrm{i}k\pi}{m}}c_{k,i}^{(m)}, \end{align*} \] para \(k=0,\ldots,m-1\).

La solución del sistema anterior, es decir, el valor de \(c_k^{(2m)}\) y \(c_{m+k}^{(2m)}\) en función de \(c_{k,p}^{(m)}\) y \(c_{k,i}^{(m)}\) es la siguiente: \[ c_k^{(2m)}=c_{k,p}^{(m)}+\mathrm{e}^{\frac{\mathrm{i}k\pi}{m}}c_{k,i}^{(m)},\quad c_{m+k}^{(2m)}=c_{k,p}^{(m)}-\mathrm{e}^{\frac{\mathrm{i}k\pi}{m}}c_{k,i}^{(m)},\ k=0,1,\ldots,m-1. \] La expresión anterior permite hallar los coeficientes \(c_k^{(2m)}\), \(c_{m+k}^{(2m)}\) para la muestra \(y_0,y_1,\ldots,y_{2m-1}\) de longitud \(2m\) en función de los coeficientes \(c_{k,p}^{(m)}\) y \(c_{k,i}^{(m)}\) correspondientes a las muestras \(y_0,y_2,\ldots,y_{2m-2}\) y \(y_1,y_3,\ldots,y_{2m-1}\), respectivamente, ambas de longitud \(m\).

El algoritmo de la transformada rápida de Fourier es ir cada vez reduciendo el tamaño de la muestra hasta llegar a una muestra de tamaño \(2\), calcular los coeficientes \(c_k\) correspondiente a esta muestra e ir doblando la muestra aplicando la expresión anterior hasta llegar a la muestra inicial \(y_0,y_1,\ldots,y_{2m-1}\) de tamaño \(2m=2^{p+1}\).

Calculemos los coeficientes \(c_0^{(2)}\) y \(c_1^{(2)}\) para una muestra de tamaño \(2\), \(y_0,y_1\). \[ \begin{align*} c_0^{(2)} = & \left(y_0\mathrm{e}^{\frac{\mathrm{i}\pi\cdot 0\cdot 0}{1}}+y_1\mathrm{e}^{\frac{\mathrm{i}\pi\cdot 1\cdot 0}{1}}\right)=y_0+y_1,\\ c_1^{(2)} = & \left(y_0\mathrm{e}^{\frac{\mathrm{i}\pi\cdot 0\cdot 1}{1}}+y_1\mathrm{e}^{\frac{\mathrm{i}\pi\cdot 1\cdot 1}{1}}\right)=y_0+y_1\mathrm{e}^{\mathrm{i}\pi}=y_0-y_1. \end{align*} \] Veamos cuantas operaciones hemos de realizar si usamos el algoritmo de la transformada rápida de Fourier.

Consideremos la expresión que nos da los coeficientes \(c_{k}^{(2m)}\) y \(c_{m+k}^{(2m)}\) en función de los coeficientes \(c_{k,p}^{(m)}\) y \(c_{k,i}^{(m)}\), \(k=0,1,\ldots,m-1\):

• Número de multiplicaciones/divisiones: \(2m\).
• Número de sumas/restas: \(2m\).

La tabla siguiente nos muestra el número de operaciones en todos los pasos:

Como puede verse en la tabla anterior hay en total \(p\) pasos y en el último paso, es decir, el paso para \(m=1\) en el que hay \(m\) muestras de tamaño \(2\), hay \(0\cdot m=m\) multiplicaciones y \(2\cdot m=2m\) sumas/restas.

Por tanto,

• Número total de multiplicaciones/divisiones: \(0+2m\cdot p=2m\cdot \ln_2 m=O(m\ln_2 m)\).
• Número total de sumas/restas: \(2m+2m\cdot p=2m(1+\ln_2 m)=O(m\ln_2 m)\).

En resumen, hemos reducido el número total de operaciones de \(O(m^2)\) a \(O(m\ln_2 m)\).

Demos el algoritmo detallado para \(m=4=2^2\), y por tanto, \(p=2\).

• Paso I: lo primero que tenemos que hacer es reordenar la muestra de \(y_k\), \(k=0\ldots,7\) para poder aplicar el algoritmo:
  • Paso \(\frac{m}{2}=2\longrightarrow m=4\):
    • muestra 1: \(y_0,y_2,y_4,y_6\),
    • muestra 2: \(y_1,y_3,y_5,y_7\).
  • Muestras de longitud \(2\):
    • muestra 1: \(y_0,y_4\),
    • muestra 2: \(y_2,y_6\),
    • muestra 3: \(y_1,y_5\),
    • muestra 4: \(y_3,y_7\).

• Paso II: calculamos los coeficientes para las muestras \(1\), \(2\), \(3\) y \(4\). \[ \begin{align*} c_{0,1}^{(2)}=&y_0+y_4,\quad c_{1,1}^{(2)}=y_0-y_4,\\ c_{0,2}^{(2)}=&y_2+y_6,\quad c_{1,2}^{(2)}=y_2-y_6,\\ c_{0,3}^{(2)}=&y_1+y_5,\quad c_{1,3}^{(2)}=y_1-y_5,\\ c_{0,4}^{(2)}=&y_3+y_7,\quad c_{1,4}^{(2)}=y_3-y_7. \end{align*} \]

• Paso III: calculamos los coeficientes para las muestras \(1\) y \(2\) en el paso \(\frac{m}{2}=2\longrightarrow m=4\):
  • Muestra 1: \(y_0,y_2,y_4,y_6\). Usamos las expresiones siguientes: \[ \begin{align*} c_{0,1}^{(4)}=&c_{0,1}^{(2)}+\mathrm{e}^{\frac{\mathrm{i}\cdot 0\pi}{2}}c_{0,2}^{(2)}=c_{0,1}^{(2)}+c_{0,2}^{(2)}=y_0+y_4+y_2+y_6\\ = & y_0+y_2+y_4+y_6,\\ c_{1,1}^{(4)}=&c_{1,1}^{(2)}+\mathrm{e}^{\frac{\mathrm{i}\pi}{2}}c_{1,2}^{(2)}=c_{1,1}^{(2)}+\mathrm{i}c_{1,2}^{(2)}=y_0-y_4+\mathrm{i}(y_2-y_6)\\ = & y_0+\mathrm{i}y_2-y_4-\mathrm{i}y_6,\\ c_{2,1}^{(4)}=&c_{0,1}^{(2)}-\mathrm{e}^{\frac{\mathrm{i}\cdot 0\pi}{2}}c_{0,2}^{(2)}=c_{0,1}^{(2)}-c_{0,2}^{(2)}=y_0+y_4-(y_2+y_6)\\ = & y_0-y_2+y_4-y_6,\\ c_{3,1}^{(4)}=&c_{1,1}^{(2)}-\mathrm{e}^{\frac{\mathrm{i}\pi}{2}}c_{1,2}^{(2)}=c_{1,1}^{(2)}-\mathrm{i}c_{1,2}^{(2)}=y_0-y_4-\mathrm{i}(y_2-y_6)\\ = & y_0-\mathrm{i}y_2-y_4+\mathrm{i}y_6, \end{align*} \]

• Paso III (continuación): calculamos los coeficientes para las muestras \(1\) y \(2\) en el paso \(\frac{m}{2}=2\longrightarrow m=4\):
  • Muestra 2: \(y_1,y_3,y_5,y_7\). Usamos las expresiones siguientes: \[ \begin{align*} c_{0,2}^{(4)}=&c_{0,3}^{(2)}+\mathrm{e}^{\frac{\mathrm{i}\cdot 0\pi}{2}}c_{0,4}^{(2)}=c_{0,3}^{(2)}+c_{0,4}^{(2)}=y_1+y_5+y_3+y_7\\ = & y_1+y_3+y_5+y_7,\\ c_{1,2}^{(4)}=&c_{1,3}^{(2)}+\mathrm{e}^{\frac{\mathrm{i}\pi}{2}}c_{1,4}^{(2)}=c_{1,3}^{(2)}+\mathrm{i}c_{1,4}^{(2)}=y_1-y_5+\mathrm{i}(y_3-y_7)\\ = & y_1+\mathrm{i}y_3-y_5-\mathrm{i}y_7,\\ c_{2,2}^{(4)}=&c_{0,3}^{(2)}-\mathrm{e}^{\frac{\mathrm{i}\cdot 0\pi}{2}}c_{0,4}^{(2)}=c_{0,3}^{(2)}-c_{0,4}^{(2)}=y_1+y_5-(y_3+y_7)\\ = & y_1-y_3+y_5-y_7,\\ c_{3,2}^{(4)}=&c_{1,3}^{(2)}-\mathrm{e}^{\frac{\mathrm{i}\pi}{2}}c_{1,4}^{(2)}=c_{1,3}^{(2)}-\mathrm{i}c_{1,4}^{(2)}=y_1-y_5-\mathrm{i}(y_3-y_7)\\ = & y_1-\mathrm{i}y_3-y_5 +\mathrm{i}y_7. \end{align*} \]

• Paso IV: calculamos los coeficientes para la muestra original \(y_0,y_1,y_2,y_3,y_4,y_5,y_6,y_7\) en el paso \(m=4\longrightarrow 2m=8\): \[ \begin{align*} c_{0}^{(8)}=& c_{0,1}^{(4)}+\mathrm{e}^{\frac{\mathrm{i}\cdot 0\pi}{4}}c_{0,2}^{(4)}= c_{0,1}^{(4)}+c_{0,2}^{(4)}\\ = &y_0+y_2+y_4+y_6+y_1+y_5+y_3+y_7\\ =& y_0+y_1+y_2+y_3+y_4+y_5+y_6+y_7,\\ c_{1}^{(8)}=& c_{1,1}^{(4)}+\mathrm{e}^{\frac{\mathrm{i}\cdot 1\pi}{4}}c_{1,2}^{(4)}= c_{1,1}^{(4)}+\left(\frac{\sqrt{2}}{2}+\frac{\sqrt{2}}{2}\mathrm{i}\right)c_{1,2}^{(4)}\\ = &y_0+\mathrm{i}y_2-y_4-\mathrm{i}y_6+\left(\frac{\sqrt{2}}{2}+\frac{\sqrt{2}}{2}\mathrm{i}\right)(y_1+\mathrm{i}y_3-y_5-\mathrm{i}y_7)\\ =& y_0+\left(\frac{\sqrt{2}}{2}+\frac{\sqrt{2}}{2}\mathrm{i}\right)y_1+\mathrm{i}y_2+\left(\frac{-\sqrt{2}}{2}+\frac{\sqrt{2}}{2}\mathrm{i}\right)y_3-y_4\\ & \quad -\left(\frac{\sqrt{2}}{2}+\frac{\sqrt{2}}{2}\mathrm{i}\right)y_5-\mathrm{i}y_6+\left(\frac{\sqrt{2}}{2}-\frac{\sqrt{2}}{2}\mathrm{i}\right)y_7,\\ \end{align*} \]

• Paso IV (continuación): \[ \begin{align*} c_{2}^{(8)}=& c_{2,1}^{(4)}+\mathrm{e}^{\frac{\mathrm{i}\cdot 2\pi}{4}}c_{2,2}^{(4)}= c_{2,1}^{(4)}+\mathrm{i}c_{2,2}^{(4)}\\ = &y_0-y_2+y_4-y_6+\mathrm{i}(y_1-y_3+y_5-y_7)\\ =& y_0+\mathrm{i}y_1-y_2-\mathrm{i}y_3+y_4+\mathrm{i}y_5-y_6-\mathrm{i}y_7,\\ c_{3}^{(8)}=& c_{3,1}^{(4)}+\mathrm{e}^{\frac{\mathrm{i}\cdot 3\pi}{4}}c_{3,2}^{(4)}= c_{3,1}^{(4)}+\left(-\frac{\sqrt{2}}{2}+\frac{\sqrt{2}}{2}\mathrm{i}\right)c_{3,2}^{(4)}\\ = &y_0-\mathrm{i}y_2-y_4+\mathrm{i}y_6+\left(-\frac{\sqrt{2}}{2}+\frac{\sqrt{2}}{2}\mathrm{i}\right)(y_1-\mathrm{i}y_3-y_5+\mathrm{i}y_7)\\ =& y_0+\left(-\frac{\sqrt{2}}{2}+\frac{\sqrt{2}}{2}\mathrm{i}\right)y_1-\mathrm{i}y_2+\left(\frac{\sqrt{2}}{2}+\frac{\sqrt{2}}{2}\mathrm{i}\right)y_3-y_4\\ & \quad +\left(\frac{\sqrt{2}}{2}-\frac{\sqrt{2}}{2}\mathrm{i}\right)y_5+\mathrm{i}y_6-\left(\frac{\sqrt{2}}{2}+\frac{\sqrt{2}}{2}\mathrm{i}\right)y_7,\\ \end{align*} \]

• Paso IV (continuación): \[ \begin{align*} c_{4}^{(8)}=& c_{0,1}^{(4)}-\mathrm{e}^{\frac{\mathrm{i}\cdot 0\pi}{4}}c_{0,2}^{(4)}= c_{0,1}^{(4)}-c_{0,2}^{(4)}\\ = &y_0+y_2+y_4+y_6-(y_1+y_5+y_3+y_7)\\ =& y_0-y_1+y_2-y_3+y_4-y_5+y_6-y_7,\\ c_{5}^{(8)}=& c_{1,1}^{(4)}-\mathrm{e}^{\frac{\mathrm{i}\cdot 1\pi}{4}}c_{1,2}^{(4)}= c_{1,1}^{(4)}-\left(\frac{\sqrt{2}}{2}+\frac{\sqrt{2}}{2}\mathrm{i}\right)c_{1,2}^{(4)}\\ = &y_0+\mathrm{i}y_2-y_4-\mathrm{i}y_6-\left(\frac{\sqrt{2}}{2}+\frac{\sqrt{2}}{2}\mathrm{i}\right)(y_1+\mathrm{i}y_3-y_5-\mathrm{i}y_7)\\ =& y_0-\left(\frac{\sqrt{2}}{2}+\frac{\sqrt{2}}{2}\mathrm{i}\right)y_1+\mathrm{i}y_2+\left(\frac{\sqrt{2}}{2}-\frac{\sqrt{2}}{2}\mathrm{i}\right)y_3-y_4\\ & \quad +\left(\frac{\sqrt{2}}{2}+\frac{\sqrt{2}}{2}\mathrm{i}\right)y_5-\mathrm{i}y_6+\left(\frac{-\sqrt{2}}{2}+\frac{\sqrt{2}}{2}\mathrm{i}\right)y_7,\\ \end{align*} \]

• Paso IV (continuación): \[ \begin{align*} c_{6}^{(8)}=& c_{2,1}^{(4)}-\mathrm{e}^{\frac{\mathrm{i}\cdot 2\pi}{4}}c_{2,2}^{(4)}= c_{2,1}^{(4)}-\mathrm{i}c_{2,2}^{(4)}\\ = &y_0-y_2+y_4-y_6-\mathrm{i}(y_1-y_3+y_5-y_7)\\ =& y_0-\mathrm{i}y_1-y_2+\mathrm{i}y_3+y_4-\mathrm{i}y_5-y_6+\mathrm{i}y_7,\\ c_{7}^{(8)}=& c_{3,1}^{(4)}-\mathrm{e}^{\frac{\mathrm{i}\cdot 3\pi}{4}}c_{3,2}^{(4)}= c_{3,1}^{(4)}-\left(-\frac{\sqrt{2}}{2}+\frac{\sqrt{2}}{2}\mathrm{i}\right)c_{3,2}^{(4)}\\ = &y_0-\mathrm{i}y_2-y_4+\mathrm{i}y_6+\left(\frac{\sqrt{2}}{2}-\frac{\sqrt{2}}{2}\mathrm{i}\right)(y_1-\mathrm{i}y_3-y_5+\mathrm{i}y_7)\\ =& y_0+\left(\frac{\sqrt{2}}{2}-\frac{\sqrt{2}}{2}\mathrm{i}\right)y_1-\mathrm{i}y_2-\left(\frac{\sqrt{2}}{2}+\frac{\sqrt{2}}{2}\mathrm{i}\right)y_3-y_4\\ & \quad +\left(-\frac{\sqrt{2}}{2}+\frac{\sqrt{2}}{2}\mathrm{i}\right)y_5+\mathrm{i}y_6+\left(\frac{\sqrt{2}}{2}+\frac{\sqrt{2}}{2}\mathrm{i}\right)y_7,\\ \end{align*} \]

El gráfico siguiente muestra esquematizado el paso III:

El gráfico siguiente muestra esquematizado el paso IV:

Caso 1: Supondremos que el lenguaje de programación en el que se va a programar el código empieza a inicializar en \(1\):

En primer lugar, tenemos que definir la rutina que nos realiza el paso \(\frac{m}{2^{j+1}}\longrightarrow \frac{m}{2^{j}}\).

Le llamamos SubirNivel que depende de un vector \(y\), paso \(j\) y del vector c donde se guardan los coeficientes \(c_k\):

• SubirNivel = Rutina(y,j,c)
  • Set \(L=\mathrm{longitud}(\mathbf{c})\) (Calculamos la longitud del vector donde se guardan los coeficientes \(c_k\))
  • Set \(\mathbf{cn} = \mathbf{c}\) (En el vector \(\mathbf{cn}\) guardaremos los nuevos coeficientes \(c_k\))

• • For k=1,...,L/2 (Aplicamos la expresión hallada para los nuevos coeficientes \(c_k\) que, recordemos guardaremos en el vector \(\mathbf{cn}\))
    • Set \(cn_k = c_k+\mathrm{e}^{\frac{(k-1)\mathrm{i}\cdot\pi}{2^j}}\cdot c_{k+L/2}\)
    • Set \(cn_{L/2+k}=c_k-\mathrm{e}^{\frac{(k-1)\mathrm{i}\cdot\pi}{2^j}}\cdot c_{k+L/2}\)
  • Set \(\mathbf{z}=\mathbf{0}\) (Inicializamos el vector \(\mathbf{z}\))
  • Set \(z[\mathrm{componentes\ pares}]=\mathbf{y}[1:L/2]\) (Las componentes pares del vector \(z\) son las primeras \(L/2\) componentes del vector \(\mathbf{y}\))
  • Set \(z[\mathrm{componentes\ impares}]=\mathbf{y}[(L/2+1):L]\) (Las componentes impares del vector \(z\) son las últimas \(L/2\) componentes del vector \(\mathbf{y}\))
  • Return \(z\), \(\mathbf{cn}\) (Devolvemos el vector \(y\) arreglado y los nuevos coeficientes guardados en el vector \(\mathbf{cn}\))

• INPUT m, valores \(y_k\), \(k=1,\ldots,2 m\).
• Set \(p=\ln_2(m)\) (Calculamos \(p\) tal que \(m=2^p\).)
• Set \(z=y\). (Empezamos a reordenar los valores \(y_k\), \(k=1,\ldots,2m\) para poder aplicar el algoritmo. Primeramente hacemos una copia de \(y\) en el vector \(z\) que será el vector que manipularemos.)

• For i=1,...p
  • For j=1,...,2^(p-1)
    • Set \(aux=z[(1+(j-1)\cdot 2^{p+2-i}):(j\cdot 2^{p+2-i})]\) (Guardamos las componentes \(1+(j-1)\cdot 2^{p+2-i},\ldots,j\cdot 2^{p+2-i}\) del vector \(z\) en el vector auxiliar \(aux\).)
    • Set \(aux=aux[\mathrm{componentes\ pares}],aux[\mathrm{componentes\ impares}]\) (Separamos las componentes pares de las impares en el vector \(aux\) escribiendo primero las componentes pares seguidas de las impares.)
    • Set \(z[(1+(j-1)\cdot 2^{p+2-i}):(j\cdot 2^{p+2-i})]=aux[1:2^{p+2-i}]\) (Guardamos las primeras \(2^{p+2-i}\) componentes del vector \(aux\) en las correspondientes componentes del vector \(z\).)
• Set \(y=z\) (Guardamos el nuevo vector \(z\) en el vector permutado \(y\).)

• Set \(\mathbf{c}^{(2)}=\mathbf{0}\) (Inicializamos los coeficientes \(c_k\) para muestras de longitud \(2\), \(c_k^{(2)}\), a \(0\) donde el vector \(\mathbf{0}\) tiene \(2m\) componentes)
• For i=0,...,m-1 (Calculamos los \(c_k\) iniciales, es decir, los \(c_k^{(2)}\) para las muestras de longitud \(2\))
  • Set \(c_{2i+1}^{(2)}=y_{2i+1}+y_{2i+2}\)
  • Set \(c_{2i+2}^{(2)}=y_{2i+1}-y_{2i+2}\)
• For j=1,...,2^(p-1) (Vamos a realizar el primer paso \(2\longrightarrow 4\))
  • Set \(\mathrm{Resultado}_j=SubirNivel(y[(1+(j-1)\cdot 2^2):(j\cdot 2^2)],1,\mathbf{c}^{(2)})\) (La lista resultado tiene \(2^{p-1}\) componentes y el valor de cada una queda indicado en la expresión. Es decir, la componente \(\mathrm{Resultado}_j\) tiene dos componentes, una correspondiente al vector \(y\) y la otra correspondiente al vector \(\mathbf{c}\). )

• For i=2,...,p (Realizamos todos los demás pasos)
  • For j=1,...,2^(p-i)
    • Set \(\mathbf{yj}=[\mathrm{Resultado}_{2j-1}^{y},\mathrm{Resultado}_{2j}^{y}]\) (Definimos el vector \(\mathbf{yj}\) como la unión los vectores \(y\) de las listas \(\mathrm{Resultado}_{2j-1}\) y \(\mathrm{Resultado}_{2j}\))
    • Set \(\mathbf{cj}=[Resultado_{2j-1}^{\mathbf{c}},Resultado_{2j}^{\mathbf{c}}]\) (Idem con los vectores \(\mathbf{c}\))
    • Set \(\mathrm{Resultadon}_j =SubirNivel(\mathbf{yj},i,\mathbf{cj})\)
  • Set \(\mathrm{Resultado}=\mathrm{Resultadon}\) (Actualizamos todo)
• Print \(\mathrm{Resultado}\) ( es una lista de dos componentes: por un lado tenemos nuestra muestra \(y\) original y por otro tenemos calculados los coeficientes \(c_k\))

Caso 2: Supondremos que el lenguaje de programación en el que se va a programar el código empieza a inicializar en \(0\) como Python:

En primer lugar, tenemos que definir la rutina que nos realiza el paso \(\frac{m}{2^{j+1}}\longrightarrow \frac{m}{2^{j}}\).

Le llamamos SubirNivel que depende de un vector \(y\), paso \(j\) y del vector c donde se guardan los coeficientes \(c_k\):

• SubirNivel = Rutina(y,j,c)
  • Set \(L=\mathrm{longitud}(\mathbf{c})\) (Calculamos la longitud del vector donde se guardan los coeficientes \(c_k\))
  • Set \(\mathbf{cn} = \mathbf{c}\) (En el vector \(\mathbf{cn}\) guardaremos los nuevos coeficientes \(c_k\))

• • For k=0,...,L/2-1 (Aplicamos la expresión hallada para los nuevos coeficientes \(c_k\) que, recordemos guardaremos en el vector \(\mathbf{cn}\))
    • Set \(cn_k = c_k+\mathrm{e}^{\frac{k\mathrm{i}\cdot\pi}{2^j}}\cdot c_{k+L/2}\)
    • Set \(cn_{L/2+k}=c_k-\mathrm{e}^{\frac{k\mathrm{i}\cdot\pi}{2^j}}\cdot c_{k+L/2}\)
  • Set \(\mathbf{z}=\mathbf{0}\) (Inicializamos el vector \(\mathbf{z}\))
  • Set \(z[\mathrm{componentes\ pares}]=\mathbf{y}[0:(L/2-1)]\) (Las componentes pares del vector \(z\) son las primeras \(L/2\) componentes del vector \(\mathbf{y}\))
  • Set \(z[\mathrm{componentes\ impares}]=\mathbf{y}[(L/2):(L-1)]\) (Las componentes impares del vector \(z\) son las últimas \(L/2\) componentes del vector \(\mathbf{y}\))
  • Return \(z\), \(\mathbf{cn}\) (Devolvemos el vector \(y\) arreglado y los nuevos coeficientes guardados en el vector \(\mathbf{cn}\))

• INPUT m, valores \(y_k\), \(k=0,\ldots,2 m-1\).
• Set \(p=\ln_2(m)\) (Calculamos \(p\) tal que \(m=2^p\).)
• Set \(z=y\). (Empezamos a reordenar los valores \(y_k\), \(k=0,\ldots,2m-1\) para poder aplicar el algoritmo. Primeramente hacemos una copia de \(y\) en el vector \(z\) que será el vector que manipularemos.)

• For i=1,...p
  • For j=1,...,2^(p-1)
    • Set \(aux=z[((j-1)\cdot 2^{p+2-i}):(j\cdot 2^{p+2-i})]\) (Guardamos las componentes \((j-1)\cdot 2^{p+2-i},\ldots,j\cdot 2^{p+2-i}-1\) del vector \(z\) en el vector auxiliar \(aux\).)
    • Set \(aux=aux[\mathrm{componentes\ pares}],aux[\mathrm{componentes\ impares}]\) (Separamos las componentes pares de las impares en el vector \(aux\) escribiendo primero las componentes pares seguidas de las impares.)
    • Set \(z[((j-1)\cdot 2^{p+2-i}):(j\cdot 2^{p+2-i}-1)]=aux[0:(2^{p+2-i}-1)]\) (Guardamos las primeras \(2^{p+2-i}\) componentes del vector \(aux\) en las correspondientes componentes del vector \(z\).)
• Set \(y=z\) (Guardamos el nuevo vector \(z\) en el vector permutado \(y\).)

• Set \(\mathbf{c}^{(2)}=\mathbf{0}\) (Inicializamos los coeficientes \(c_k\) para muestras de longitud \(2\), \(c_k^{(2)}\), a \(0\) donde el vector \(\mathbf{0}\) tiene \(2m\) componentes)
• For i=0,...,m-1 (Calculamos los \(c_k\) iniciales, es decir, los \(c_k^{(2)}\) para las muestras de longitud \(2\))
  • Set \(c_{2i}^{(2)}=y_{2i}+y_{2i+1}\)
  • Set \(c_{2i+1}^{(2)}=y_{2i}-y_{2i+1}\)
• For j=1,...,2^(p-1) (Vamos a realizar el primer paso \(2\longrightarrow 4\))
  • Set \(\mathrm{Resultado}_{j-1}=SubirNivel(y[((j-1)\cdot 2^2):(j\cdot 2^2)],1,\mathbf{c}^{(2)}[((j-1)\cdot 2^2):(j\cdot 2^2)])\) (La lista resultado tiene \(2^{p-1}\) componentes y el valor de cada una queda indicado en la expresión. Es decir, la componente \(\mathrm{Resultado}_j\) tiene dos componentes, una correspondiente al vector \(y\) y la otra correspondiente al vector \(\mathbf{c}\). )

• For i=2,...,p (Realizamos todos los demás pasos)
  • For j=0,...,2^(p-i)-1
    • Set \(\mathbf{yj}=[\mathrm{Resultado}_{2j}^{y},\mathrm{Resultado}_{2j+1}^{y}]\) (Definimos el vector \(\mathbf{yj}\) como la unión los vectores \(y\) de las listas \(\mathrm{Resultado}_{2j-1}\) y \(\mathrm{Resultado}_{2j}\))
    • Set \(\mathbf{cj}=[Resultado_{2j}^{\mathbf{c}},Resultado_{2j+1}^{\mathbf{c}}]\) (Idem con los vectores \(\mathbf{c}\))
    • Set \(\mathrm{Resultadon}_j =SubirNivel(\mathbf{yj},i,\mathbf{cj})\)
  • Set \(\mathrm{Resultado}=\mathrm{Resultadon}\) (Actualizamos todo)
• Print \(\mathrm{Resultado}\) ( es una lista de dos componentes: por un lado tenemos nuestra muestra \(y\) original y por otro tenemos calculados los coeficientes \(c_k\))

Supongamos que nuestra nube de puntos \((t_k,y_k)\), \(k=0,1,\ldots,2m-1\) está en un intervalo \([a,b]\) y queremos hallar la correspondiente transformada de Fourier de dicha nube para este intervalo.

En este caso, hemos de tener en cuenta la aplicación siguiente que nos transforma valores del intervalo \([-\pi,\pi]\) al intervalo \([a,b]\) y viceversa: \[ \begin{align*} [-\pi,\pi] & \longrightarrow [a,b]\\ x & \longrightarrow t=\frac{(b-a)}{2\pi}x+\frac{a+b}{2},\\ x=\frac{2\pi}{(b-a)}\left(t-\frac{(a+b)}{2}\right) &\longleftarrow t. \end{align*} \]

Entonces la función transformada de Fourier en el intervalo \([a,b]\) será de la forma: \[ \begin{align*} S_m(t)= & \frac{a_0+a_m\cos\left(\frac{2\pi m}{(b-a)}\left(t-\frac{(a+b)}{2}\right)\right)}{2}\\ & \quad +\sum_{j=1}^{m-1} a_j\cos\left(\frac{2\pi j}{(b-a)}\left(t-\frac{(a+b)}{2}\right)\right)\\ & \quad +\sum_{j=1}^{m-1}b_j\sin\left(\frac{2\pi j}{(b-a)}\left(t-\frac{(a+b)}{2}\right)\right), \end{align*} \]

donde: \[ \begin{align*} a_k = & \frac{1}{m}\sum_{i=0}^{2m-1}y_i\cos\left(\frac{2\pi k}{(b-a)}\left(t_i-\frac{(a+b)}{2}\right)\right),\ k=0,1,\ldots,m,\\ b_k = & \frac{1}{m}\sum_{i=0}^{2m-1}y_i\sin\left(\frac{2\pi k}{(b-a)}\left(t_i-\frac{(a+b)}{2}\right)\right),\ k=0,1,\ldots,m-1.\\ \end{align*} \]

En la práctica se transforma la nube de puntos \((t_k,y_k)\), \(k=0,1,\ldots,2m-1\) en el intervalo \([-\pi,\pi]\), \[ (t_k,y_k)\longrightarrow \left(x_k=\frac{2\pi}{(b-a)}\left(t_k -\frac{(a+b)}{2}\right),y_k\right),\ k=0,1,\ldots,2m-1, \] se calcula la función transformada de Fourier en el intervalo \([-\pi,\pi]\) hallando los coeficientes \(a_k\) y \(b_k\) y por último se considera la expresión anterior \(S_m(t)\).