FFT python scipy

La librería scipy de python incluye algoritmos para la transformada rápida de Fourier FFT.

Como ejemplo, el resultado para

x(t)=\cos(2\pi fs t)

con fs=20 se presenta en un gráfico con:
– la señal en el tiempo,
– el componente real de FFT
– el componente imaginario de FFT
– la magnitud de FFT- la fase de FFT.

Observe las escalas:la frecuencia se muestra en Hz.

Se recomienda probar el script de python con otras funciones.

# Transformadas Rapida de Fourier
# entrada funcion par
# Propuesta: edelros@espol.edu.ec
import numpy as np
import scipy.fftpack as fourier
import matplotlib.pyplot as plt

# Definir la funcion de Entrada par
def entradax(t,fs):
    x=np.cos(2*np.pi*fs*t)
    return(x)

# PROGRAMA
# INGRESO
t0=0
tn=0.5 #float(input('rango segundos [0,tn]:'))
n= 1024 #int(input(' n muestras en rango:'))

# PROCEDIMIENTO
dt=(tn-t0)/n    # intervalo de muestreo
# Analógica Referencia para mostrar que es par
fs=20
t=np.arange(-tn,tn,dt)  # eje tiempo analógica
m=len(t)
xanalog=np.zeros(m, dtype=float)
for i in range(0,m):
    xanalog[i]=entradax(t[i],fs)

# FFT: Transformada Rapida de Fourier
# Analiza la parte t>=0 de xnalog muestras[n:2*n]
xf=fourier.fft(xanalog[n:2*n])
xf=fourier.fftshift(xf)
# Rango de frecuencia para eje
frq=fourier.fftfreq(n, dt)
frq=fourier.fftshift(frq)

# x[w] real
xfreal=(1/n)*np.real(xf)
# x[w] imaginario
xfimag=(1/n)*np.imag(xf)
# x[w] magnitud
xfabs=(1/n)*np.abs(xf)
# x[w] angulo
xfangle=(1/n)*np.unwrap(np.angle(xf))

#SALIDA
plt.figure(1)       # define la grafica
plt.suptitle('Transformada Rápida Fourier FFT')

plt.subplot(321)    # grafica de 3x2, subgrafica 1
plt.ylabel('xanalog[t]')
plt.xlabel('tiempo')
plt.plot(t,xanalog)
plt.margins(0,0.05)
plt.grid()

ventana=0.2 # ventana de frecuencia a observar alrededor f=0
ra=int(len(frq)*(0.5-ventana))
rb=int(len(frq)*(0.5+ventana))
plt.subplot(323)    # grafica de 3x2, subgrafica 3
plt.ylabel('x[f] real')
plt.xlabel(' frecuencia (Hz)')
plt.plot(frq[ra:rb],xfreal[ra:rb])
plt.margins(0,0.05)
plt.grid()

plt.subplot(325)    # grafica de 3x2, subgrafica 5
plt.ylabel('x[f] imag')
plt.xlabel(' frecuencia (Hz)')
plt.plot(frq[ra:rb],xfimag[ra:rb])
plt.margins(0,0.05)
plt.grid()

plt.subplot(324)    # grafica de 3x2, subgrafica 4
plt.ylabel('x[f] magnitud')
plt.xlabel(' frecuencia (Hz)')
plt.plot(frq[ra:rb],xfabs[ra:rb])
plt.margins(0,0.05)
plt.grid()

plt.subplot(326)    # grafica de 3x2, subgrafica 6
plt.ylabel('x[f] fase')
plt.xlabel(' frecuencia (Hz)')
plt.plot(frq[ra:rb],xfangle[ra:rb])
plt.margins(0,0.05)
plt.grid()

plt.show()

Para el cálculo de la transformada se han usado los valores de t≥0

la transformada calcula los valores positivos y luego los negativos, por lo que se usa fftshift(xf) antes de graficar para poner el orden apropiado.

Las k muestras debe ser un valor de 2n tal como 512, 1024, 2048, …

