# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: math
#     language: python
#     name: math
# ---

# # Transformada Discreta de Fourier
#
# La Transformada Discreta de Fourier (*Discrete Fourier Transform* - DFT) es una técnica utilizada para convertir datos espaciales o temporales al *dominio de la frecuencia*.
#
# Podemos dividir la transformada de Fourier en diferentes tipos. La subdivisión más básica es la basada en el tipo de datos sobre los que la transformada opera: funciones continuas o discretas. En este laboratorio trabajaremos solo con la transformada Discreta de Fourier (DFT). La librería `scipy` provee el módulo `scipy.fft` para calcular la transformada discreta de Fourier utilozando el algoritmo de la transformada rápida de Fourier (*Fast Fourier Transform* - FFT).
#
# Es usual que los términos "Transformada Discreta" y "Transformada Rápida" se usen de forma intercambiable aunque no son precisamente lo mismo. La transformada rápida de Fourier (FFT) es un algoritmo para calcular la transformad discreta de Fourier (DFT).
#
# En general, necesitamso usar la transformada de Fourier si queremos obtener las frecuencias presentes en un señal. En algunas situaciones puede resultar más conveniente estudiar un problema en el dominio de la frecuencia que en el dominio del tiempo.
#
# > If you want to find the secrets of the universe, think in terms of energy, frequency and vibration. (Nikola Tesla)
#

# ### Analizando señales simples
#
# Empecemos analizando señales con una sola componente de frecuencia, i.e. señales de la forma $sin(2 \pi  f t)$.

# +
import numpy as np
from matplotlib import pyplot as plt

FRECUENCIA_MUESTREO = 100  # Frecuencia de muestreo
DURACION = 5  # Duración de la señal generada

def generar_onda_sinusoidal(freq, freq_muestreo, duracion):
    x = np.linspace(0, duracion, freq_muestreo * duracion, endpoint=False)
    frequencies = x * freq
    # 2pi porque sin recibe radianes
    y = np.sin((2 * np.pi) * frequencies)
    return x, y

# Generamos una señal sinusoidal con una frecuencia de 2 Hz y 
# con duración de 5 segundos.
x, y = generar_onda_sinusoidal(5, FRECUENCIA_MUESTREO, DURACION)
plt.plot(x, y)

plt.xlabel('Tiempo [s]')
plt.ylabel('Amplitud de la Señal')
plt.grid()

plt.show()

# +
from scipy import fftpack

# Obtenemos la transformada de Fourier de la señal y
Y = fftpack.fft(y)
# Obtenemos las frecuencias asociadas, note que esto depende
# de las muestras tomadas.
freqs = fftpack.fftfreq(len(y)) * FRECUENCIA_MUESTREO

fig, ax = plt.subplots()

ax.stem(freqs, np.abs(Y), use_line_collection=True)
ax.set_xlabel('Frequencia en Hertz [Hz]')
ax.set_ylabel('Magnitud en el Dominio de la Frecuencia')
ax.set_xlim(-FRECUENCIA_MUESTREO / 2, FRECUENCIA_MUESTREO / 2)

plt.grid(linestyle = "--")
plt.show()
# -

# Note que en la gráfica anterior se observan picos para las frecuencias $f = 5$ Hz y $f = -5$ Hz. Recuerde que la transformada **continua** de Fourier para $f(t) = \sin (\omega t )$ es $F(\omega) = i \pi [\delta(\omega + \omega_0) - \delta (\omega  - \omega_0)]$.

# +
# Generemos dos señales con frecuencias f1 = 5 Hz y f2 = 13 Hz respectivamente
t, x1 = generar_onda_sinusoidal(5, FRECUENCIA_MUESTREO, DURACION)
t, x2 = generar_onda_sinusoidal(13, FRECUENCIA_MUESTREO, DURACION)

# Generemos una señal a partir de las anteriores
x = x1 + 0.5*x2

fig, ax = plt.subplots()
ax.plot(t, x)
ax.set_xlabel('Tiempo [s]')
ax.set_ylabel('Amplitud de la Señal')
plt.grid()
plt.show();

# +
from scipy import fftpack

# Obtenemos la transformada de Fourier de la señal y
X = fftpack.fft(x)
# Obtenemos las frecuencias asociadas, note que esto depende
# de las muestras tomadas.
freqs = fftpack.fftfreq(len(x)) * FRECUENCIA_MUESTREO

fig, ax = plt.subplots()

ax.stem(freqs, np.abs(X), use_line_collection=True)
ax.set_xlabel('Frequencia en Hertz [Hz]')
ax.set_ylabel('Magnitud en el Dominio de la Frecuencia')

plt.grid(linestyle = "--")
plt.show()
# -

