Categorías

## Basis in a vector space

A vector space basis is the skeleton from which a vector space is built. It allows to decompose any signal into a linear combination of simple building blocks, namely, the basis vectors. The Fourier Transform is just a change of basis.

A vector basis is the linear combination of a set of vector that can write any vector of the space.

$w^{(k)} \leftarrow \text{basis}$

The canonical basis in $$\mathbb{R}^2$$ are:

$$e^{(0)} = [1, 0]^T \ \ e^{(1) } = [0,1]^T$$

Nevertheless, there are more basis for $$\mathbb{R}^2$$:

$$e^{(0)} = [1, 0]^T \ \ e^{(1) } = [1,1]^T$$

This former basis is not linearly independent as information of $$e^{(0)}$$ is inside $$e^{(1)}$$.

## Formal definition

H is a vector space.

W is a set of vectors from H such that $$W = \left\{ w^{(k)} \right\}$$

W is a basis of H if:

1. We can write $$\forall x \in H$$: $$x = \sum_{k=0}^{K-1} \alpha_k w^{(k)}, \ \ \alpha_k \in \mathbb{C}$$
2. $$\alpha_k$$ are unique, namely, there is linear independence in the basis, as a given point can only be expressed in a unique combination of the basis.

Orthogonal basis are those which inner product returns 0:

$$\left \langle w^{(k)}, w^{(n)} \right \rangle = 0, \ \ \text{for } k \neq n$$

In addition, if the self inner product of every basis element return 1, the basis are orthonormal.

## How to change the basis?

An element in the vector space can be represented with a new basis computing the projection of the current basis in the new basis. If $$x$$ is a vector element and is represented with the vector basis $$w^{(K)}$$ with the coefficients $$a_k$$, it can also be represented as a linear combination of the basis $$v^{(k)}$$ with the coefficients $$\beta_k$$. In a mathematical notation:

$x = \sum_{k=0}^{K-1} \alpha_k w^{(k)} = \sum_{k=0}^{K-1} \beta_k v^{(k)}$

If $$\left\{ v^{(k)} \right\}$$ is orthonormal, the new coefficients $$\beta_k$$ can be computed as a linear combination of the previous coefficients and the projection of the new basis over the original one:

$\beta_h = \left \langle v^{(h)}, x \right \rangle = \left \langle v^{(h)}, \sum_{k=0}^{K-1} \alpha_k w^{(k)} \right \rangle = \sum_{k=0}^{K-1} \alpha_k \left\langle v^{(h)}, w^{(k)} \right \rangle$

This operation can also be represented in a matrix form as follows:

$\beta_h = \begin{bmatrix} c_{00} & c_{01} & \cdots & c_{0\left(K-1 \right )}\\ & & \vdots & \\ c_{\left(K-1 \right )0} & c_{\left(K-1 \right )01} & \cdots & c_{\left(K-1 \right )\left(K-1 \right )} \end{bmatrix}\begin{bmatrix} \alpha_0 \\ \vdots \\ \alpha_{K-1} \end{bmatrix}$

This operation is widely used in algebra. A well-known example of a change of basis could be the Discrete Fourier Transform (DFT).

Categorías

## Energy and power of a signal

$E_x = \sum_{n=-\infty}^{+\infty} \left| x\left[ n \right] \right|^2$
$P_x = \lim_{N\to\infty} \frac{1}{2N + 1} \sum_{n = -N}^{N} \left| x[n] \right|^2$

Categorías

## Not all sinusoids are periodic in discrete time

For a signal being periodic, it must fulfill the following condition:

$x[n] = x[n+M]$
where $$M \in \mathbb{Z}$$.

If $$x[n] = e^{j\left(\omega n + \phi\right)}$$, then:

$e^{j\left(\omega n + \phi\right)} = e^{j\left(\omega \left(n+M\right) + \phi\right)}$

$e^{j\omega n } e^{j\phi} = e^{j\omega n }e^{j\omega M} e^{j\phi}$

$1 =e^{j\omega M}$

