In the following Python script it is shown how the aliasing effect (having signals whose frequency is beyond half of the sampling frequency, i.e., \(\frac{f_s}{2}\)) can be intuitively understood as a folding of the frequencies outside the \(\left[0, \frac{f_s}{2}\right]\) range.

In the script, the sampling frequency is 10 kHz. Then, 3 sinusoids of frequency \(f_{s1} = 500~Hz\), \(f_{s2} = 6~kHz\) and \(f_{s3} = 11~kHz\) are generated. The FFT of the sum of s1, s2 and s3 is computed and plotted. Given that the Nyquist frequency in the system is 5 kHz, it can be seen on the FFT plot how the frequencies wrap around \(\frac{f_s}{2}\) and 0.

Just as a reminder, the way the X-axis on the FFT is explained in FFT resolution.

import numpy as np
import numpy.fft as fft
import matplotlib.pyplot as plt

# Define sampling frequency 
fs = 10e3
Ts = 1/fs

# Define time array
n_periods = 1000
total_time = 1/fs * n_periods
t = np.arange(0, total_time, Ts)


# %% Define signal 1
f1 = 500
s1 = np.sin(2*np.pi*f1*t)


# %% Define signal 2
# Since fs = 5_000, f2 will fold to 4_000 Hz bin
f2 = 6_000
s2 = 0.5 * np.sin(2*np.pi*f2*t)

# %% Define signal 3
# Since fs = 5_000, f3 will fold first on 5_000 Hz and then on 0 again. Therefore, f3 will be shown at 1_000 Hz bin
f3 = 11_000
s3 = 0.25 * np.sin(2*np.pi*f3*t)

# %% Define signal as a sum of s1, s2 and s3
s = s1 + s2 + s3

# %% Compute FFT 
# Total number of samples 
N_samples = len(s)
print(f"Number of samples: {N_samples}")
# Floor in case division is not integer
half_N_samples = int(np.floor(N_samples/2))
# Define frequency bin width 
delta_f = fs/N_samples

# Compute frequency bins until Nyquist frequency
freq_bins = np.arange(0, fs/2, delta_f)

Y = fft.fft(s)

plt.figure(figsize=(8,6), dpi=300)

# Plot only half of the FFT result. Plot absolute value and scale it to represent amplitude
plt.plot(freq_bins, 2.0 / N_samples * np.abs(Y)[:half_N_samples])
plt.grid(True, which='both')
plt.minorticks_on()
plt.title('Signal FFT')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Amplitude')
plt.show()

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:

sol-flauta

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

Una manera sencilla de realizar FFT desde cualquier sistema operativo sin necesidad de disponer de programas como MATLAB, Octave o Mathematica, es utilizando un script de Python. Este lenguaje cuenta con la ventaja de ser muy sencillo y además de disponer una ingente cantidad de librerías para poder hacer casi cualquier cosa.
El objetivo de este pequeño programa es realizar una FFT de un trozo de audio grabado directamente desde un micrófono. Para ello se realiza una escucha pasiva a la espera de un sonido cuyo volumen esté por encima de un determinado umbral. Este umbral se puede modificar con la variable threshold, y dependiendo del micrófono y su sensibilidad debe de ser ajustado de manera diferente. Una vez detectado, se graba unos pocos segundos de audio, cuya longitud puede ser modificada en el segundo bucle for. Actualmente graba aproximadamente 1 segundo, pero si se cambia el 7 por un 100 la longitud es de 13 segundos.
El audio se guarda como un archivo para ser posteriormente utilizado en el cálculo de la FFT y ser graficada al momento. Todo este procedimiento está dentro de un bucle infinito, por lo que, aunque la ventana de la gráfica es bloqueante, está continuamente escuchando.

Las librerías necesarias para poder ejecutar este script son: numpy, matplotlib, pyaudio y scipy.

#!/usr/bin/python
# -*- coding: utf-8 -*-
 
import pyaudio
import struct
import math
import wave
import numpy as np
import matplotlib.pyplot as plt
from scipy.io import wavfile
#import urllib  
#import urllib2 
import os
 
class Microphone:
 
    def rms(self,frame):
        count = len(frame)/2
        format = "%dh"%(count)
        shorts = struct.unpack( format, frame )
        sum_squares = 0.0
        for sample in shorts:
            n = sample * (1.0/32768.0)
            sum_squares += n*n
        rms = math.pow(sum_squares/count,0.5);
        return rms * 1000
 
 
    def passiveListen(self,persona):
 
        CHUNK = 1024; RATE = 8000; THRESHOLD = 200; LISTEN_TIME = 5
 
        didDetect = False
 
        # prepare recording stream
        p = pyaudio.PyAudio()
        stream = p.open(format=pyaudio.paInt16, channels=1, rate=RATE, input=True, frames_per_buffer=CHUNK)
 
        # stores the audio data
        all =[]
 
        # starts passive listening for disturbances 
        print RATE / CHUNK * LISTEN_TIME
        for i in range(0, RATE / CHUNK * LISTEN_TIME):
            input = stream.read(CHUNK)
            rms_value = self.rms(input)
            print rms_value
            if (rms_value < THRESHOLD):
                didDetect = True
                print "Listening...\n"
                break
 
        if not didDetect:
            stream.stop_stream()
            stream.close()
            return False
 
 
        # append all the chunks
        all.append(input)
        for i in range(0, 7):
            data = stream.read(CHUNK)
            all.append(data)
 
        # save the audio data   
        data = ''.join(all)
        stream.stop_stream()
        stream.close()
        wf = wave.open('audio.wav', 'wb')
        wf.setnchannels(1)
        wf.setsampwidth(p.get_sample_size(pyaudio.paInt16))
        wf.setframerate(RATE)
        wf.writeframes(data)
        wf.close()
 
        return True
 
if __name__ == '__main__':
    mic = Microphone()
    while True:
        if  mic.passiveListen('ok Google'):
            fs, data = wavfile.read('audio.wav')
            L = len(data)
            c = np.fft.fft(data) # create a list of complex number
            freq = np.fft.fftfreq(L)
            #freq = np.linspace(0, 1/(2L), L/2)
            print freq
            freq_in_hertz = abs(freq * fs)
            plt.plot(freq_in_hertz, abs(c))
            plt.show()