commit 2d0e20368adafedb0ce9a7df5ad56e1fda5e61d5 parent a39703bb68ff3523d9ea6cec3a5e8a762d1b34ac Author: Martin Kloeckner <mjkloeckner@gmail.com> Date: Thu, 27 Nov 2025 12:53:21 -0300 split `main.py` into multiple files Diffstat:
78 files changed, 902 insertions(+), 864 deletions(-) diff --git a/tp/Makefile b/tp/Makefile @@ -14,4 +14,4 @@ open: main nohup zathura $(DOCNAME).pdf > /dev/null 2>&1 & clean: - rm -r *.blg *.bbl *.aux *.log __pycache__ + rm -ff *.blg *.bbl *.aux *.lof *.out *.log diff --git a/tp/main.pdf b/tp/main.pdf Binary files differ. diff --git a/tp/main.py b/tp/main.py @@ -1,431 +0,0 @@ -import numpy as np -from scipy.io import wavfile -from scipy.signal import firwin, freqz, tf2zpk, get_window -from utils import * -import librosa - -## Datos -file1_path = 'data/cancion1.wav' -file2_path = 'data/cancion2.wav' -filter1_h_file_path = 'data/respuesta_impulso_1.txt' -filter2_h_file_path = 'data/respuesta_impulso_2.txt' - -a4_flauta_file_path = 'data/a4_flauta.wav' -a4_clarinete_file_path = 'data/a4_clarinete.wav' -a4_violin_file_path = 'data/a4_violin.wav' - -file1_fs, file1_data = wavfile.read(file1_path) -file2_fs, file2_data = wavfile.read(file2_path) -filter1_h = np.loadtxt(filter1_h_file_path) -filter2_h = np.loadtxt(filter2_h_file_path) - -a4_flauta_fs, a4_flauta_data = wavfile.read(a4_flauta_file_path) -a4_clarinete_fs, a4_clarinete_data = wavfile.read(a4_clarinete_file_path) -a4_violin_fs, a4_violin_data = wavfile.read(a4_violin_file_path) - -file1_filter1_output = np.convolve(file1_data, filter1_h, mode='same') -file1_filter2_output = np.convolve(file1_data, filter2_h, mode='same') -file2_filter1_output = np.convolve(file2_data, filter1_h, mode='same') -file2_filter2_output = np.convolve(file2_data, filter2_h, mode='same') - -############################## Primera Parte ################################### - -def time_domain_cancion1(): - ### grafico completo - time_plot(file1_fs, file1_data, file1_path) - - ### porciones cuasi-periodicas 'cancion1' - time_plot(file1_fs, file1_data, "cancion1_0_248s_a_0_256s", - t=0.248, dt=0.008, a=0.24978, da=0.003) - - time_plot(file1_fs, file1_data, "cancion1_0_520s_a_0_528s", - t=0.520, dt=0.008, a=0.5208, da=0.003) - - ### salida de filtro 'cancion1' - save_to_wav(file1_fs, file1_filter1_output, "file1_filter1_output.wav") - save_to_wav(file1_fs, file1_filter2_output, "file1_filter2_output.wav") - - ### grafico comparando la muestra 1 original y filtrada 1 - data_arr = [normalize(file1_data), normalize(file1_filter1_output)] - leg_arr = ['Señal de audio', 'Señal de audio filtrada'] - - time_plot_multiple(file1_fs, data_arr, leg_arr, - "cancion1_filter1_output_compare") - - time_plot_multiple(file1_fs, data_arr, leg_arr, - "cancion1_filter1_output_compare_0_248_a_0_256", - t=0.248, dt=0.008) - - ### grafico comparando la muestra 1 original y filtrada 2 - data_arr = [normalize(file1_data), normalize(file1_filter2_output)] - leg_arr = ['Señal de audio', 'Señal de audio filtrada'] - - time_plot_multiple(file1_fs, data_arr, leg_arr, - "cancion1_filter2_output_compare") - time_plot_multiple(file1_fs, data_arr, leg_arr, - "cancion1_filter2_output_compare_0_248_a_0_256", - t=0.248, dt=0.008) - -def time_domain_cancion2(): - ### grafico completo - time_plot(file2_fs, file2_data, "cancion2", t=6) - - ### porciones cuasi-periodicas 'cancion2' - time_plot(file2_fs, file2_data, "cancion2_14_72s_a_14_73s", t=14.720, dt=0.01) - time_plot(file2_fs, file2_data, "cancion2_26_57s_a_26_58s", t=26.570, dt=0.01) - - save_to_wav(file2_fs, file2_filter1_output, "file2_filter1_output.wav") - save_to_wav(file2_fs, file2_filter2_output, "file2_filter2_output.wav") - - ### grafico comparando la muestra 2 original y filtrada 2 - data_arr = [normalize(file2_data), normalize(file2_filter1_output)] - leg_arr = ['Señal original', 'Señal filtrada'] - - time_plot_multiple(file2_fs, data_arr, leg_arr, - "cancion2_6s_filter1_output_compare", t=6) - time_plot_multiple(file2_fs, data_arr, leg_arr, - "cancion2_6s_filter1_output_compare_26_57_a_26_58", - t=26.57, dt=0.01) - - ### grafico comparando la muestra 2 original y filtrada 2 - data_arr = [normalize(file2_data), normalize(file2_filter2_output)] - leg_arr = ['Señal original', 'Señal filtrada'] - - time_plot_multiple(file2_fs, data_arr, leg_arr, - "cancion2_6s_filter2_output_compare", t=6) - time_plot_multiple(file2_fs, data_arr, leg_arr, - "cancion2_6s_filter2_output_compare_26_57_a_26_58", - t=26.57, dt=0.01) - -def time_domain_music_instruments(): - ### grafico de los instrumentos musicales - time_plot(a4_flauta_fs, a4_flauta_data, "a4_flauta", t=0.25, dt=0.010) - time_plot(a4_clarinete_fs, a4_clarinete_data, "a4_clarinete", t=0.25, dt=0.010) - time_plot(a4_violin_fs, a4_violin_data, "a4_violin", t=0.25, dt=0.010) - - -############################## Segunda parte ################################## - -def freq_domain_cancion1(): - freq_plot(file1_fs, file1_data, "cancion1_fft", f_max=8000) - freq_plot(file1_fs, file1_filter1_output, "cancion1_filter1_output_fft", - f_max=8000) - freq_plot(file1_fs, file1_filter2_output, "cancion1_filter2_output_fft", - f_max=8000) - -def freq_domain_cancion2(): - freq_plot(file2_fs, file2_data, "cancion2_fft", - f_max=8000) - freq_plot(file2_fs, file2_filter1_output, "cancion2_filter1_output_fft", - f_max=8000) - freq_plot(file2_fs, file2_filter2_output, "cancion2_filter2_output_fft", - f_max=8000) - -def freq_domain_spectograms(): - # Formas de funcion ventana (en tiempo, en frecuencia tienen otra forma) - # 'boxcar': rectangular - # 'bartlett': triangular - # 'hann': similar a medio ciclo de seno - # https://en.wikipedia.org/wiki/Window_function - - for i in [512, 1024, 2048]: - for window in ['boxcar', 'bartlett', 'hann', 'hamming']: - spectogram_plot(file1_fs, file1_data, - f"cancion1_espectograma_{window}_{i:04d}", N=i, - win=window, ylim=[0, 17500]) - - spectogram_plot(file2_fs, file2_data, - f"cancion2_espectograma_{window}_{i:04d}", t=6, N=i, - win=window, ylim=[0, 8000]) - -def a4_flauta_cutoff(): - a4_flauta_cutoff_fft, a4_flauta_cutoff_freqs = freq_compute_fft( - a4_flauta_fs, a4_flauta_data) - - # armonicos mayores a 'cutoff_freq' Hz son descartadas, analogo a aplicar un - # filtro pasabajos ideal - - cutoff_freq = 1000 - a4_flauta_cutoff_fft[np.abs(a4_flauta_cutoff_freqs) > cutoff_freq] = 0.0 - - # espectro de la señal filtrada - N = len(a4_flauta_cutoff_fft) - x = a4_flauta_cutoff_freqs[:N // 2] - y = np.abs(a4_flauta_cutoff_fft[:N // 2]) - fig, axis = freq_graph_data(x, y, f_max=2000, show=False) - save_plot(fig, f"a4_flauta_cutoff_{cutoff_freq}Hz_fft") - - # señal temporal reconstruida - a4_flauta_cutoff = ifft(a4_flauta_cutoff_fft).real - - time_plot(a4_flauta_fs, a4_flauta_cutoff, - f"a4_flauta_cutoff_{cutoff_freq}Hz", 0.25, 0.010) - - save_to_wav(a4_flauta_fs, a4_flauta_cutoff, - f"a4_flauta_cutoff_{cutoff_freq}Hz.wav") - - data_arr = [normalize(a4_flauta_data), normalize(a4_flauta_cutoff)] - leg_arr = [ - 'Señal de nota A4 de flauta', - 'Señal de nota A4 de flauta filtrada' - ] - time_plot_multiple(a4_flauta_fs, data_arr, leg_arr, - "a4_flauta_cutoff_time_comparison", t=0.25, dt=0.008) - - data_arr = [a4_flauta_cutoff, a4_flauta_data] - leg_arr = ["Nota musical A4 con flauta filtrada", "Nota musical A4 con flauta"] - freq_plot_multiple(a4_flauta_fs, data_arr, leg_arr, f_max=4000, t=0.253, - dt=8*0.002272727, save_name="a4_flauta_comparison", show=False) - -def a4_clarinete_cutoff(): - a4_clarinete_cutoff_fft, a4_clarinete_cutoff_freqs = freq_compute_fft( - a4_clarinete_fs, a4_clarinete_data) - - cutoff_freq = 3000 - a4_clarinete_cutoff_fft[np.abs(a4_clarinete_cutoff_freqs) > cutoff_freq] = 0.0 - - # espectro de la señal filtrada - N = len(a4_clarinete_cutoff_fft) - x = a4_clarinete_cutoff_freqs[:N // 2] - y = np.abs(a4_clarinete_cutoff_fft[:N // 2]) - fig, axis = freq_graph_data(x, y, f_max=3000, show=False) - save_plot(fig, f"a4_clarinete_cutoff_{cutoff_freq}Hz_fft") - - # señal temporal reconstruida - a4_clarinete_cutoff = ifft(a4_clarinete_cutoff_fft).real - time_plot(a4_clarinete_fs, a4_clarinete_cutoff, - f"a4_clarinete_cutoff_{cutoff_freq}Hz", 0.25, 0.010) - save_to_wav(a4_clarinete_fs, a4_clarinete_cutoff, - f"a4_clarinete_cutoff_{cutoff_freq}Hz.wav") - - data_arr = [normalize(a4_clarinete_data), normalize(a4_clarinete_cutoff)] - leg_arr = [ - 'Señal de nota A4 de clarinete', - 'Señal de nota A4 de clarinete filtrada' - ] - time_plot_multiple(a4_clarinete_fs, data_arr, leg_arr, - "a4_clarinete_cutoff_time_comparison", t=0.25, dt=0.008) - - data_arr = [a4_clarinete_cutoff, a4_clarinete_data] - leg_arr = [ - "Nota musical A4 con clarinete filtrada", - "Nota musical A4 con clarinete" - ] - - freq_plot_multiple(a4_clarinete_fs, data_arr, leg_arr, f_max=8000, t=0.253, - dt=8*0.002272727, save_name="a4_clarinete_comparison", show=False) - -def a4_violin_cutoff(): - a4_violin_cutoff_fft, a4_violin_cutoff_freqs = freq_compute_fft( - a4_violin_fs, a4_violin_data) - - cutoff_freq = 4000 - a4_violin_cutoff_fft[np.abs(a4_violin_cutoff_freqs) > cutoff_freq] = 0.0 - - # espectro - N = len(a4_violin_cutoff_fft) - x = a4_violin_cutoff_freqs[:N // 2] - y = np.abs(a4_violin_cutoff_fft[:N // 2]) - fig, axis = freq_graph_data(x, y, f_max=4000, show=False) - save_plot(fig, f"a4_violin_cutoff_{cutoff_freq}Hz_fft") - - # señal temporal reconstruida - a4_violin_cutoff = ifft(a4_violin_cutoff_fft).real - time_plot(a4_violin_fs, a4_violin_cutoff, - f"a4_violin_cutoff_{cutoff_freq}Hz", 0.25, 0.010) - save_to_wav(a4_violin_fs, a4_violin_cutoff, f"a4_violin_cutoff_{cutoff_freq}Hz.wav") - - data_arr = [normalize(a4_violin_data), normalize(a4_violin_cutoff)] - leg_arr = ['Señal de nota A4 de violin', 'Señal de nota A4 de violin filtrada'] - time_plot_multiple(a4_violin_fs, data_arr, leg_arr, - "a4_violin_cutoff_time_comparison", t=0.25, dt=0.008) - - data_arr = [a4_violin_cutoff, a4_violin_data] - leg_arr = ["Nota musical A4 con violin filtrada", "Nota musical A4 con violin"] - freq_plot_multiple(a4_violin_fs, data_arr, leg_arr, f_max=8000, t=0.253, - dt=8*0.002272727, save_name="a4_violin_comparison", show=False) - -def a4_flauta_fseries(): - fft_freqs_arr = [] - - for i in [8, 4, 1]: - fft, freqs = freq_compute_fft(a4_flauta_fs, a4_flauta_data, t=0.253, dt=i*0.002272727) - fft_freqs_arr.append([fft, freqs, - f"Serie de Fourier {i} periodo{'s' if i != 1 else ''}"]) - - fig, axis = freq_graph_multiple_data(fft_freqs_arr, show=False, f_max=4500) - save_plot(fig, "a4_flauta_fseries_comparison") - -def a4_clarinete_fseries(): - fft_freqs_arr = [] - - for i in [8, 4, 1]: - fft, freqs = freq_compute_fft( - a4_clarinete_fs, a4_clarinete_data, t=0.253, dt=i*0.002272727) - - fft_freqs_arr.append([fft, freqs, - f"Serie de Fourier {i} periodo{'s' if i != 1 else ''}"]) - - fig, axis = freq_graph_multiple_data(fft_freqs_arr, show=False, f_max=8000) - save_plot(fig, "a4_clarinete_fseries_comparison") - -def a4_violin_fseries(): - fft_freqs_arr = [] - - for i in [8, 4, 1]: - fft, freqs = freq_compute_fft(a4_violin_fs, a4_violin_data, - t=0.253, dt=i*0.002272727) - fft_freqs_arr.append([fft, freqs, - f"Serie de Fourier {i} periodo{'s' if i != 1 else ''}"]) - - fig, axis = freq_graph_multiple_data(fft_freqs_arr, show=False, f_max=8000) - save_plot(fig, "a4_violin_fseries_comparison") - - -############################# Tercera parte ################################### - -cutoff = 2650 # frecuencia de corte -M = 700 # orden FIR (número de coeficientes) - -def filtro_fir(): - # Diseño FIR pasabajos con ventana - b = firwin(M, cutoff, fs=fs, window='hamming') - - w, H = freqz(b, worN=2048, fs=fs) - fase = np.unwrap(np.angle(H))*180/np.pi - - fig, ax1, ax2 = freq_response_plot(w, H, fase, show=False) - save_plot(fig, "respuesta_en_frecuencia_pasa-bajos_fir") - -# `a` son los coeficientes de la respuesta al impulso (coinciden con los -# coeficientes de respuesta en frecuencia) -def filtro_fir_polos_y_ceros(a): - zeros, poles, gain = tf2zpk(a, [1]) - - # Crear figura - fig, axis = plt.subplots(figsize=(8, 4)) - - axis.scatter(np.real(zeros), np.imag(zeros), - s=25, facecolors='none', edgecolors='tab:blue', zorder=10, - label='Ceros', linewidth=1.25) - - axis.scatter(np.real(poles), np.imag(poles), - s=25, marker='x', color='tab:red', - label='Polos') - - axis.set_xlabel("Real", color="black") - axis.set_ylabel("Imaginario", color="black") - - # Unidad círculo para referencia - # theta = np.linspace(0, 2*np.pi, 100) - # plt.plot(np.cos(theta), np.sin(theta)) # círculo unitario - - # axis.yaxis.set_major_locator(MaxNLocator(nbins=5)) - - axis.grid(True, which='major', color='black', linestyle=':', linewidth=1.00) - axis.grid(True, which='minor', color='black', linestyle=':', linewidth=0.50) - axis.xaxis.set_minor_locator(AutoMinorLocator(2)) - - plt.grid(True) - plt.axis('equal') - axis.legend() - save_plot(fig, "polos_y_ceros_pasa-bajos_fir") - - -def filtro_fir_deducido(): - fs = 44100 - - # respuesta ideal pasabajos: sinc centrada en M/2 - n = np.arange(M + 1) - wc = 2*np.pi*cutoff / fs - - # h_ideal = sinc(wc*n)/(pi n); wc = 2pi*fc/fs - # se normaliza la ganancia a 1 multiplicando por 2.0*(fc/fs) - h_ideal = np.sinc(2.0 * (cutoff/fs) * (n - M/2)) - - # ventana de Hamming - v_hamming = 0.54 - 0.46 * np.cos(2*np.pi*n/M) - - v_rectangular = [] - for i in n: - v_rectangular.append(1 if i < 350 else 0) - - # respuesta del filtro FIR (version acotada de la sinc) - h = h_ideal * v_hamming - - # se normaliza para tener ganancia unitaria para frecuancias <= fc - h = h / np.sum(h) - - # fig, ax = dtime_plot(M, h, "respuesta_al_impulso_filtro_fir", - # f'Respuesta al impulso filtro FIR grado {M}') - - # respuesta en frecuencia del filtro - w, H = freqz(h, worN=2048, fs=fs) - fase = np.unwrap(np.angle(H)) * 180 / np.pi - - # fig, ax1, ax2 = freq_response_plot(w, H, fase, show=False, fc=5e3) - # save_plot(fig, "respuesta_en_frecuencia_pasa-bajos_fir") - - # polos y ceros - # filtro_fir_polos_y_ceros(h) - - file_path = 'canciones/000002.mp3' - - file_data, file_fs = librosa.load(file_path, sr=None, mono=True) - file_fs = int(fs) - - # spectogram_plot(file_fs, file_data, - # f"espectograma_fs_original_44100Hz", N=1024, ylim=[0, 20000]) - - file_filter_output = np.convolve(file_data, h, mode='same') - - # spectogram_plot(file_fs, file_filter_output, - # f"espectograma_fs_{cutoff}Hz", t=0, N=1024, ylim=[0, 20000]) - - freq_plot(44100, v_hamming, "v_hamming_freq", f_max=8000) - freq_plot(44100, v_rectangular, "v_rectangular_freq", f_max=8000) - - freq_compute_fft(44100, v_hamming) - - for i in [512, 1024, 2048]: - # for window in ['boxcar', 'bartlett', 'hamming']: - - # spectogram_plot(file1_fs, file_filter_output, - # f"espectograma_submuestreado_{window}_{i:04d}", N=i, - # win=window, ylim=[0, 3000], t=5, dt=1) - - N = 1024 - beta = 8.6 - kaiser_window = get_window(("kaiser", beta), N) - spectogram_plot(file1_fs, file_filter_output, - f"espectograma_submuestreado_kaiser_window_{i:04d}", N=i, - win=kaiser_window, ylim=[0, 3000], t=5, dt=1) - -############################# Llamados a funciones ############################ - -def primer_y_segunda_parte(): - # time_domain - time_domain_cancion1() - time_domain_cancion2() - time_domain_music_instruments() - - # freq_domain - freq_domain_cancion1() - freq_domain_cancion2() - freq_domain_spectograms() - - freq_plot(48000, filter1_h, "filter1_h_fft", f_max=2000) - freq_plot(48000, filter2_h, "filter2_h_fft", f_max=8000) - - a4_flauta_fseries() - a4_clarinete_fseries() - a4_violin_fseries() - - # obs: para realizar el filtrado se toma toda la señal no solo un periodo - a4_flauta_cutoff() - a4_clarinete_cutoff() - a4_violin_cutoff() - - -# primer_y_segunda_parte() -filtro_fir_deducido() diff --git a/tp/plot/a4_clarinete.png b/tp/plot/a4_clarinete.png Binary files differ. diff --git a/tp/plot/a4_clarinete_comparison.png b/tp/plot/a4_clarinete_comparison.png Binary files differ. diff --git a/tp/plot/a4_clarinete_cutoff_3000Hz.png b/tp/plot/a4_clarinete_cutoff_3000Hz.png Binary files differ. diff --git a/tp/plot/a4_clarinete_cutoff_3000Hz_fft.png b/tp/plot/a4_clarinete_cutoff_3000Hz_fft.png Binary files differ. diff --git a/tp/plot/a4_clarinete_cutoff_time_comparison.png b/tp/plot/a4_clarinete_cutoff_time_comparison.png Binary files differ. diff --git a/tp/plot/a4_clarinete_fseries_comparison.png b/tp/plot/a4_clarinete_fseries_comparison.png Binary files differ. diff --git a/tp/plot/a4_flauta.png b/tp/plot/a4_flauta.png Binary files differ. diff --git a/tp/plot/a4_flauta_comparison.png b/tp/plot/a4_flauta_comparison.png Binary files differ. diff --git a/tp/plot/a4_flauta_cutoff_1000Hz.png b/tp/plot/a4_flauta_cutoff_1000Hz.png Binary files differ. diff --git a/tp/plot/a4_flauta_cutoff_1000Hz_fft.png b/tp/plot/a4_flauta_cutoff_1000Hz_fft.png Binary files differ. diff --git a/tp/plot/a4_flauta_cutoff_time_comparison.png b/tp/plot/a4_flauta_cutoff_time_comparison.png Binary files differ. diff --git a/tp/plot/a4_flauta_fseries_comparison.png b/tp/plot/a4_flauta_fseries_comparison.png Binary files differ. diff --git a/tp/plot/a4_violin.png b/tp/plot/a4_violin.png Binary files differ. diff --git a/tp/plot/a4_violin_comparison.png b/tp/plot/a4_violin_comparison.png Binary files differ. diff --git a/tp/plot/a4_violin_cutoff_4000Hz.png b/tp/plot/a4_violin_cutoff_4000Hz.png Binary files differ. diff --git a/tp/plot/a4_violin_cutoff_4000Hz_fft.png b/tp/plot/a4_violin_cutoff_4000Hz_fft.png Binary files differ. diff --git a/tp/plot/a4_violin_cutoff_time_comparison.png b/tp/plot/a4_violin_cutoff_time_comparison.png Binary files differ. diff --git a/tp/plot/a4_violin_fseries_comparison.png b/tp/plot/a4_violin_fseries_comparison.png Binary files differ. diff --git a/tp/plot/cancion1.png b/tp/plot/cancion1.png Binary files differ. diff --git a/tp/plot/cancion1_0_248s_a_0_256s.png b/tp/plot/cancion1_0_248s_a_0_256s.png Binary files differ. diff --git a/tp/plot/cancion1_0_520s_a_0_528s.png b/tp/plot/cancion1_0_520s_a_0_528s.png Binary files differ. diff --git a/tp/plot/cancion1_espectograma_bartlett_0512.png b/tp/plot/cancion1_espectograma_bartlett_0512.png Binary files differ. diff --git a/tp/plot/cancion1_espectograma_bartlett_1024.png b/tp/plot/cancion1_espectograma_bartlett_1024.png Binary files differ. diff --git a/tp/plot/cancion1_espectograma_bartlett_2048.png b/tp/plot/cancion1_espectograma_bartlett_2048.png Binary files differ. diff --git a/tp/plot/cancion1_espectograma_boxcar_0512.png b/tp/plot/cancion1_espectograma_boxcar_0512.png Binary files differ. diff --git a/tp/plot/cancion1_espectograma_boxcar_1024.png b/tp/plot/cancion1_espectograma_boxcar_1024.png Binary files differ. diff --git a/tp/plot/cancion1_espectograma_boxcar_2048.png b/tp/plot/cancion1_espectograma_boxcar_2048.png Binary files differ. diff --git a/tp/plot/cancion1_espectograma_hamming_0512.png b/tp/plot/cancion1_espectograma_hamming_0512.png Binary files differ. diff --git a/tp/plot/cancion1_espectograma_hamming_1024.png b/tp/plot/cancion1_espectograma_hamming_1024.png Binary files differ. diff --git a/tp/plot/cancion1_espectograma_hamming_2048.png b/tp/plot/cancion1_espectograma_hamming_2048.png Binary files differ. diff --git a/tp/plot/cancion1_espectograma_hann_0512.png b/tp/plot/cancion1_espectograma_hann_0512.png Binary files differ. diff --git a/tp/plot/cancion1_espectograma_hann_1024.png b/tp/plot/cancion1_espectograma_hann_1024.png Binary files differ. diff --git a/tp/plot/cancion1_espectograma_hann_2048.png b/tp/plot/cancion1_espectograma_hann_2048.png Binary files differ. diff --git a/tp/plot/cancion1_fft.png b/tp/plot/cancion1_fft.png Binary files differ. diff --git a/tp/plot/cancion1_filter1_output_compare.png b/tp/plot/cancion1_filter1_output_compare.png Binary files differ. diff --git a/tp/plot/cancion1_filter1_output_compare_0_248_a_0_256.png b/tp/plot/cancion1_filter1_output_compare_0_248_a_0_256.png Binary files differ. diff --git a/tp/plot/cancion1_filter1_output_fft.png b/tp/plot/cancion1_filter1_output_fft.png Binary files differ. diff --git a/tp/plot/cancion1_filter2_output_compare.png b/tp/plot/cancion1_filter2_output_compare.png Binary files differ. diff --git a/tp/plot/cancion1_filter2_output_compare_0_248_a_0_256.png b/tp/plot/cancion1_filter2_output_compare_0_248_a_0_256.png Binary files differ. diff --git a/tp/plot/cancion1_filter2_output_fft.png b/tp/plot/cancion1_filter2_output_fft.png Binary files differ. diff --git a/tp/plot/cancion2.png b/tp/plot/cancion2.png Binary files differ. diff --git a/tp/plot/cancion2_14_72s_a_14_73s.png b/tp/plot/cancion2_14_72s_a_14_73s.png Binary files differ. diff --git a/tp/plot/cancion2_26_57s_a_26_58s.png b/tp/plot/cancion2_26_57s_a_26_58s.png Binary files differ. diff --git a/tp/plot/cancion2_6s_filter1_output_compare.png b/tp/plot/cancion2_6s_filter1_output_compare.png Binary files differ. diff --git a/tp/plot/cancion2_6s_filter1_output_compare_26_57_a_26_58.png b/tp/plot/cancion2_6s_filter1_output_compare_26_57_a_26_58.png Binary files differ. diff --git a/tp/plot/cancion2_6s_filter2_output_compare.png b/tp/plot/cancion2_6s_filter2_output_compare.png Binary files differ. diff --git a/tp/plot/cancion2_6s_filter2_output_compare_26_57_a_26_58.png b/tp/plot/cancion2_6s_filter2_output_compare_26_57_a_26_58.png Binary files differ. diff --git a/tp/plot/cancion2_espectograma_bartlett_0512.png b/tp/plot/cancion2_espectograma_bartlett_0512.png Binary files differ. diff --git a/tp/plot/cancion2_espectograma_bartlett_1024.png b/tp/plot/cancion2_espectograma_bartlett_1024.png Binary files differ. diff --git a/tp/plot/cancion2_espectograma_bartlett_2048.png b/tp/plot/cancion2_espectograma_bartlett_2048.png Binary files differ. diff --git a/tp/plot/cancion2_espectograma_boxcar_0512.png b/tp/plot/cancion2_espectograma_boxcar_0512.png Binary files differ. diff --git a/tp/plot/cancion2_espectograma_boxcar_1024.png b/tp/plot/cancion2_espectograma_boxcar_1024.png Binary files differ. diff --git a/tp/plot/cancion2_espectograma_boxcar_2048.png b/tp/plot/cancion2_espectograma_boxcar_2048.png Binary files differ. diff --git a/tp/plot/cancion2_espectograma_hamming_0512.png b/tp/plot/cancion2_espectograma_hamming_0512.png Binary files differ. diff --git a/tp/plot/cancion2_espectograma_hamming_1024.png b/tp/plot/cancion2_espectograma_hamming_1024.png Binary files differ. diff --git a/tp/plot/cancion2_espectograma_hamming_2048.png b/tp/plot/cancion2_espectograma_hamming_2048.png Binary files differ. diff --git a/tp/plot/cancion2_espectograma_hann_0512.png b/tp/plot/cancion2_espectograma_hann_0512.png Binary files differ. diff --git a/tp/plot/cancion2_espectograma_hann_1024.png b/tp/plot/cancion2_espectograma_hann_1024.png Binary files differ. diff --git a/tp/plot/cancion2_espectograma_hann_2048.png b/tp/plot/cancion2_espectograma_hann_2048.png Binary files differ. diff --git a/tp/plot/cancion2_fft.png b/tp/plot/cancion2_fft.png Binary files differ. diff --git a/tp/plot/cancion2_filter1_output_fft.png b/tp/plot/cancion2_filter1_output_fft.png Binary files differ. diff --git a/tp/plot/cancion2_filter2_output_fft.png b/tp/plot/cancion2_filter2_output_fft.png Binary files differ. diff --git a/tp/plot/filter1_h_fft.png b/tp/plot/filter1_h_fft.png Binary files differ. diff --git a/tp/plot/filter2_h_fft.png b/tp/plot/filter2_h_fft.png Binary files differ. diff --git a/tp/plot/polos_y_ceros_pasa-bajos_fir.png b/tp/plot/polos_y_ceros_pasa-bajos_fir.png Binary files differ. diff --git a/tp/plot/polos_y_ceros_pasa-bajos_fir_cero_marcado.png b/tp/plot/polos_y_ceros_pasa-bajos_fir_cero_marcado.png Binary files differ. diff --git a/tp/plot/respuesta_al_impulso_filtro_fir.png b/tp/plot/respuesta_al_impulso_filtro_fir.png Binary files differ. diff --git a/tp/plot/respuesta_en_frecuencia_pasa-bajos_fir.png b/tp/plot/respuesta_en_frecuencia_pasa-bajos_fir.png Binary files differ. diff --git a/tp/scripts/data.py b/tp/scripts/data.py @@ -0,0 +1,45 @@ +from scipy.io import wavfile +import numpy as np +import librosa +import os + +data_dir = '../data/' +plot_dir = '../plot/' +out_dir = '../out/' + +## Datos +file1_path = data_dir + 'cancion1.wav' +file2_path = data_dir + 'cancion2.wav' +filter1_h_file_path = data_dir + 'respuesta_impulso_1.txt' +filter2_h_file_path = data_dir + 'respuesta_impulso_2.txt' + +a4_flauta_file_path = data_dir + 'a4_flauta.wav' +a4_clarinete_file_path = data_dir + 'a4_clarinete.wav' +a4_violin_file_path = data_dir + 'a4_violin.wav' + +file1_fs, file1_data = wavfile.read(file1_path) +file2_fs, file2_data = wavfile.read(file2_path) +filter1_h = np.loadtxt(filter1_h_file_path) +filter2_h = np.loadtxt(filter2_h_file_path) + +a4_flauta_fs, a4_flauta_data = wavfile.read(a4_flauta_file_path) +a4_clarinete_fs, a4_clarinete_data = wavfile.read(a4_clarinete_file_path) +a4_violin_fs, a4_violin_data = wavfile.read(a4_violin_file_path) + +file1_filter1_output = np.convolve(file1_data, filter1_h, mode='same') +file1_filter2_output = np.convolve(file1_data, filter2_h, mode='same') +file2_filter1_output = np.convolve(file2_data, filter1_h, mode='same') +file2_filter2_output = np.convolve(file2_data, filter2_h, mode='same') + +## Datos tercera parte + +# todas las canciones tienen formato mp3 y 44100Hz de frecuencia de muestreo +canciones_dataset_dir = data_dir + 'canciones/' +canciones_dataset = [] + +for filename in os.listdir(canciones_dataset_dir): + filename = canciones_dataset_dir + filename + print(filename) + file_data, file_fs = librosa.load(data_dir + filename, sr=None, mono=True) + canciones_dataset.append(file_data) + file_fs = int(file_fs) diff --git a/tp/scripts/main.py b/tp/scripts/main.py @@ -0,0 +1,41 @@ +import numpy as np + +# archivos locales +from utils import * +from primera_parte import * +from segunda_parte import * +from tercera_parte import * +from data import * + +############################# Llamados a funciones ############################ + +def primera_parte(): + # time_domain + time_domain_cancion1() + time_domain_cancion2() + time_domain_music_instruments() + +def segunda_parte(): + freq_domain + freq_domain_cancion1() + freq_domain_cancion2() + freq_domain_spectograms() + + freq_plot(48000, filter1_h, "filter1_h_fft", f_max=2000) + freq_plot(48000, filter2_h, "filter2_h_fft", f_max=8000) + + a4_flauta_fseries() + a4_clarinete_fseries() + a4_violin_fseries() + + # obs: para realizar el filtrado se toma toda la señal no solo un periodo + a4_flauta_cutoff() + a4_clarinete_cutoff() + a4_violin_cutoff() + +def tercera_parte(): + filtro_fir_deducido() + +# primera_parte() +# segunda_parte() +tercera_parte() diff --git a/tp/scripts/primera_parte.py b/tp/scripts/primera_parte.py @@ -0,0 +1,79 @@ +from utils import * +from data import * + +############################## Primera Parte ################################### + +def time_domain_cancion1(): + ### grafico completo + time_plot(file1_fs, file1_data, file1_path) + + ### porciones cuasi-periodicas 'cancion1' + time_plot(file1_fs, file1_data, "cancion1_0_248s_a_0_256s", + t=0.248, dt=0.008, a=0.24978, da=0.003) + + time_plot(file1_fs, file1_data, "cancion1_0_520s_a_0_528s", + t=0.520, dt=0.008, a=0.5208, da=0.003) + + ### salida de filtro 'cancion1' + save_to_wav(file1_fs, file1_filter1_output, "file1_filter1_output.wav") + save_to_wav(file1_fs, file1_filter2_output, "file1_filter2_output.wav") + + ### grafico comparando la muestra 1 original y filtrada 1 + data_arr = [normalize(file1_data), normalize(file1_filter1_output)] + leg_arr = ['Señal de audio', 'Señal de audio filtrada'] + + time_plot_multiple(file1_fs, data_arr, leg_arr, + "cancion1_filter1_output_compare") + + time_plot_multiple(file1_fs, data_arr, leg_arr, + "cancion1_filter1_output_compare_0_248_a_0_256", + t=0.248, dt=0.008) + + ### grafico comparando la muestra 1 original y filtrada 2 + data_arr = [normalize(file1_data), normalize(file1_filter2_output)] + leg_arr = ['Señal de audio', 'Señal de audio filtrada'] + + time_plot_multiple(file1_fs, data_arr, leg_arr, + "cancion1_filter2_output_compare") + time_plot_multiple(file1_fs, data_arr, leg_arr, + "cancion1_filter2_output_compare_0_248_a_0_256", + t=0.248, dt=0.008) + +def time_domain_cancion2(): + ### grafico completo + time_plot(file2_fs, file2_data, "cancion2", t=6) + + ### porciones cuasi-periodicas 'cancion2' + time_plot(file2_fs, file2_data, "cancion2_14_72s_a_14_73s", t=14.720, dt=0.01) + time_plot(file2_fs, file2_data, "cancion2_26_57s_a_26_58s", t=26.570, dt=0.01) + + save_to_wav(file2_fs, file2_filter1_output, "file2_filter1_output.wav") + save_to_wav(file2_fs, file2_filter2_output, "file2_filter2_output.wav") + + ### grafico comparando la muestra 2 original y filtrada 2 + data_arr = [normalize(file2_data), normalize(file2_filter1_output)] + leg_arr = ['Señal original', 'Señal filtrada'] + + time_plot_multiple(file2_fs, data_arr, leg_arr, + "cancion2_6s_filter1_output_compare", t=6) + time_plot_multiple(file2_fs, data_arr, leg_arr, + "cancion2_6s_filter1_output_compare_26_57_a_26_58", + t=26.57, dt=0.01) + + ### grafico comparando la muestra 2 original y filtrada 2 + data_arr = [normalize(file2_data), normalize(file2_filter2_output)] + leg_arr = ['Señal original', 'Señal filtrada'] + + time_plot_multiple(file2_fs, data_arr, leg_arr, + "cancion2_6s_filter2_output_compare", t=6) + time_plot_multiple(file2_fs, data_arr, leg_arr, + "cancion2_6s_filter2_output_compare_26_57_a_26_58", + t=26.57, dt=0.01) + +def time_domain_music_instruments(): + ### grafico de los instrumentos musicales + time_plot(a4_flauta_fs, a4_flauta_data, "a4_flauta", t=0.25, dt=0.010) + time_plot(a4_clarinete_fs, a4_clarinete_data, "a4_clarinete", t=0.25, dt=0.010) + time_plot(a4_violin_fs, a4_violin_data, "a4_violin", t=0.25, dt=0.010) + + diff --git a/tp/scripts/segunda_parte.py b/tp/scripts/segunda_parte.py @@ -0,0 +1,179 @@ +from utils import * + +############################## Segunda parte ################################## + +def freq_domain_cancion1(): + freq_plot(file1_fs, file1_data, "cancion1_fft", f_max=8000) + freq_plot(file1_fs, file1_filter1_output, "cancion1_filter1_output_fft", + f_max=8000) + freq_plot(file1_fs, file1_filter2_output, "cancion1_filter2_output_fft", + f_max=8000) + +def freq_domain_cancion2(): + freq_plot(file2_fs, file2_data, "cancion2_fft", + f_max=8000) + freq_plot(file2_fs, file2_filter1_output, "cancion2_filter1_output_fft", + f_max=8000) + freq_plot(file2_fs, file2_filter2_output, "cancion2_filter2_output_fft", + f_max=8000) + +def freq_domain_spectograms(): + # Formas de funcion ventana (en tiempo, en frecuencia tienen otra forma) + # 'boxcar': rectangular + # 'bartlett': triangular + # 'hann': similar a medio ciclo de seno + # https://en.wikipedia.org/wiki/Window_function + + for i in [512, 1024, 2048]: + for window in ['boxcar', 'bartlett', 'hann', 'hamming']: + spectogram_plot(file1_fs, file1_data, + f"cancion1_espectograma_{window}_{i:04d}", N=i, + win=window, ylim=[0, 17500]) + + spectogram_plot(file2_fs, file2_data, + f"cancion2_espectograma_{window}_{i:04d}", t=6, N=i, + win=window, ylim=[0, 8000]) + +def a4_flauta_cutoff(): + a4_flauta_cutoff_fft, a4_flauta_cutoff_freqs = freq_compute_fft( + a4_flauta_fs, a4_flauta_data) + + # armonicos mayores a 'cutoff_freq' Hz son descartadas, analogo a aplicar un + # filtro pasabajos ideal + cutoff_freq = 1000 + a4_flauta_cutoff_fft[np.abs(a4_flauta_cutoff_freqs) > cutoff_freq] = 0.0 + + # espectro de la señal filtrada + N = len(a4_flauta_cutoff_fft) + x = a4_flauta_cutoff_freqs[:N // 2] + y = np.abs(a4_flauta_cutoff_fft[:N // 2]) + fig, axis = freq_graph_data(x, y, f_max=2000, show=False) + save_plot(fig, f"a4_flauta_cutoff_{cutoff_freq}Hz_fft") + + # señal temporal reconstruida + a4_flauta_cutoff = ifft(a4_flauta_cutoff_fft).real + + time_plot(a4_flauta_fs, a4_flauta_cutoff, + f"a4_flauta_cutoff_{cutoff_freq}Hz", 0.25, 0.010) + + save_to_wav(a4_flauta_fs, a4_flauta_cutoff, + f"a4_flauta_cutoff_{cutoff_freq}Hz.wav") + + data_arr = [normalize(a4_flauta_data), normalize(a4_flauta_cutoff)] + leg_arr = [ + 'Señal de nota A4 de flauta', + 'Señal de nota A4 de flauta filtrada' + ] + time_plot_multiple(a4_flauta_fs, data_arr, leg_arr, + "a4_flauta_cutoff_time_comparison", t=0.25, dt=0.008) + + data_arr = [a4_flauta_cutoff, a4_flauta_data] + leg_arr = ["Nota musical A4 con flauta filtrada", "Nota musical A4 con flauta"] + freq_plot_multiple(a4_flauta_fs, data_arr, leg_arr, f_max=4000, t=0.253, + dt=8*0.002272727, save_name="a4_flauta_comparison", show=False) + +def a4_clarinete_cutoff(): + a4_clarinete_cutoff_fft, a4_clarinete_cutoff_freqs = freq_compute_fft( + a4_clarinete_fs, a4_clarinete_data) + + cutoff_freq = 3000 + a4_clarinete_cutoff_fft[np.abs(a4_clarinete_cutoff_freqs) > cutoff_freq] = 0.0 + + # espectro de la señal filtrada + N = len(a4_clarinete_cutoff_fft) + x = a4_clarinete_cutoff_freqs[:N // 2] + y = np.abs(a4_clarinete_cutoff_fft[:N // 2]) + fig, axis = freq_graph_data(x, y, f_max=3000, show=False) + save_plot(fig, f"a4_clarinete_cutoff_{cutoff_freq}Hz_fft") + + # señal temporal reconstruida + a4_clarinete_cutoff = ifft(a4_clarinete_cutoff_fft).real + time_plot(a4_clarinete_fs, a4_clarinete_cutoff, + f"a4_clarinete_cutoff_{cutoff_freq}Hz", 0.25, 0.010) + save_to_wav(a4_clarinete_fs, a4_clarinete_cutoff, + f"a4_clarinete_cutoff_{cutoff_freq}Hz.wav") + + data_arr = [normalize(a4_clarinete_data), normalize(a4_clarinete_cutoff)] + leg_arr = [ + 'Señal de nota A4 de clarinete', + 'Señal de nota A4 de clarinete filtrada' + ] + time_plot_multiple(a4_clarinete_fs, data_arr, leg_arr, + "a4_clarinete_cutoff_time_comparison", t=0.25, dt=0.008) + + data_arr = [a4_clarinete_cutoff, a4_clarinete_data] + leg_arr = [ + "Nota musical A4 con clarinete filtrada", + "Nota musical A4 con clarinete" + ] + + freq_plot_multiple(a4_clarinete_fs, data_arr, leg_arr, f_max=8000, t=0.253, + dt=8*0.002272727, save_name="a4_clarinete_comparison", show=False) + +def a4_violin_cutoff(): + a4_violin_cutoff_fft, a4_violin_cutoff_freqs = freq_compute_fft( + a4_violin_fs, a4_violin_data) + + cutoff_freq = 4000 + a4_violin_cutoff_fft[np.abs(a4_violin_cutoff_freqs) > cutoff_freq] = 0.0 + + # espectro + N = len(a4_violin_cutoff_fft) + x = a4_violin_cutoff_freqs[:N // 2] + y = np.abs(a4_violin_cutoff_fft[:N // 2]) + fig, axis = freq_graph_data(x, y, f_max=4000, show=False) + save_plot(fig, f"a4_violin_cutoff_{cutoff_freq}Hz_fft") + + # señal temporal reconstruida + a4_violin_cutoff = ifft(a4_violin_cutoff_fft).real + time_plot(a4_violin_fs, a4_violin_cutoff, + f"a4_violin_cutoff_{cutoff_freq}Hz", 0.25, 0.010) + save_to_wav(a4_violin_fs, a4_violin_cutoff, f"a4_violin_cutoff_{cutoff_freq}Hz.wav") + + data_arr = [normalize(a4_violin_data), normalize(a4_violin_cutoff)] + leg_arr = ['Señal de nota A4 de violin', 'Señal de nota A4 de violin filtrada'] + time_plot_multiple(a4_violin_fs, data_arr, leg_arr, + "a4_violin_cutoff_time_comparison", t=0.25, dt=0.008) + + data_arr = [a4_violin_cutoff, a4_violin_data] + leg_arr = ["Nota musical A4 con violin filtrada", "Nota musical A4 con violin"] + freq_plot_multiple(a4_violin_fs, data_arr, leg_arr, f_max=8000, t=0.253, + dt=8*0.002272727, save_name="a4_violin_comparison", show=False) + +def a4_flauta_fseries(): + fft_freqs_arr = [] + + for i in [8, 4, 1]: + fft, freqs = freq_compute_fft(a4_flauta_fs, a4_flauta_data, t=0.253, dt=i*0.002272727) + fft_freqs_arr.append([fft, freqs, + f"Serie de Fourier {i} periodo{'s' if i != 1 else ''}"]) + + fig, axis = freq_graph_multiple_data(fft_freqs_arr, show=False, f_max=4500) + save_plot(fig, "a4_flauta_fseries_comparison") + +def a4_clarinete_fseries(): + fft_freqs_arr = [] + + for i in [8, 4, 1]: + fft, freqs = freq_compute_fft( + a4_clarinete_fs, a4_clarinete_data, t=0.253, dt=i*0.002272727) + + fft_freqs_arr.append([fft, freqs, + f"Serie de Fourier {i} periodo{'s' if i != 1 else ''}"]) + + fig, axis = freq_graph_multiple_data(fft_freqs_arr, show=False, f_max=8000) + save_plot(fig, "a4_clarinete_fseries_comparison") + +def a4_violin_fseries(): + fft_freqs_arr = [] + + for i in [8, 4, 1]: + fft, freqs = freq_compute_fft(a4_violin_fs, a4_violin_data, + t=0.253, dt=i*0.002272727) + fft_freqs_arr.append([fft, freqs, + f"Serie de Fourier {i} periodo{'s' if i != 1 else ''}"]) + + fig, axis = freq_graph_multiple_data(fft_freqs_arr, show=False, f_max=8000) + save_plot(fig, "a4_violin_fseries_comparison") + + diff --git a/tp/scripts/tercera_parte.py b/tp/scripts/tercera_parte.py @@ -0,0 +1,124 @@ +from utils import * +from scipy.signal import firwin, freqz, tf2zpk, get_window + +############################# Tercera parte ################################### + +cutoff = 2650 # frecuencia de corte +M = 700 # orden FIR (número de coeficientes) + +def filtro_fir(): + # Diseño FIR pasabajos con ventana + b = firwin(M, cutoff, fs=fs, window='hamming') + + w, H = freqz(b, worN=2048, fs=fs) + fase = np.unwrap(np.angle(H))*180/np.pi + + fig, ax1, ax2 = freq_response_plot(w, H, fase, show=False) + save_plot(fig, "respuesta_en_frecuencia_pasa-bajos_fir") + +# `a` son los coeficientes de la respuesta al impulso (coinciden con los +# coeficientes de respuesta en frecuencia) +def filtro_fir_polos_y_ceros(a): + zeros, poles, gain = tf2zpk(a, [1]) + + # Crear figura + fig, axis = plt.subplots(figsize=(8, 4)) + + axis.scatter(np.real(zeros), np.imag(zeros), + s=25, facecolors='none', edgecolors='tab:blue', zorder=10, + label='Ceros', linewidth=1.25) + + axis.scatter(np.real(poles), np.imag(poles), + s=25, marker='x', color='tab:red', + label='Polos') + + axis.set_xlabel("Real", color="black") + axis.set_ylabel("Imaginario", color="black") + + # Unidad círculo para referencia + # theta = np.linspace(0, 2*np.pi, 100) + # plt.plot(np.cos(theta), np.sin(theta)) # círculo unitario + + # axis.yaxis.set_major_locator(MaxNLocator(nbins=5)) + + axis.grid(True, which='major', color='black', linestyle=':', linewidth=1.00) + axis.grid(True, which='minor', color='black', linestyle=':', linewidth=0.50) + axis.xaxis.set_minor_locator(AutoMinorLocator(2)) + + plt.grid(True) + plt.axis('equal') + axis.legend() + save_plot(fig, "polos_y_ceros_pasa-bajos_fir") + + +fs = 44100 + +# a diferencia de la funcion `filtro_fir` se genera el filtro mediante +# operaciones elementales, como la multiplicacion por ventana, en lugar de usar +# una funcion de libreria externa como `firwin` +def filtro_fir_deducido(): + # respuesta ideal pasabajos: sinc centrada en M/2 + n = np.arange(M + 1) + wc = 2*np.pi*cutoff / fs + + # h_ideal = sinc(wc*n)/(pi n); wc = 2pi*fc/fs + # se normaliza la ganancia a 1 multiplicando por 2.0*(fc/fs) + h_ideal = np.sinc(2.0 * (cutoff/fs) * (n - M/2)) + + # ventana de Hamming + v_hamming = 0.54 - 0.46 * np.cos(2*np.pi*n/M) + + v_rectangular = [1 if i < 350 else 0 for i in range(0, M)] + + # respuesta del filtro FIR (version acotada de la sinc) + h = h_ideal * v_hamming + + # se normaliza para tener ganancia unitaria para frecuancias <= fc + h = h / np.sum(h) + return h + +def filtro_fir_analisis(h, fs): + fig, ax = dtime_plot(M, h, "respuesta_al_impulso_filtro_fir", + f'Respuesta al impulso filtro FIR grado {M}') + + # respuesta en frecuencia del filtro + w, H = freqz(h, worN=2048, fs=fs) + fase = np.unwrap(np.angle(H)) * 180 / np.pi + + fig, ax1, ax2 = freq_response_plot(w, H, fase, show=False, fc=5e3) + save_plot(fig, "respuesta_en_frecuencia_pasa-bajos_fir") + + # polos y ceros + filtro_fir_polos_y_ceros(h) + + + +""" +def ejemplo_cancion_filtrado_con_filtro_fir(): + # spectogram_plot(file_fs, file_data, + # f"espectograma_fs_original_44100Hz", N=1024, ylim=[0, 20000]) + + file_filter_output = np.convolve(file_data, h, mode='same') + + # spectogram_plot(file_fs, file_filter_output, + # f"espectograma_fs_{cutoff}Hz", t=0, N=1024, ylim=[0, 20000]) + + freq_plot(44100, v_hamming, "v_hamming_freq", f_max=8000) + freq_plot(44100, v_rectangular, "v_rectangular_freq", f_max=8000) + + freq_compute_fft(44100, v_hamming) + + for i in [512, 1024, 2048]: + # for window in ['boxcar', 'bartlett', 'hamming']: + + # spectogram_plot(file1_fs, file_filter_output, + # f"espectograma_submuestreado_{window}_{i:04d}", N=i, + # win=window, ylim=[0, 3000], t=5, dt=1) + + N = 1024 + beta = 8.6 + kaiser_window = get_window(("kaiser", beta), N) + spectogram_plot(file1_fs, file_filter_output, + f"espectograma_submuestreado_kaiser_window_{i:04d}", N=i, + win=kaiser_window, ylim=[0, 3000], t=5, dt=1) +""" diff --git a/tp/scripts/utils.py b/tp/scripts/utils.py @@ -0,0 +1,433 @@ +import matplotlib.pyplot as plt +from matplotlib.ticker import AutoMinorLocator +from matplotlib.ticker import MultipleLocator +from matplotlib.ticker import MaxNLocator +from matplotlib.ticker import FuncFormatter +from cycler import cycler + +import matplotlib +import numpy as np +import os + +from scipy.io import wavfile +from scipy.fft import fft, ifft, fftfreq +from scipy.signal import spectrogram + +from data import * + +matplotlib.rcParams['font.family'] = 'Inter' +matplotlib.rcParams['font.size'] = 12 +matplotlib.rcParams['axes.prop_cycle'] = cycler( + color=['#1f77b4', '#ff0000', '#ff5f1f', 'green']) +matplotlib.use("TkAgg") + +def ticks_label_format(x, pos): + # 3 decimales, se eliminan los ceros y puntos + return f"{x:.3f}".rstrip("0").rstrip(".") + +def time_graph_multiple_data(x, y_arr, y_lab, t=0, dt=0, a=0, da=0, show=True): + figure, axis = plt.subplots(figsize=(8, 4)) + + for i, y in enumerate(y_arr): + axis.plot(x, y, label=y_lab[i], alpha=0.75) + + axis.set(xlabel='Tiempo [s]', ylabel='Amplitud normalizada') + + axis.minorticks_on() + axis.grid(True, which='major', color='black', linestyle=':', linewidth=1.00) + axis.grid(True, which='minor', color='black', linestyle=':', linewidth=0.50) + + # configuracion de ticks del eje x + axis.xaxis.set_major_locator(MaxNLocator(nbins=5)) + axis.xaxis.set_minor_locator(AutoMinorLocator(5)) + + axis.yaxis.set_major_locator(MaxNLocator(nbins=5)) + axis.yaxis.set_minor_locator(AutoMinorLocator(4)) + + plt.tight_layout() + + # max 3 decimales + axis.xaxis.set_major_formatter(FuncFormatter(ticks_label_format)) + + axis.set_xlim([t, t+dt if dt > 0 else x[-1]]) + axis.set_ylim(-1.1, 1.1) + + # resaltado de parte de la señal (solo si a != 0) + axis.axvspan(a, a+da, color='skyblue', + alpha=0 if a == 0 else 0.50, + label=f"Un periodo T={da}s" if da != 0 else "") + axis.legend(loc='upper left') + + if show: + plt.show() + + return figure, axis + +# todos deben la misma cantidad de elementos que el primero +def time_plot_multiple(fs, data_arr, leg_arr, save_name="", t=0, dt=0, a=0, da=0, show=False): + x = np.arange(len(data_arr[0])) / fs + fig, ax = time_graph_multiple_data(x, data_arr, leg_arr, t, dt, a=a, da=da, show=show) + + if show == False: + save_plot(fig, save_name) + + return fig, ax + +def time_graph_data(x, y, t=0, dt=0, a=0, da=0, show=True): + figure, axis = plt.subplots(figsize=(8, 4)) + + axis.plot(x, y, label='Señal de audio') + axis.set(xlabel='Tiempo [s]', ylabel='Amplitud normalizada') + + axis.minorticks_on() + axis.grid(True, which='major', color='black', linestyle=':', linewidth=1.00) + axis.grid(True, which='minor', color='black', linestyle=':', linewidth=0.50) + + # configuracion de ticks del eje x + axis.xaxis.set_major_locator(MaxNLocator(nbins=5)) + axis.xaxis.set_minor_locator(AutoMinorLocator(5)) + + axis.yaxis.set_major_locator(MaxNLocator(nbins=5)) + axis.yaxis.set_minor_locator(AutoMinorLocator(4)) + + plt.tight_layout() + + # max 3 decimales + axis.xaxis.set_major_formatter(FuncFormatter(ticks_label_format)) + + axis.set_xlim(t, t+dt if dt > 0 else x[-1]) + axis.set_ylim(-1.1, 1.1) + + # resaltado de parte de la señal (solo si a != 0) + axis.axvspan(a, a+da, color='skyblue', + alpha=0 if a == 0 else 0.50, + label=f"Un periodo T={da}s" if da != 0 else "") + axis.legend(loc='upper left') + + if show: + plt.show() + + return figure, axis + +def normalize(data): + data = data.astype(np.float32) + data /= np.max(np.abs(data)) + return data + +def time_plot(fs, data, save_name="", t=0, dt=0, a=0, da=0): + show = True if save_name == "" else False + + # normaliza la amplitud dividiendo por el valor maximo del tipo de dato + data = normalize(data) + + x = np.arange(len(data)) / fs + fig, ax = time_graph_data(x, data, t, dt, a, da, show) + + if show == False: + save_plot(fig, save_name) + + return fig, ax + +def save_plot(fig, name): + base_name = os.path.basename(name) + file_name, ext = os.path.splitext(base_name) + file_path_no_ext = f'{plot_dir}{file_name}' + + save_name = f'{file_path_no_ext}.png' + print(save_name) + + # crea carpeta para plots + os.makedirs(plot_dir, exist_ok=True) + fig.savefig(save_name, dpi=250, bbox_inches="tight") + plt.close(fig) # liberar memoria + +def save_to_wav(fs, data, save_name): + # normalizar para prevenir clipping + data = data / np.max(np.abs(data)) + + # convertir a 16-bit PCM para WAV + data_as_int16 = np.int16(data * 32767) + + # crea carpeta para wavs + os.makedirs(out_dir, exist_ok=True) + + file_path = f'{out_dir}{save_name}' + print(file_path) + wavfile.write(file_path, fs, data_as_int16) + +# frecuencia + +# data = [[fft], [freqs], [legends]] +def freq_graph_multiple_data(data, f_min=0, f_max=0, y_min=0, y_max=0, show=True): + fig, axis = plt.subplots(figsize=(8, 4)) + + for i, (fft, freqs, label) in enumerate(data): + # print(label) + N = len(freqs) + x = freqs[:N // 2] + y = np.abs(fft[:N // 2]) + axis.plot(x, y, label=label, alpha=0.90, linewidth=((len(data)-i-1)*0.5 + 1.5)) + + axis.set(xlabel='Frecuencia [Hz]', ylabel='Magnitud') + + axis.minorticks_on() + axis.grid(True, which='major', color='black', linestyle=':', linewidth=1.00) + axis.grid(True, which='minor', color='black', linestyle=':', linewidth=0.50) + + # configuracion de ticks del eje x + axis.xaxis.set_major_locator(MaxNLocator(nbins=5)) + axis.xaxis.set_minor_locator(AutoMinorLocator(5)) + + axis.yaxis.set_major_locator(MaxNLocator(nbins=5)) + axis.yaxis.set_minor_locator(AutoMinorLocator(4)) + + plt.tight_layout() + + # max 3 decimales + axis.xaxis.set_major_formatter(FuncFormatter(ticks_label_format)) + plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0)) + + axis.set_xlim([f_min, f_max if f_max != 0 else 20000]) + axis.set_ylim([y_min, y_max if y_max != 0 else 1.05*max(y)]) + + axis.legend(loc='upper right') + + if show: + plt.show() + + return fig, axis + + +def freq_graph_data(x, y, f_min=0, f_max=0, y_min=0, y_max=0, show=True): + fig, axis = plt.subplots(figsize=(8, 4)) + + axis.plot(x, y) + axis.set(xlabel='Frecuencia [Hz]', ylabel='Magnitud') + + axis.minorticks_on() + axis.grid(True, which='major', color='black', linestyle=':', linewidth=1.00) + axis.grid(True, which='minor', color='black', linestyle=':', linewidth=0.50) + + # configuracion de ticks del eje x + axis.xaxis.set_major_locator(MaxNLocator(nbins=5)) + axis.xaxis.set_minor_locator(AutoMinorLocator(5)) + + axis.yaxis.set_major_locator(MaxNLocator(nbins=5)) + axis.yaxis.set_minor_locator(AutoMinorLocator(4)) + + plt.tight_layout() + + # max 3 decimales + axis.xaxis.set_major_formatter(FuncFormatter(ticks_label_format)) + plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0)) + + axis.set_xlim([f_min, f_max if f_max != 0 else 20000]) + axis.set_ylim([y_min, y_max if y_max != 0 else 1.05*max(y)]) + + if show: + plt.show() + + return fig, axis + +def freq_compute_fft(fs, data, t=0, dt=0, N=0): + i = 0 + di = fs*len(data) + if t != 0 or dt != 0: + i = int(t*fs) + di = int((t+dt)*fs) + + interval_data = data[i:di] + + # puntos de la fft + if N == 0: + N = len(interval_data) + + interval_fft = fft(interval_data, N) + interval_freqs = fftfreq(N, d=1/fs) + + return interval_fft, interval_freqs + +# hace la transformacion a frecuencias y pasa lo transformado a `freq_graph_data` +def freq_plot(fs, data, save_name="", f_min=0, f_max=0, y_min=0, y_max=0, + t=0, dt=0, a=0, da=0, show=False): + + interval_fft, interval_freqs = freq_compute_fft(fs, data, t, dt) + N = len(interval_fft) + + # se toma la parte positiva en ambos casos (primer parte del arreglo) + x = interval_freqs[:N // 2] + y = np.abs(interval_fft[:N // 2]) + + fig, ax = freq_graph_data(x, y, f_min, f_max, y_min, y_max, show=show) + + if save_name != "": + save_plot(fig, save_name) + + return fig, ax + +# frecuencia de muestreo comun +# computa y grafica en una figura la fft the los datos en `data_arr` +def freq_plot_multiple(fs, data_arr, leg_arr, save_name="", + f_min=0, f_max=0, y_min=0, y_max=0, t=0, dt=0, show=True): + + fft_freqs_arr = [] + for i, data in enumerate(data_arr): + fft, freqs = freq_compute_fft(fs, data, t, dt) + fft_freqs_arr.append([fft, freqs, leg_arr[i]]) + + fig, axis = freq_graph_multiple_data(fft_freqs_arr, f_min, f_max, y_min, y_max, show) + + if save_name != "": + save_plot(fig, save_name) + +def spectogram_plot(fs, data, save_name="", t=0, dt=0, N=1024, overlp=16, win='hamm', xlim=[], ylim=[], show=False): + if dt == 0: + dt = (len(data)/fs)-t + + i = int(t*fs) + di = int((t+dt)*fs) + interval_data = data[i:di] + + # `nperseg` tamaño de ventana (número de muestras por segmento) + # `noverlap` cantidad de solapamiento entre ventanas + f, time, Sxx = spectrogram(interval_data, fs=fs, nperseg=N, noverlap=overlp, + window=win) + fig, axis = plt.subplots(figsize=(8, 4)) + + # plt.pcolormesh(time, f, Sxx**0.10, shading='gouraud') + plt.pcolormesh(time, f, 10*np.log10(Sxx + 1e-12), shading='gouraud') + + plt.ylabel('Frecuencia [Hz]') + plt.xlabel('Tiempo [s]') + + if len(xlim) != 0: + plt.xlim(xlim) + + if len(ylim) != 0: + plt.ylim(ylim) + else: + plt.ylim(1, 20000) + + if show == True: + plt.show() + else: + if save_name != "": + save_plot(fig, save_name) + + return fig, axis + +def bode_plot(w, H, show=True): + figure, axis = plt.subplots(figsize=(8, 4)) + + axis.plot(w, 20*np.log10(np.abs(H))) + + axis.set(xlabel='Frecuencia [Hz]', ylabel='Magnitud [dB]') + axis.minorticks_on() + axis.grid(True, which='major', color='black', linestyle=':', linewidth=1.00) + axis.grid(True, which='minor', color='black', linestyle=':', linewidth=0.50) + plt.tight_layout() + + axis.set_xlim(0.0, 20e3) + + axis.legend(loc='upper left') + + if show: + plt.show() + + return figure, axis + +def freq_response_plot(w, H, phase, show=True, fc=20e3): + fig, ax1 = plt.subplots(figsize=(8, 4)) + + H_db = 20*np.log10(np.abs(H)) + + line1, = ax1.plot(w, H_db) + ax1.set(xlabel='Frecuencia [Hz]', ylabel='Magnitud [dB]') + ax1.minorticks_on() + ax1.grid(True, which='major', color='black', linestyle=':', linewidth=1.00) + ax1.grid(True, which='minor', color='black', linestyle=':', linewidth=0.50) + ax1.set_xlim(0.0, fc) + + ax2 = ax1.twinx() + line2, = ax2.plot(w, phase, color="tab:red") + ax2.set_ylabel("Fase [grados]", color="black") + ax2.tick_params(axis='y', labelcolor="black") + + # axis.yaxis.set_major_locator(MaxNLocator(nbins=5)) + ax1.yaxis.set_minor_locator(AutoMinorLocator(2)) + + # Add ONE point + + # Find index closest to -3 dB + idx = np.argmin(np.abs(H_db + 3)) # H_db = -3 => H_db +3 = 0 + w_3db = w[idx] + H_3db = H_db[idx] + line3 = ax1.scatter(w_3db, H_3db, color='tab:green', s=50, zorder=10) + + # nyquist + w_nyquist = 2756.25 + idx = np.argmin(np.abs(w - w_nyquist)) # H_db = -3 => H_db +3 = 0 + H_nyquist = H_db[idx] + line4 = ax1.scatter(w_nyquist, H_nyquist, color='tab:orange', s=50, zorder=10) + + ax1.legend([line1, line2, line3, line4], + ["Magnitud [dB]", + "Fase [grados]", + r'-3dB $\approx$ %0.0f Hz'%w_3db, + r'Nyquist $\approx$ %0.0f Hz'%w_nyquist], + loc='upper right') + + plt.tight_layout() + + if show: + plt.show() + + return fig, ax1, ax2 + +def dtime_plot(N, f, save_name="", legend="", n=0, dn=0, a=0, da=0): + show = True if save_name == "" else False + + n = np.arange(N + 1) + + fig, axis = plt.subplots(figsize=(8,4)) + axis.set(xlabel='Tiempo discreto', ylabel='Amplitud') + + markerline, stemlines, baseline = axis.stem( + n, f, + markerfmt='o', # tipo de marcador en la cabeza + basefmt="k-", + ) + + markerline.set_markersize(2.0) + stemlines.set_linewidth(0.35) + baseline.set_linewidth(0.5) + + stemlines.set_zorder(2) + markerline.set_zorder(3) + baseline.set_zorder(1) + + axis.grid(True, which='major', color='black', linestyle=':', linewidth=1.00) + axis.grid(True, which='minor', color='black', linestyle=':', linewidth=0.50) + axis.set_xlim(0, N+1) + axis.set_ylim(-0.03, 0.15) + + # configuracion de ticks del eje x + axis.xaxis.set_major_locator(MaxNLocator(nbins=15)) + axis.xaxis.set_minor_locator(AutoMinorLocator(2)) + + axis.yaxis.set_minor_locator(AutoMinorLocator(2)) + + # axis.yaxis.set_major_locator(MaxNLocator(nbins=5)) + # axis.yaxis.set_minor_locator(AutoMinorLocator(4)) + + if legend != "": + axis.legend([markerline], [legend], loc='upper right') + + if show == False: + save_plot(fig, save_name) + else: + plt.show() + + return fig, axis + +# np.linspace(start, stop, num).astype(int) diff --git a/tp/utils.py b/tp/utils.py @@ -1,432 +0,0 @@ -import matplotlib.pyplot as plt -from matplotlib.ticker import AutoMinorLocator -from matplotlib.ticker import MultipleLocator -from matplotlib.ticker import MaxNLocator -from matplotlib.ticker import FuncFormatter -from cycler import cycler - -import matplotlib -import numpy as np -import os - -from scipy.io import wavfile -from scipy.fft import fft, ifft, fftfreq -from scipy.signal import spectrogram - -plot_dir_name = 'plot' -out_dir_name = 'out' - -matplotlib.rcParams['font.family'] = 'Inter' -matplotlib.rcParams['font.size'] = 12 -matplotlib.rcParams['axes.prop_cycle'] = cycler( - color=['#1f77b4', '#ff0000', '#ff5f1f', 'green']) -matplotlib.use("TkAgg") - -def ticks_label_format(x, pos): - # 3 decimales, se eliminan los ceros y puntos - return f"{x:.3f}".rstrip("0").rstrip(".") - -def time_graph_multiple_data(x, y_arr, y_lab, t=0, dt=0, a=0, da=0, show=True): - figure, axis = plt.subplots(figsize=(8, 4)) - - for i, y in enumerate(y_arr): - axis.plot(x, y, label=y_lab[i], alpha=0.75) - - axis.set(xlabel='Tiempo [s]', ylabel='Amplitud normalizada') - - axis.minorticks_on() - axis.grid(True, which='major', color='black', linestyle=':', linewidth=1.00) - axis.grid(True, which='minor', color='black', linestyle=':', linewidth=0.50) - - # configuracion de ticks del eje x - axis.xaxis.set_major_locator(MaxNLocator(nbins=5)) - axis.xaxis.set_minor_locator(AutoMinorLocator(5)) - - axis.yaxis.set_major_locator(MaxNLocator(nbins=5)) - axis.yaxis.set_minor_locator(AutoMinorLocator(4)) - - plt.tight_layout() - - # max 3 decimales - axis.xaxis.set_major_formatter(FuncFormatter(ticks_label_format)) - - axis.set_xlim([t, t+dt if dt > 0 else x[-1]]) - axis.set_ylim(-1.1, 1.1) - - # resaltado de parte de la señal (solo si a != 0) - axis.axvspan(a, a+da, color='skyblue', - alpha=0 if a == 0 else 0.50, - label=f"Un periodo T={da}s" if da != 0 else "") - axis.legend(loc='upper left') - - if show: - plt.show() - - return figure, axis - -# todos deben la misma cantidad de elementos que el primero -def time_plot_multiple(fs, data_arr, leg_arr, save_name="", t=0, dt=0, a=0, da=0, show=False): - x = np.arange(len(data_arr[0])) / fs - fig, ax = time_graph_multiple_data(x, data_arr, leg_arr, t, dt, a=a, da=da, show=show) - - if show == False: - save_plot(fig, save_name) - - return fig, ax - -def time_graph_data(x, y, t=0, dt=0, a=0, da=0, show=True): - figure, axis = plt.subplots(figsize=(8, 4)) - - axis.plot(x, y, label='Señal de audio') - axis.set(xlabel='Tiempo [s]', ylabel='Amplitud normalizada') - - axis.minorticks_on() - axis.grid(True, which='major', color='black', linestyle=':', linewidth=1.00) - axis.grid(True, which='minor', color='black', linestyle=':', linewidth=0.50) - - # configuracion de ticks del eje x - axis.xaxis.set_major_locator(MaxNLocator(nbins=5)) - axis.xaxis.set_minor_locator(AutoMinorLocator(5)) - - axis.yaxis.set_major_locator(MaxNLocator(nbins=5)) - axis.yaxis.set_minor_locator(AutoMinorLocator(4)) - - plt.tight_layout() - - # max 3 decimales - axis.xaxis.set_major_formatter(FuncFormatter(ticks_label_format)) - - axis.set_xlim(t, t+dt if dt > 0 else x[-1]) - axis.set_ylim(-1.1, 1.1) - - # resaltado de parte de la señal (solo si a != 0) - axis.axvspan(a, a+da, color='skyblue', - alpha=0 if a == 0 else 0.50, - label=f"Un periodo T={da}s" if da != 0 else "") - axis.legend(loc='upper left') - - if show: - plt.show() - - return figure, axis - -def normalize(data): - data = data.astype(np.float32) - data /= np.max(np.abs(data)) - return data - -def time_plot(fs, data, save_name="", t=0, dt=0, a=0, da=0): - show = True if save_name == "" else False - - # normaliza la amplitud dividiendo por el valor maximo del tipo de dato - data = normalize(data) - - x = np.arange(len(data)) / fs - fig, ax = time_graph_data(x, data, t, dt, a, da, show) - - if show == False: - save_plot(fig, save_name) - - return fig, ax - -def save_plot(fig, name): - base_name = os.path.basename(name) - file_name, ext = os.path.splitext(base_name) - file_path_no_ext = f'{plot_dir_name}/{file_name}' - - save_name = f'{file_path_no_ext}.png' - print(save_name) - - # crea carpeta para plots - os.makedirs(plot_dir_name, exist_ok=True) - fig.savefig(save_name, dpi=250, bbox_inches="tight") - plt.close(fig) # liberar memoria - -def save_to_wav(fs, data, save_name): - # normalizar para prevenir clipping - data = data / np.max(np.abs(data)) - - # convertir a 16-bit PCM para WAV - data_as_int16 = np.int16(data * 32767) - - # crea carpeta para wavs - os.makedirs(out_dir_name, exist_ok=True) - - file_path = f'{out_dir_name}/{save_name}' - print(file_path) - wavfile.write(file_path, fs, data_as_int16) - -# frecuencia - -# data = [[fft], [freqs], [legends]] -def freq_graph_multiple_data(data, f_min=0, f_max=0, y_min=0, y_max=0, show=True): - fig, axis = plt.subplots(figsize=(8, 4)) - - for i, (fft, freqs, label) in enumerate(data): - # print(label) - N = len(freqs) - x = freqs[:N // 2] - y = np.abs(fft[:N // 2]) - axis.plot(x, y, label=label, alpha=0.90, linewidth=((len(data)-i-1)*0.5 + 1.5)) - - axis.set(xlabel='Frecuencia [Hz]', ylabel='Magnitud') - - axis.minorticks_on() - axis.grid(True, which='major', color='black', linestyle=':', linewidth=1.00) - axis.grid(True, which='minor', color='black', linestyle=':', linewidth=0.50) - - # configuracion de ticks del eje x - axis.xaxis.set_major_locator(MaxNLocator(nbins=5)) - axis.xaxis.set_minor_locator(AutoMinorLocator(5)) - - axis.yaxis.set_major_locator(MaxNLocator(nbins=5)) - axis.yaxis.set_minor_locator(AutoMinorLocator(4)) - - plt.tight_layout() - - # max 3 decimales - axis.xaxis.set_major_formatter(FuncFormatter(ticks_label_format)) - plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0)) - - axis.set_xlim([f_min, f_max if f_max != 0 else 20000]) - axis.set_ylim([y_min, y_max if y_max != 0 else 1.05*max(y)]) - - axis.legend(loc='upper right') - - if show: - plt.show() - - return fig, axis - - -def freq_graph_data(x, y, f_min=0, f_max=0, y_min=0, y_max=0, show=True): - fig, axis = plt.subplots(figsize=(8, 4)) - - axis.plot(x, y) - axis.set(xlabel='Frecuencia [Hz]', ylabel='Magnitud') - - axis.minorticks_on() - axis.grid(True, which='major', color='black', linestyle=':', linewidth=1.00) - axis.grid(True, which='minor', color='black', linestyle=':', linewidth=0.50) - - # configuracion de ticks del eje x - axis.xaxis.set_major_locator(MaxNLocator(nbins=5)) - axis.xaxis.set_minor_locator(AutoMinorLocator(5)) - - axis.yaxis.set_major_locator(MaxNLocator(nbins=5)) - axis.yaxis.set_minor_locator(AutoMinorLocator(4)) - - plt.tight_layout() - - # max 3 decimales - axis.xaxis.set_major_formatter(FuncFormatter(ticks_label_format)) - plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0)) - - axis.set_xlim([f_min, f_max if f_max != 0 else 20000]) - axis.set_ylim([y_min, y_max if y_max != 0 else 1.05*max(y)]) - - if show: - plt.show() - - return fig, axis - -def freq_compute_fft(fs, data, t=0, dt=0): - i = 0 - di = fs*len(data) - if t != 0 or dt != 0: - i = int(t*fs) - di = int((t+dt)*fs) - - interval_data = data[i:di] - - # puntos de la fft - N = len(interval_data) - interval_fft = fft(interval_data, 20000) - interval_freqs = fftfreq(20000, d=1/fs) - - return interval_fft, interval_freqs - -# hace la transformacion a frecuencias y pasa lo transformado a `freq_graph_data` -def freq_plot(fs, data, save_name="", f_min=0, f_max=0, y_min=0, y_max=0, - t=0, dt=0, a=0, da=0, show=False): - - interval_fft, interval_freqs = freq_compute_fft(fs, data, t, dt) - N = len(interval_fft) - - # se toma la parte positiva en ambos casos (primer parte del arreglo) - x = interval_freqs[:N // 2] - y = np.abs(interval_fft[:N // 2]) - - fig, ax = freq_graph_data(x, y, f_min, f_max, y_min, y_max, show=show) - - if save_name != "": - save_plot(fig, save_name) - - return fig, ax - -# frecuencia de muestreo comun -# computa y grafica en una figura la fft the los datos en `data_arr` -def freq_plot_multiple(fs, data_arr, leg_arr, save_name="", - f_min=0, f_max=0, y_min=0, y_max=0, t=0, dt=0, show=True): - - fft_freqs_arr = [] - for i, data in enumerate(data_arr): - fft, freqs = freq_compute_fft(fs, data, t, dt) - fft_freqs_arr.append([fft, freqs, leg_arr[i]]) - - fig, axis = freq_graph_multiple_data(fft_freqs_arr, f_min, f_max, y_min, y_max, show) - - if save_name != "": - save_plot(fig, save_name) - -def spectogram_plot(fs, data, save_name="", t=0, dt=0, N=1024, overlp=16, win='hamm', xlim=[], ylim=[], show=False): - if dt == 0: - dt = (len(data)/fs)-t - - i = int(t*fs) - di = int((t+dt)*fs) - interval_data = data[i:di] - - # `nperseg` tamaño de ventana (número de muestras por segmento) - # `noverlap` cantidad de solapamiento entre ventanas - f, time, Sxx = spectrogram(interval_data, fs=fs, nperseg=N, noverlap=overlp, - window=win) - fig, axis = plt.subplots(figsize=(8, 4)) - - # plt.pcolormesh(time, f, Sxx**0.10, shading='gouraud') - plt.pcolormesh(time, f, 10*np.log10(Sxx + 1e-12), shading='gouraud') - - plt.ylabel('Frecuencia [Hz]') - plt.xlabel('Tiempo [s]') - - if len(xlim) != 0: - plt.xlim(xlim) - - if len(ylim) != 0: - plt.ylim(ylim) - else: - plt.ylim(1, 20000) - - if show == True: - plt.show() - else: - if save_name != "": - save_plot(fig, save_name) - - return fig, axis - -def bode_plot(w, H, show=True): - figure, axis = plt.subplots(figsize=(8, 4)) - - axis.plot(w, 20*np.log10(np.abs(H))) - - axis.set(xlabel='Frecuencia [Hz]', ylabel='Magnitud [dB]') - axis.minorticks_on() - axis.grid(True, which='major', color='black', linestyle=':', linewidth=1.00) - axis.grid(True, which='minor', color='black', linestyle=':', linewidth=0.50) - plt.tight_layout() - - axis.set_xlim(0.0, 20e3) - - axis.legend(loc='upper left') - - if show: - plt.show() - - return figure, axis - -def freq_response_plot(w, H, phase, show=True, fc=20e3): - fig, ax1 = plt.subplots(figsize=(8, 4)) - - H_db = 20*np.log10(np.abs(H)) - - line1, = ax1.plot(w, H_db) - ax1.set(xlabel='Frecuencia [Hz]', ylabel='Magnitud [dB]') - ax1.minorticks_on() - ax1.grid(True, which='major', color='black', linestyle=':', linewidth=1.00) - ax1.grid(True, which='minor', color='black', linestyle=':', linewidth=0.50) - ax1.set_xlim(0.0, fc) - - ax2 = ax1.twinx() - line2, = ax2.plot(w, phase, color="tab:red") - ax2.set_ylabel("Fase [grados]", color="black") - ax2.tick_params(axis='y', labelcolor="black") - - # axis.yaxis.set_major_locator(MaxNLocator(nbins=5)) - ax1.yaxis.set_minor_locator(AutoMinorLocator(2)) - - # Add ONE point - - # Find index closest to -3 dB - idx = np.argmin(np.abs(H_db + 3)) # H_db = -3 => H_db +3 = 0 - w_3db = w[idx] - H_3db = H_db[idx] - line3 = ax1.scatter(w_3db, H_3db, color='tab:green', s=50, zorder=10) - - # nyquist - w_nyquist = 2756.25 - idx = np.argmin(np.abs(w - w_nyquist)) # H_db = -3 => H_db +3 = 0 - H_nyquist = H_db[idx] - line4 = ax1.scatter(w_nyquist, H_nyquist, color='tab:orange', s=50, zorder=10) - - ax1.legend([line1, line2, line3, line4], - ["Magnitud [dB]", - "Fase [grados]", - r'-3dB $\approx$ %0.0f Hz'%w_3db, - r'Nyquist $\approx$ %0.0f Hz'%w_nyquist], - loc='upper right') - - plt.tight_layout() - - if show: - plt.show() - - return fig, ax1, ax2 - -def dtime_plot(N, f, save_name="", legend="", n=0, dn=0, a=0, da=0): - show = True if save_name == "" else False - - n = np.arange(N + 1) - - fig, axis = plt.subplots(figsize=(8,4)) - axis.set(xlabel='Tiempo discreto', ylabel='Amplitud') - - markerline, stemlines, baseline = axis.stem( - n, f, - markerfmt='o', # tipo de marcador en la cabeza - basefmt="k-", - ) - - markerline.set_markersize(2.0) - stemlines.set_linewidth(0.35) - baseline.set_linewidth(0.5) - - stemlines.set_zorder(2) - markerline.set_zorder(3) - baseline.set_zorder(1) - - axis.grid(True, which='major', color='black', linestyle=':', linewidth=1.00) - axis.grid(True, which='minor', color='black', linestyle=':', linewidth=0.50) - axis.set_xlim(0, N+1) - axis.set_ylim(-0.03, 0.15) - - # configuracion de ticks del eje x - axis.xaxis.set_major_locator(MaxNLocator(nbins=15)) - axis.xaxis.set_minor_locator(AutoMinorLocator(2)) - - axis.yaxis.set_minor_locator(AutoMinorLocator(2)) - - # axis.yaxis.set_major_locator(MaxNLocator(nbins=5)) - # axis.yaxis.set_minor_locator(AutoMinorLocator(4)) - - if legend != "": - axis.legend([markerline], [legend], loc='upper right') - - if show == False: - save_plot(fig, save_name) - else: - plt.show() - - return fig, axis - -# np.linspace(start, stop, num).astype(int)