For $$e^{j \omega M}$$ to be 1, the angle must be multiple of $$2\pi$$. Then:

$\omega M = 2\pi N$

where $$N \in \mathbb{Z}$$ and $$2 \pi N$$ represents any integer multiple of $$2\pi$$.
Then, the frequency must meet:
$\omega = \frac{2\pi N}{M}$

This is, $$\omega$$ must be a rational multiple of $$2\pi$$.

What does this mean?

Basically, that a periodic signal can’t take any arbitrary frequency in discrete time.

Categorías

## Onda estacionaria

Onda estacionaria creada por dos ondas viajeras.

Categorías

## Resolución FFT

$\Delta f = \frac{F_s}{N_{samples}}$

Si por ejemplo, $$F_s = 5~GSa/s = 5\cdot10^9~Sa/s$$ y $$N_{samples} = 25000$$, la resolución frecuencial es de:

$\Delta f = \frac{5 \cdot 10^9}{25000} = 200~kHz$

Categorías

## Sistemas temporales discretos

La definición de sistema discreto en el tiempo corresponde a una transformación o operador que relaciona un conjunto de valor de entrada $$x[n]$$ con un conjunto de valores de salida $$y[n]$$.
$y[n] = T\{x[n]\}$

## Sistemas sin memoria

Son aquellos en los que los valores de salida $$y[n]$$ solo dependen de valores presentes de la entrada $$x[n]$$. Es decir, no dependen de $$x[n-1]$$, $$x[n-2]$$, …
Un sistema sin memoria podría ser:
$y[n] = \left( x[n] \right)^2$

## Sistemas lineales

Los sistemas lineales cumplen la propiedad de superposición. Esto significa que la salida será  igual a la suma proporcional de las entradas.

$T\{a\cdot x_1[n] + b\cdot x_2[n] \} = a\cdot T\{x_1[n]\} + b\cdot T\{x_[n]\}$

## Sistemas invariantes en tiempo

Un sistema invariante en el tiempo es aquel cuya salida se ve desplazada en tiempo si así lo hace la entrada. Por tanto, un desplazamiento temporal de la señal de entrada, provocará un desplazamiento en la señal de salida. De esta manera, un sistema será invariante si para todo $$n_0$$, la señal $$x_1 = x[n-n_0]$$ produce una señal a la salida cuyo valor se $$y_1[n] = y[n-n_0]$$.

Un ejemplo de un sistema no invariante en tiempo es:
$y[n] = x[Mn]$
La respuesta $$y_1[n-n_0]$$ a la entrada $$x_1[n-n_0]$$ es:
$y_1[n] = x_1[Mn] = x[Mn-n_0]$
Como vemos, $y[n-n_0] = x[M(n-n_0)] \neq y_1[n]$
Por tanto, el sistema no el invariante.

## Sistemas causales

Un sistema causal es aquel que para determinar los valores de salida del sistema, solo se necesiten los valores presentes o pasados de la entrada.

## Sistemas estables

Son aquellos que para cada una de las posibles entradas finitas del sistema, la salida siempre es finita. Cada una de las posibles entradas finitas del sistema significa:
$|x[n]| \leq B_x < \infty$ Una salida finitia signfica: $|y[n]| \leq B_y < \infty$

## Sistemas lineales invariantes en el tiempo

Los sistemas que son invariantes en tiempo y a la vez lineales se pueden estudiar fácilmente a través del operador convolución. De esta manera, es posible obtener cuál va a ser la salida de un sistema en función de su entrada.

Un sistema queda completamente descrito a través de su función de transferencia $$h[n]$$. Esta función de transferencia es la salida que obtenemos cuando a la entrada aplicamos un impulso $$\delta[n]$$.

## Autofunciones en sistemas lineales invariantes

Como hemos visto, los sistemas lineales invariantes son de especial interés debido a que tienen ventajas sobre los sistemas no lineal o no invariantes y es que se puede describir mediante su respuesta impulsional.