# Nuevamente, note que podemos diferenciar las dos componentes de frecuencia de la señal anterior (5 Hz y 13 Hz). Ademas se observa que la amplitud de la señal de 13 Hz es la mitad de la de 5 Hz.
#
# Estes ejemplos ilustran el uso de la transformada de Fourier Discreta. Es muy fácil interpretar la transformada para señales que tienen comonentes de frecuencia específicos, las señales de interés usualmente tienen muchos componentes de frecuencia (posiblemente infinitos componentes) y es donde la transformada discreta resulta extremadamente útil.
#
# <img src="https://upload.wikimedia.org/wikipedia/commons/5/50/Fourier_transform_time_and_frequency_domains.gif"/>

# #### <ins> Problema 1: </ins>
#
# Identifique las componentes de frecuencia presentes en la siguiente señal que fue muestreada a `FRECUENCIA_MUESTREO = 200 Hz` para $0 \leq t \leq 5$.

# +
import numpy as np
from matplotlib import pyplot as plt

# Señal para el problema 1
x = np.load('senial_problema1.npy')
# -

# Resuelva aquí el problema 1.


# ### DFT y filtrado
#
# Veremos como funciona la transformada de Fourier Discreta (`scipy.fftpack.fft()`, `scipy.fftpack.fftfreq()`) y la su inversa `scipy.fftpack.ifft()`.

# +
# Seed para el generador de numeros aleatorios
np.random.seed(31415)

# Frecuencia de la señal
f = 1.0/5
# Frecuencia de muestreo
fs = 50
t = np.arange(0, 20, 1/fs)
# Generemos una señal
senial = np.sin(2 * np.pi * f * t)
# Modelemos algo de ruido
ruido = 0.5 * np.random.randn(t.size)
# Generemos una señal con ruido
senial_ruido = senial + ruido
plt.plot(t, senial_ruido, label='Original signal')
plt.grid()

# +
# Calculemos la transformada de Fourier de la señal con ruido
senial_fft = fftpack.fft(senial_ruido)

# Espectro de potencia de la transformada
espectro_potencia = np.abs(senial_fft)**2

# Obtenemos las correspondientes frecuencias
freq = fftpack.fftfreq(senial_ruido.size, d=1/fs)

# Grafiquemos el espectro de potencia
plt.figure(figsize=(6, 5))
plt.plot(freq, espectro_potencia)
plt.xlabel('Frecuencia [Hz]')
plt.ylabel('$|F|^2$')

# Quedemos solo con las frecuencias positivas
pos_freq_mask = np.where(freq > 0)
pos_freq = freq[pos_freq_mask]

# Encontremos la frecuencia mas alta
peak_freq = pos_freq[espectro_potencia[pos_freq_mask].argmax()]

print("Frecuencia Pico", peak_freq)

# Veamos en donde está la frecuencia mas alta (Frecuencia pico)
axes = plt.axes([0.55, 0.3, 0.3, 0.5])
plt.title('Frecuencia Pico')
plt.plot(pos_freq[:8], espectro_potencia[:8])
plt.setp(axes, yticks=[])
plt.show()

# +
high_freq_fft = senial_fft.copy()
high_freq_fft[np.abs(freq) > peak_freq] = 0
filtered_sig = fftpack.ifft(high_freq_fft)

plt.figure(figsize=(6, 5))
#plt.plot(t, sig, label='Original signal')
plt.plot(t, filtered_sig, label='Filtered signal')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.grid()
plt.show()


# -

# #### <ins> Problema 2: </ins>
#
# Implementar un filtro de paso bajo, i.e. dada una señal `x` deberá elimar todas las frecuencias superiores a una frecuencia dada `f_c`. El ejemplo anterior puede servirle de ejemplo, deberá elimar (hacer que su amplitud sea `0`) todas las frecuencias no deseadas.

def filtro_paso_bajo(x, fs, f_c):
    """
    :param x: (np.array) es la señal quese desea filtrar
    :param fs: es la frecuencia a la que fue muestreada la señal x
    :param f_c: es la frecuencia de corte
    
    :return: returna la señal filtrada como un arreglo de numpy con la misma dimensión que x
    """
    pass


# +
from matplotlib import rcParams
rcParams['figure.figsize'] = 15, 10
import warnings
warnings.filterwarnings("ignore")

# Señal de prueba pra el filtro
fs = 50
t = np.linspace(0, 5, fs*5, endpoint=False)
all_freq = [2, 5, 15, 20]
x = sum([np.sin(2*np.pi*fi*t) for fi in all_freq])

plt.subplot(2, 2, 1)
# Señal orginal
plt.plot(t, x, label='Señal original')
plt.grid()
plt.legend()

plt.subplot(2, 2, 2)
# Filtrar las señales mayores a 4 Hz
y = filtro_paso_bajo(x, fs, 4)
plt.plot(t, y, linewidth=2, label='$f_c = 4 Hz$', color='red')
plt.ylim(-3, 3)
plt.grid()
plt.legend()

plt.subplot(2, 2, 3)
# Filtrar las señales mayores a 10 Hz
y = filtro_paso_bajo(x, fs, 10)
plt.plot(t, y, linewidth=2, label='$f_c = 10 Hz$', color='red')
plt.ylim(-3, 3)
plt.grid()
plt.legend()