Audio en formato .wav

Lectura de archivo .wav en python

Para procesar un archivo de audio en formato .wav se usarán las instrucciones de la libreria scipy.
La instruccion usa el archivo con ‘nombre.wav’ dado como variable tipo texto y se obtienen dos variables que representan:
fsonido: frecuencia de muestreo del sonido en PCM y
sonido: que es un arreglo con las muestras del sonido.

El archivo de audio debe encontrarse en el mismo directorio que el script de python, por ejemplo, dado el archivo ‘Alarm.wav’ , se puede procesar con las instrucciones mostradas:

Alarm01.wav

import scipy.io.wavfile as waves

# INGRESO
archivo='Alarm01.wav'
fsonido, sonido = waves.read(archivo)

con lo que se obtiene:

>>> fsonido
22050
>>> sonido
array([[0, 0],
       [0, 0],
       [0, 0],
       ..., 
       [0, 0],
       [0, 0],
       [0, 0]], dtype=int16)
>>> np.shape(sonido)
(122868, 2)

En el ejemplo, la frecuencia de muestreo es de 22050 Hz. El sonido es estéreo al tener dos columnas que corresponden a los canales izquierdo y derecho.
Para usar un solo canal, se copian los datos a un nuevo arreglo. Para separar el canal izquierdo por ejemplo, se usan las instrucciones:

>>> izquierdo=sonido[:,0].copy()
>>> izquierdo
array([0, 0, 0, ..., 0, 0, 0], dtype=int16

Con lo que se tienen los datos listos para ser procesados.

Guardar un archivo de audio .wav

Luego de procesar los datos de audio, y guardar el resultado en un archivo con ‘nombre.wav’ se usa la instruccion waves.write() de la libreria scipy, que requiere:
archivo: el nombre del archivo resultante, con extension.wav
fsonido: la frecuencia de muestreo del sonido (entero)
sonido: el arreglo de la señal de audio como entero de 16 bits (dtype=’int16′).

Recuerde haber invocado a las librerias numpy y scipy que se presentan como referencia en el script.

import numpy as np
import scipy.io.wavfile as waves

# PROCEDIMIENTO
# Arreglos para datos con k muestras
sonidofinal=np.zeros(k, dtype='int16')

# Salida
archivo='audiofinal.wav'
waves.write(archivo, int(fsonido),sonidofinal)

El archivo de audio resultante se escucharà usando un programa como «windows media player»

Notas:
Se puede añadir al nombre la ruta de ubicación del archivo en el disco duro. Ejemplo para windows: ‘C:\Users\mis documentos\archivo.wav

Gráfica de un canal de audio

Para mostrar en un gráfico un canal de audio se usa la libreria matplotlib. Para muestra, usando los datos del ejemplo anterior:

import matplotlib.pyplot as plt

izquierdo=sonido[:,0].copy()

# Salida a gràfico
plt.plot(izquierdo)
plt.show()

con lo que se obtiene la siguiente gráfica:

Con lo que se puede revisar la forma de la señal de audio

scipy Resumen

la librería de funciones scipy (scientific python), dispone de funciones para el tratamiento de señales, estadíasticas, audio, entre otras.

El orden de las instrucciones es el que aparece en los post del blog.

instrucciones
import scipy.io.wavfile as waves librerias de audio en formato wac para lectura y escritura de archivos.
fsonido, sonido = waves.read(archivo) lectura de datos de un archivo de audio en formato wav. Se obtiene la frecuencia de muestreo en fsonido y los datos en sonido
import scipy.integrate as integrate importar metodos de integración de scipy
integrate.simps(valrores, ejex) integral de muestras de señal usando el método de Simpson.
. .
. .

numpy Resumen

la librería de funciones numpy facilitan el cálculo numérico, con datos de vectores y matrices en forma de «arrays».

El orden de las instrucciones es el que aparece en los post del blog.

instrucciones
import numpy as np Importar librerias de funciones Numpy, python numérico, usando un álias de dos letras «np»:
vector=np.arange(a,b,dt) crea un vector con valores en el rango [a,b) y espaciados dt.

t=np.arange(0,10,2)>>> t
array([0, 2, 4, 6, 8])
np.pi constante con valor π

>>> np.pi
3.141592653589793
np.sin(t)
np.cos(t)
función trigonométrica en radianes. La variable t puede ser un escalar o un arreglo.

>>> t=0.65
>>> np.sin(0.65)
0.60518640573603955
>>> t=[0, 0.3, 0.6]
>>> np.sin(t)
array([ 0. , 0.29552021, 0.56464247])
np.abs() obtiene el valor absoluto de un número. En el caso de un número complejo obtiene la parte real.
np.real(complejo)
np.imag(complejo)
obtiene la parte real de los números complejos en un vector. Se aplica lo mismo para la parte imaginaria del número complejo.
complex(a,b) crea el número complejo a partir de los valores de a y b.
a=2
b=3
el resultado es: 2+3j
np.piecewise(t, t>=donde, [1,0]) función que crea a partir de t, los valores de la condicion t>=donde, ubicando los valores de 1, para otro caso es 0. Usada en la funcion escalón.
np.roots([a,b,c]) obtiene las raíces del polinomio:
ax2+bx+c

siendo:

x2 + 3 x + 2 = (x+1)(x+2)
>>> np.roots([1,3,2])
array([-2., -1.])
np.linalg.solve(A,B) Resuelve el sistema de ecuaciones dado por una matriz A y un vector B. siendo, por ejemplo:

 0 =  c1 +  c2
-5 = -c1 - 2c2

c1 = -5 
c2 = 5
>>> A = [[ 1, 1],
	 [-1,-2]]
>>> B = [0,-5]
>>> np.linalg.solve(A,B)
array([-5.,  5.])
>>>

matplotlib Resumen

La librería de funciones matplotlib permite realizar las gráficas de datos. En los ejercicios se usan como datos los vectores y matrices «arrays».

El orden de las instrucciones es el que aparece en las entradas del blog.

instrucciones
%matplotlib inline intrucción de IPython para que las gráficas se incluyan en la página.
De no usarla, las gráficas aparecen en ventanas aparte.
import matplotlib.pyplot as plt Importar librerias de funciones matplotlib.pyplot, usando un álias de tres letras «plt»:
plt.plot(x,y) genera un gráfico de línea con los valores de x,y
plt.show() muestra el gráfico creado con las intrucciones.
Es la última instrucción a usar luego de crear el gráfico.
plt.xlabel(‘textoejex’)

plt.ylabel(‘textoejey’)

Asigna nombres a los ejes de abscisas y ordenas. El nombre se escribe entre ‘apóstrofes’ o «comillas».
plt.stem(x,y) gráfico de líneas verticales y un punto. Usado para mostrar señales discretas en los libros de texto de la bibliografía.
plt.figure(k) permite generar varias gráficas, numeradas cada una por el valor de k. En Python simple se muestran en ventanas separadas.
plt.title(‘texto’) escribe el título del gráfico, definido por ‘texto’
plt.fill_between(rangox, 0, valores, color=’green’ dibuja en el rango un área entre 0 y los valores, al color descrito: ‘green’, ‘ligthgreen’, ‘red’, ‘magenta’, ‘blue’, ‘yellow’, etc
plt.axvline(x=0, color=’red’)
plt.axhline(0, color=’red’)
ubica una linea vertical en el valor dado en x. ejemplo x=0.
El otro caso se refiere a una línea horizontal
plt.grid(True)
plt.grid(False)
dibuja líneas de división en la gráfica.
plt.legend() incluye una leyenda para cada curva a graficar. Trabaja en conjunto con
plt.plot(x,y, label=’nombre curva’)
plt.margins(0.1) crea un margen en la gráfica de valor 0.1 expandiendo los bordes sobre los máximos y mínimos. Se usa para ver claramente las variaciones en los extremos o cuando los valores máximos son constantes en un periodo.
. .
. .