Por propiedades de la función $$\delta [n]$$, es posible describir la secuencia $$x[n]$$ como:
$x[n] = \sum_{k=-\infty}^{+\infty} x[k]\cdot \delta[n-k]$

De este modo, al aplicar la secuencia $$x[k]$$ a la entrada de un sistema, su salida será:
$y[n] = T\left\{x[n] \right\} = T\left\{\sum_{k=-\infty}^{+\infty} x[k]\cdot \delta[n-k] \right\} =\sum_{k=-\infty}^{+\infty} x[k]\cdot T\left\{ \delta[n-k] \right\}$

Si identificamos la transformación del tren de deltas $$T\left\{ \delta[n-k] \right\}$$ como su respuesta impulsional $$h[n]$$, obtenemos lo que se conoce como su ecuación de convolución:
$y[n] = \sum_{k=-\infty}^{+\infty} x[k]\cdot h[n-k]$

Como ya se explicó en otra entrada, las autofunciones son útiles para la caracterización de sistemas. Las autofunciones son aquellas que al aplicarse en la entrada de un sistema, su salida es la misma que la función de la entrada multiplicada por una constante. Es decir, si $$x[n]$$ es una autofunción de $$y[n]$$:

$y[n] = T\left\{x[n] \right\} = a \cdot x[n]$

donde $$a$$ es una constante.

Una de las autofunciones de los sistemas invariantes es la exponencial compleja $$x[n] = z^n$$ donde $$z$$ es un número complejo cualquiera. Esta función es autofunción de cualquier sistema lineal invariante ya que:
$y[n] = \sum_{k=-\infty}^{+\infty} z^{n-k} h[k] = z^n\sum_{k=-\infty}^{+\infty} z^{-k} h[k]$

Si esta serie converge, ya que el sumatorio no depende de $$n$$, $$y[n]$$ puede escribirse como:
$y[n] = H(z) z^n$

Este es lo mismo que decir que si a la entrada de un sistema lineal invariante se aplica una secuencia exponencial $$z^n$$, a la salida se obtiene la misma secuencia multiplicada por una constante. Por tanto, queda demostrado que las exponenciales complejas $$z^n$$ son autofunciones de un sistema lineal invariante cuyos autovalores están determiandos por:
$H(z) = \sum_{k=-\infty}^{+\infty} z^{-k} h[k]$

Si se interpreta $$z$$ como una variable, $$H(z)$$ es una función compleja de variable independiente compleja que se conoce como función de transferencia y describe el comportamiento del sistema frente a cualquier entrada.

## Sistemas definidos por ecuaciones en diferencias finitas

Los sistemas discretos se pueden caracterizar, además de por su respuesta impulsional, por ecuaciones en diferencias finitas lineales con coeficientes constantes. Una ecuación en diferencias finitas lineal con coeficientes constantes corresponde a:
$\sum_{k=0}^P a_k y[n-k] = \sum_{k=0}^{Q} b_k x[n-k]$
donde $$x[n]$$ es una dato (ya que es la secuencia de entrada), $$y[n]$$ es la incógnita de la ecuación y $$a_k$$ y $$b_k$$ son los coeficientes independientes de $$n$$. El orden de la ecuación corresponde al mayor entre $$P$$ y $$Q$$.

No obstante, si se describe un sistema mediante ecuaciones en diferencias finitas, este no queda totalmente caracterizado debido a la falta de información sobre el estado del sistema antes de aplicar ninguna secuencia a la entrada. Es por eso, que para poder determinar cuál será la salida del sistema a una determinada secuencia, es necesario conocer también las condiciones iniciales del sistema. Por ejemplo, en el sistema que define un sumador:
$y[n] = x[n] + y[n-1]$
si se le aplica una secuencia constante $$x[n] = 1$$ para $$n\geq 0$$, el valor de $$y[0]$$ será:
$y[0] = x[0] + y[-1] = 1 + y[-1]$
Solo se podrá determinar el valor de $$y[0]$$ si conocemos $$y[-1]$$. Normalmente, el sistema se supone que estaba en reposo y que las condiciones iniciales son $$y[n] = 0$$ para $$n<0$$. De esta manera, $$y[0] = 1$$
\subsubsection{Sistemas recurrentes y no recurrentes}