plt.subplot(2, 2, 4)
# Filtrar las señales mayores a 6 Hz
y = filtro_paso_bajo(x, fs, 16)
plt.plot(t, y, linewidth=2, label='$f_c = 16 Hz$', color='red')
plt.ylim(-3, 3)
plt.grid()
plt.legend()

plt.show()
# -

# ### Analizando señales de audio con la DFT
#
# Definimos el espectro de un sonido como la representación de la distribución de energía sonora de dicho sonido en función de la **frecuencia**. Utilizaremos la librería [librosa](https://librosa.org/doc/latest/index.html) para generar espectros.

# Instalar librosa.
# !pip install librosa

# +
from IPython.display import Audio # para reproducir audio en este notebook
import librosa # para leer un archivo de audio
from librosa import display # para hacer plot de una señal de audio

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

# Configuración de las gráficas
from matplotlib import rcParams
rcParams['figure.figsize'] = 16, 6
rcParams['axes.grid'] = True
# -

# Leemos un archivo de audio, librosa.load nos devuelve las muestras y la frecuencia de muestreo
file_path = 'sample_sound.wav'
samples, sampling_rate = librosa.load(file_path, sr = None, mono=True, offset= 0.0, duration = None)

duration = len(samples) / sampling_rate
print("El audio tiene una duración de %.2f y está muestreado a fs = %.2f Hz" % (duration, sampling_rate))

# Podemos reproducir el auido
Audio(file_path)

# Graficamos la señal en el dominio del tiempo
librosa.display.waveplot(y = samples, sr = sampling_rate)
plt.xlabel('Tiempo (segundos)')
plt.ylabel('Amplitud')
plt.show()

# Para observar su espectro de frecuencias calculamos la transformada discreta de Fourier con las muestras de la señal de audio.

# +
N = len(samples) # número de muestras
T = 1/sampling_rate # duración - espaciado de cada muestra
yf = fftpack.fft(samples) # transformamos
xf = fftpack.fftfreq(N, d=T)
plt.plot(xf[:N//2], np.abs(yf[:N//2]), color='green') # solo la primera parte contiene frecuencias positivas

plt.xlabel("Frecuencia (en Hz)")
plt.ylabel("Magnitud")
plt.show()
# -

# ### Espectrograma
#
# El espectrograma es una representación visual de las variaciones de la frecuencia en el eje vertical, y de la intensidad del sonido mediante niveles de colores a lo largo del tiempo que se representa en el eje horizontal. 
#
# **Para la obtención del espectrograma se aplica una transformada de Fourier inicialmente a la señal, mediante el algoritmo de la transformada rápida de Fourier o FF**. 
#
# Dependiendo del tamaño de la ventana que utilizamos para el análisis de Fourier tendremos diferentes niveles de resolución del espectrograma. 
#
#  - Si se aplica una ventana muy grande obtendremos un espectrograma muy detallado pero a costa de incrementar el tiempo de cálculo necesario para esta operación. 
#  - Para el caso de una ventana demasiado pequeña el efecto es el inverso y no seremos capaces de distinguir los diferentes armónicos si están muy juntos en el espectrograma.

# +
_, _, _, im = plt.specgram(samples, Fs=sampling_rate)
plt.xlabel("Tiempo (en segundos)")
plt.ylabel("Frecuencia (en Hz)")
plt.colorbar(im).set_label("Intensidad (en dB)")

plt.show()
# -

# El siguiente es un espectrograma del sonido de un violín. Note como las líneas brillantes en la parte de abajo corresponden a armónicos fundamentales de cada nota y las otras líneas brillantes cercanas son los sobretonos armónicos.

file_path = 'violin.ogg'
samples, sampling_rate = librosa.load(file_path, sr = None, mono=True, offset= 0.0, duration = None)
Audio(file_path)

# +
_, _, _, im = plt.specgram(samples, Fs=sampling_rate)
plt.xlabel("Tiempo (en segundos)")
plt.ylabel("Frecuencia (en Hz)")
plt.colorbar(im).set_label("Intensidad (en dB)")

plt.show()
# -

# #### <ins> Problema 3: </ins>
#
# Grafique el espectrograma del archivo de audio `problema_3.wav` adjunto con los archivos de este laboratorio, i.e. cargue el archivo de audio, identifique la frecuencia de muestreo, la duración y finalmente grafique su espectrograma.



# <hr>
#
# # **Instrucciones Generales**
# 1. Este laboratorio será realizado de manera *individual*.
# 2. **Fecha de entrega**: Lunes 18 de Octubre de 2021. Su solución debe subirla en un archivo ZIP enviado por GES y debe contener el archivo .ipynb con sus respuesta a cada inciso solicitado en cada uno de los *Problemas* indicados en la parte superior.

StarcoderdataJupyterScriptsDedupFiltered