Los sistemas recurrentes son aquellos que proporcionan valores de $$y[n]$$ en función de valores de la propia salida calculados anteriormente. Es decir:
$y[n] = \frac{1}{a_0}\left( \sum_{k=0}^{Q} b_k x[n-k] – \sum_{k=1}^{P} a_k y[n-k] \right)$

Si la salida $$y[n]$$ solo depende de los valores presentes y pasados de $$x[n]$$, o lo que es lo mismo, P = 0, el sistema es no recurrente.

$y[n] = \sum_{k=0}^{Q} \frac{b_k}{a_0} x[n-k]$

### Respuesta impulsional

Para obtener la respuesta impulsional de un sistema lineal invariante definido por diferencias finitas, basta con aplicar a la entrada la secuencia exponencial $$x[n] = z^n$$. Como hemos visto, esta secuencia es autofunción del sistema y por tanto la salida será:
$y[n] = H(z) z^n$

$$H(z)$$ no tiene por qué existir siempre ya que el sumatorio que lo define debe converger.

Un sistema lineal invariante definido por ecuaciones de diferencias finitas, como hemos visto, corresponde a la ecuación:

$\sum_{k=0}^P a_k y[n-k] = \sum_{k=0}^{Q} b_k x[n-k]$

Al sustituir $$x[n]$$ e $$y[n]$$ en la ecuación, se obtiene:
$\sum_{k=0}^P a_k H(z)z^{n-k} = \sum_{k=0}^{Q} b_k z^{n-k}$

$$H(z)$$ no depende de k, por lo que puede salir fuera del sumatorio.

$H(z) \sum_{k=0}^P a_k z^{n-k} = \sum_{k=0}^{Q} b_k z^{n-k}$

Si despejamos $$H(z)$$ obtenemos:

$H(z) = \frac{\sum_{k=0}^{Q} b_k z^{n-k}}{\sum_{k=0}^P a_k z^{n-k}} = \frac{\sum_{k=0}^{Q} b_k z^{-k}}{\sum_{k=0}^P a_k z^{-k}}$

Por tanto, la función de transferencia de un sistema caracterizado por una ecuación en diferencias finitas es un cociente de polinomios en $$z^{-1}$$ cuyos coeficientes son directamente los coeficientes de la ecuación.

Por ejemplo, si queremos diseñar un sistema de orden 2 no recurrente que cancele un tono a la frecuencia $$\omega_0 = \frac{\pi}{4}$$, la ecuación en diferencias finitas debe seguir la siguiente expresión:
$y[n] = \sum_{k=0}^Q b_k x[n-k]$

ya que el sistema es no recurrente. Por ser de orden 2, $$Q=2$$. Por tanto, queda calcular los coeficiente $$b_k$$.

$x[n] = \cos{\left(\frac{\pi}{4}n\right)} = \frac{1}{2} \left( e^{j \frac{\pi}{4} n} + e^{-j \frac{\pi}{4} n} \right)$

Debido a la fórmula de Euler ($$e^{jx} = \cos{x} + j \sin{x}$$), la entrada puede escribirse como combinación lineal de dos exponenciales. Como se ha demostrado, este tipo de funciones son autofunciones del sistema lineal invariante y por tanto la salida será:
$y[n] = \frac{1}{2} e^{j \frac{\pi}{4} n} H\left(e^{j \frac{\pi}{4}} \right) + \frac{1}{2} e^{-j \frac{\pi}{4} n} H\left(e^{-j \frac{\pi}{4}} \right)$

Para cancelar un tono a esta frecuencia la salida debe ser nula, $$y[n] = 0$$. Esto ocurre si:
$H\left(e^{j \frac{\pi}{4}} \right) = H\left(e^{-j \frac{\pi}{4}} \right) = 0$

La función de transferencia $$H(z)$$ será:
$H(z) = \sum_{k=0}^{Q=2} b_k z^{-k}$
ya que $$P = 0$$ por ser un sistema no recurrente.

Podemos desarrollar $$H(z)$$ como:
$H(z) = b_0 + b_1 z^{-1} + b_2 z^{-2}$

Esta expresión se puede escribir como un producto de raíces y es equivalente a:
$H(z) = b_0 \prod_{i=0}^1 \left(1-z_i z^{-1}\right) = b_0 \left( 1 -z_0 z^{-1}\right) \left(1-z_1 z^{-1} \right)$

donde $$z_0 = e^{j \frac{\pi}{4}}$$ y $$z_1 = e^{-j \frac{\pi}{4}}$$.

Podemos identificar ambas expresiones si desarrollamos este producto:
$H(z) = b_0 \left(1 + \left(-z_0-z_1\right) z^{-1} + z_0 z_1 z^{-2} \right)$

El término $$b_0$$ solo determina la ganancia del sistema y puede valer simplemente 1. Por tanto:
$b_0 = b_0 = 1$
$b_1 = -z_0 -z_1 = -e^{j \frac{\pi}{4}}-e^{-j \frac{\pi}{4}} = 2\cdot \cos{\frac{\pi}{4}} = \sqrt{2}$
$b_2 = z_0 z_1 = e^{j \frac{\pi}{4}}\cdot e^{-j \frac{\pi}{4}} = 1$

Finalmente, $$H(z)$$ queda como:
$H(z) = 1 + \sqrt{2}z^{-1} + z^{-2}$

y:

$y[n] = 1 + \sqrt{2} z^{-1} + z^{-2}$

Categorías

Categorías

## Vector propio y valor propio

Los vectores propios de un operador lineal son aquellos que después haber sido transformados por el operador, tienen la misma dirección que tenían escalados un cierto valor. Este valor de reescalado es el valor propio.

By – The image of the Mona Lisa used in this image is Image:Mona_Lisa.jpeg. Arrows created using w:Inkscape, transformations made using w:GIMP, Public Domain, Link

En la imagen tenemos dos vectores, el rojo y el azul representados sobre el cuadro. Si realizamos la transformación de rotar el cuadro, es decir después de aplicar el operador lineal de rotación, obtenemos que el vector rojo sigue teniendo la misma dirección que tenía con la misma amplitud que tenía. Por tanto, el vector rojo es un vector propio (eigenvector) de esta transformación. Sin embargo, el vector azul sí cambia su dirección y no es vector propio.

Los eigenvectores tienen especial relevancia en la teoría de sistemas, ya que permiten escribir una convolución como una multiplicación.

Si la entrada a un sistema H es $$x(t)$$ y lo podemos definir como una suma ponderada de sinusoides complejas:

$x(t) = \sum_{k=1}^M{a_k e^{j2\pi f_k t}}$

Si $$e^{j2\pi f_k t}$$ es un vector propio del sistema con un valor propio igual a $$H(j2\pi f_k)$$, en lugar de escribir la salida como:

$y(t) = h(t) \ast x(t)$

Podemos escribir la salida del sistema como:

$y(t) = \sum_{k=1}^M{a_k H\left( j 2 \pi f_k\right) e^{j 2 \pi f_k t}}$

Por tanto, podemos describir el comportamiento de una señal en función de su frecuencia en lugar del tiempo. Y este es el principio del análisis de Fourier de las señales.

Categorías

## Series de Fourier

Una manera de representar una señal periódica, es mediante su serie de Fourier. La serie de Fourier de una función tiene la ventaja de expresar la misma función como una suma de senos y cosenos. Por ejemplo, una señal como un tren de pulsos, que tiene una expresión analítica poco rigurosa matematicamente hablando, se puede descomponer como suma senos y cosenos o lo que es lo mismo, sumas de exponenciales complejas.

Si tenemos una señal v(t) con un periodo $$T_0 = 1/f_0$$ , su desarrollo en serie de Fourier exponencial es
$v(t) = \sum_{n = -\infty}^{+\infty} c_n e^{j2\pi n f_0 t} ~~~~~ n = 0, 1, 2, …~~~(1)$

Los coeficientes de la serie se pueden calcular como
$$c_n = \int_{T_o} v(t) e^{-j2 \pi n f_0 t}dt$$

Una manera intuitiva de interpretar esta expresión, es pensar en el término $$c_n$$ como el peso que tiene la frecuencia $$n\cdot f$$ dentro de la señal v(t). De esta manera se promedia el producto $$v(t) e^{-j2 \pi n f_0 t}$$ para todas las frecuencias.

Cabe destacar que $$c_n$$ puede ser un (y normalmente son) números complejos, se pueden expresar en forma polar como
$$c_n = \left| c_n \right| e^{j~arg~c_n}$$

Cabe destacar que:

Todas las frecuencias son multiplos enteros o armónicos de la frecuencia fundamental $$f_0 = 1/T_0$$, por lo que la distancia entre líneas espectrales es uniforme.
La componente continua es el valor promedio de la función.
Si una función es real en tiempo, el espectro frecuencial tendrá simetría par para la amplitud e impar para la fase.

Para poner un ejemplo práctico, vamos ha calcular la serie de Fourier para un tren de pulsos rectangulares con periodicidad $$T_0$$ cuya expresión analítica es:
$v(t)= \left\{\begin{matrix} A & \left|t \right| < \tau/2\\ 0 & \left|t \right| > \tau/2 \end{matrix}\right.$

Aplicando la definición:
$c_n = \frac{1}{T_0} \int_{-T_0/2}^{T_0/2} v(t) e^{-j2 \pi n f_0 t}dt = \frac{1}{T_0} \int_{-\tau/2}^{\tau/2} A e^{-j2 \pi n f_0 t}dt =\\ = \frac{A}{-j2\pi n f_0 T_0} \left( e^{-j \pi n f_0 \tau} – e^{+j \pi n f_0 \tau}\right) = \frac{A}{T_0} \frac{\sin{\left( \pi n f_0 \tau \right)}}{n f_0 \tau} =\\ = \frac{A}{T_0} sinc \left(n f_0 \tau\right)$

Como hemos visto, si manipulamos la expresión de la serie de Fourier podemos llegar a una expresión compacta en la que aparece la función sinc.

Por tanto, si ahora fueramos capaces de generar los infinitos coeficientes $$c_n$$, sustituyendolos en (1) seriamos capaces de reconstruir el tren de pulsos cuadrados a la perfección. Si por el contrario, solo calculamos un número finito de ellos, la función resultante tendría un aspecto parecido al de la imagen, donde se aprecia el fenómeno de Gibbs:

La animación anterior está implementada con la librería JavaScript P5.js. El código fuente es el siguiente:

<code>
let time = 0;
let path = [];
let slider;

function setup(){
createCanvas(600, 400);
slider = createSlider(0, 25, 3, 1);
slider.position(20, 380);
}

function draw(){
background(124,124,124);
textSize(20);
fill(255);
text(' n = ' + slider.value(), width/2, 30);
translate(200, 200);

stroke(255);

let x = 0;
let prevx = 0;
let y = 0;
let prevy = 0;
let c = color(255);

for(let i = 0; i < slider.value(); i++){
fill(c);
stroke(255);
let n = 2*i + 1;
let amp = 100 * 4/(n * TWO_PI);
x += amp * cos(n * time);
y += amp * sin(n * time);
// Point
ellipse(x, y, 5);
line(prevx, prevy, x, y);
noFill();
// Circle
ellipse(prevx, prevy, amp * 2);
prevx = x;
prevy = y;
}

path.unshift(y);

translate(200, 0);

beginShape();
for(let i = 0; i < path.length; i++){
vertex(i, path[i]);
}
endShape();

line(x-200, y, 0, path[0]);

if(path.length > 200){
path.pop();
}

time += 0.05;
}
</code>


Fuente: elaboración propia

Categorías

## Representar el espectro frecuencial de un archivo audio con MATLAB

Puede ser bastante interesante saber cómo poder poder representar gráficamente el espectro de frecuencia un archivo de audio con MATLAB.

Para transformar una función del dominio temporal al frecuencial, es necesario hacer la transformada de Fourier de dicha función. En el caso de una señal contínua la transformada recibe el nombre de Transformada de Fourier mientras que si es de caracter discreto o digital, se le conoce como DFT (Discrete Fourier Transform).

En MATLAB, calcular la transformada de Fourier es tan sencillo como utilizar la funcion FFT (Fast Fourier Transform). Sin embargo, quizá no sea del todo claro cómo se ha de escoger el eje de coordenadas frecuenciales (eje X) para obtener una representación adecuada.

audio = audioread('audio.wav');
Fs = 16000;
L = length(audio);
NFFT = 2^nextpow2(L);
Y = fft(audio, NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
plot(f, 2*abs(Y(1:NFFT/2+1);

Comentemos línea a línea su significado.

audio = audioread('audio.wav');

Cargamos el archivo de audio en una variable llamada ‘audio’ de MATLAB.

Fs = 16000;

Introducimos la frecuencia de muestreo a la cual ha sido adquirido el audio. Es de vital importancia saber este parámetro, ya que de lo contrario no podremos identificar las frecuencias. Si por ejemplo has grabado el audio con un móvil, muchas aplicaciones permiten seleccionar la frecuencia de muestreo. Típicamente la frecuencia de muestro en un móvil es de 16 kHz, aunque para aseguraros podéis abrir el archivo con Audacity y os indicará la frecuencia de muestreo.

L = length(audio);

Calculamos el número de muestras que contiene el audio y lo guardamos en la variable L.

NFFT = 2^nextpow2(L);

La función fft() es más eficiente cuando calcula la transformada de Fourier de una función cuya longitud es una potencia de 2 (2, 4, 8, …). Por este motivo, vamos a calcular cuál es la potencia de 2 más próxima al número de muestras que contiene el audio. Si el archivo tiene 6 muestras, el número que le pasaremos a la función fft será 8, para que sea más eficiente. A este número lo llamaremos NFFT. Hay que tener en cuenta que este número será mayor que el número de muestras que contiene el archivo. Sin embargo, esto no afectará al resultado ya que las muestras que faltan son interpretadas como 0. De esta manera no añaden ningún contenido frecuencial al archivo original.

Y = fft(audio, NFFT)/L;

Calculamos la FFT con la función fft() cuyos argumentos son el audio que queremos transformar (audio), y el número de muestras sobre la que queremos que haga la transformada (NFFT). Dividimos entre L el resultado para normalizar los valores. La transformada la tendremos en la variable Y.

f = Fs/2*linspace(0,1,NFFT/2+1);

Ahora toca construir el eje de coordenadas. Las frecuencia que devuelve la FFt van de la frecuencia 0 hasta la frecuencia mitad de la frecuencia de muestreo $$\frac{F_s}{2}$$ . Por tanto, necesitamos generar un vector que vaya desde 0 hasta $$\frac{F_s}{2}$$. El número de elementos que este vector debe de contener es $$\frac{NFFT}{2}+1$$ ya que debido a que la DFT genera una imagen espejo del espectro de frecuencia, solo querremos visualizar la mitad.

plot(f, 2*abs(Y(1:NFFT/2+1)));

Por último generamos la gráfica cuyo eje de coordenadas está contenido en la variable f y cuyo eje de ordenadas es 2 veces el valor absoluto (para mostrar solo la magnitud y no la parte real e imaginaria) de la mitad del espectro.

También podéis descargar el script audio_fft.m que he creado con función llamada audio_fft(audio, Fs) donde todo esto está ya encapsulado y que podéis descargar pulsando el enlace. Un ejemplo de uso sería:

audio_fft('archivo.wav', 16000);

Como ejemplo he grabado un Sol con una flauta celta y he calculado el espectro frecuencial. El sonido es el siguiente:

Y el resultado es el siguiente:

Donde se puede apreciar claramente la frecuencia fundamental y sus armónicos.