You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
751 lines
60 KiB
Plaintext
751 lines
60 KiB
Plaintext
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"# Generacion de muestras"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"import numpy as np\n",
|
|
"\n",
|
|
"def generar_muestra_coseno(N, num_cosenos, varianza_ruido):\n",
|
|
" \"\"\"\n",
|
|
" Genera un vector de muestras como suma de cosenos con parámetros aleatorios y ruido gaussiano.\n",
|
|
" \n",
|
|
" Parámetros:\n",
|
|
" -----------\n",
|
|
" N : int\n",
|
|
" Cantidad de muestras a generar.\n",
|
|
" num_cosenos : int\n",
|
|
" Número de cosenos a sumar.\n",
|
|
" varianza_ruido : float\n",
|
|
" Varianza del ruido gaussiano agregado.\n",
|
|
" \n",
|
|
" Retorna:\n",
|
|
" --------\n",
|
|
" x : ndarray de shape (N,)\n",
|
|
" Vector de muestras resultante.\n",
|
|
" parametros : list of dict\n",
|
|
" Lista con las constantes de cada coseno (amplitud, frecuencia y fase).\n",
|
|
" varianza_ruido : float\n",
|
|
" Varianza del ruido utilizado.\n",
|
|
" \"\"\"\n",
|
|
" t = np.arange(N) / N # tiempo normalizado entre 0 y 1\n",
|
|
" \n",
|
|
" # Parámetros aleatorios de los cosenos\n",
|
|
" amplitudes = np.random.uniform(0.5, 2.0, num_cosenos) # amplitudes aleatorias\n",
|
|
" frecuencias = np.random.uniform(1, 20, num_cosenos) # frecuencias aleatorias (en Hz normalizado)\n",
|
|
" fases = np.random.uniform(0, 2*np.pi, num_cosenos) # fases aleatorias\n",
|
|
" \n",
|
|
" # Generar señal base como suma de cosenos\n",
|
|
" x_clean = np.zeros(N)\n",
|
|
" for A, f, phi in zip(amplitudes, frecuencias, fases):\n",
|
|
" x_clean += A * np.cos(2 * np.pi * f * t + phi)\n",
|
|
" \n",
|
|
" # Agregar ruido gaussiano\n",
|
|
" ruido = np.random.normal(0, np.sqrt(varianza_ruido), N)\n",
|
|
" x = x_clean + ruido\n",
|
|
" \n",
|
|
" # Empaquetar parámetros\n",
|
|
" parametros = []\n",
|
|
" for A, f, phi in zip(amplitudes, frecuencias, fases):\n",
|
|
" parametros.append({\n",
|
|
" \"amplitud\": A,\n",
|
|
" \"frecuencia\": f,\n",
|
|
" \"fase\": phi\n",
|
|
" })\n",
|
|
" \n",
|
|
" return x, parametros, varianza_ruido\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"# Estimador de "
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"from scipy.optimize import curve_fit\n",
|
|
"\n",
|
|
"def modelo_cosenos(t, *params):\n",
|
|
" \"\"\"\n",
|
|
" Modelo de suma de cosenos.\n",
|
|
" \n",
|
|
" params = [A1, f1, phi1, A2, f2, phi2, ..., An, fn, phin]\n",
|
|
" \"\"\"\n",
|
|
" num_cosenos = len(params) // 3\n",
|
|
" y = np.zeros_like(t)\n",
|
|
" for i in range(num_cosenos):\n",
|
|
" A = params[3*i]\n",
|
|
" f = params[3*i + 1]\n",
|
|
" phi = params[3*i + 2]\n",
|
|
" y += A * np.cos(2 * np.pi * f * t + phi)\n",
|
|
" return y\n",
|
|
"\n",
|
|
"def estimar_parametros(y, num_cosenos, max_iter=10000):\n",
|
|
" \"\"\"\n",
|
|
" Estima los parámetros de una señal que es suma de cosenos mediante mínimos cuadrados no lineales.\n",
|
|
" \n",
|
|
" Parámetros:\n",
|
|
" -----------\n",
|
|
" y : ndarray\n",
|
|
" Vector de muestras observadas.\n",
|
|
" num_cosenos : int\n",
|
|
" Número de cosenos que componen la señal.\n",
|
|
" max_iter : int\n",
|
|
" Número máximo de iteraciones del optimizador.\n",
|
|
" \n",
|
|
" Retorna:\n",
|
|
" --------\n",
|
|
" parametros : list of dict\n",
|
|
" Estimación de los parámetros (amplitud, frecuencia, fase).\n",
|
|
" \"\"\"\n",
|
|
" N = len(y)\n",
|
|
" t = np.arange(N) / N\n",
|
|
" \n",
|
|
" # Valores iniciales aproximados (necesarios para que converja)\n",
|
|
" A0 = np.ones(num_cosenos) * (np.std(y) / num_cosenos)\n",
|
|
" f0 = np.linspace(1, 10, num_cosenos) # inicializar frecuencias en un rango\n",
|
|
" phi0 = np.zeros(num_cosenos)\n",
|
|
" \n",
|
|
" p0 = []\n",
|
|
" for A, f, phi in zip(A0, f0, phi0):\n",
|
|
" p0 += [A, f, phi]\n",
|
|
" p0 = np.array(p0)\n",
|
|
" \n",
|
|
" # Ajuste no lineal\n",
|
|
" popt, _ = curve_fit(modelo_cosenos, t, y, p0=p0, maxfev=max_iter)\n",
|
|
" \n",
|
|
" # Empaquetar parámetros\n",
|
|
" parametros = []\n",
|
|
" for i in range(num_cosenos):\n",
|
|
" A = popt[3*i]\n",
|
|
" f = popt[3*i+1]\n",
|
|
" phi = popt[3*i+2]\n",
|
|
" parametros.append({\n",
|
|
" \"amplitud\": A,\n",
|
|
" \"frecuencia\": f,\n",
|
|
" \"fase\": phi\n",
|
|
" })\n",
|
|
" \n",
|
|
" return parametros\n",
|
|
"\n",
|
|
"def reconstruir_senal(N, parametros):\n",
|
|
" \"\"\"\n",
|
|
" Reconstruye la señal a partir de parámetros de cosenos.\n",
|
|
" \n",
|
|
" parametros : list of dict con {\"amplitud\", \"frecuencia\", \"fase\"}\n",
|
|
" \"\"\"\n",
|
|
" t = np.arange(N) / N\n",
|
|
" y = np.zeros(N)\n",
|
|
" for p in parametros:\n",
|
|
" y += p[\"amplitud\"] * np.cos(2 * np.pi * p[\"frecuencia\"] * t + p[\"fase\"])\n",
|
|
" return t, y\n",
|
|
"\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"\n",
|
|
"def estimador_de_orden(y, orden_max, C, max_iter=10000, epsilon=1e-12, devolver_parametros=False):\n",
|
|
" \"\"\"\n",
|
|
" Estima el orden óptimo (número de cosenos) usando un criterio:\n",
|
|
" Criterio(k) = N * log(var_residuo_k + epsilon) + C * k\n",
|
|
"\n",
|
|
" Parámetros\n",
|
|
" ----------\n",
|
|
" y : ndarray (N,)\n",
|
|
" Señal observada.\n",
|
|
" orden_max : int\n",
|
|
" Orden máximo a evaluar (número máximo de cosenos).\n",
|
|
" C : float\n",
|
|
" Constante que multiplica al término de penalización por el orden k.\n",
|
|
" max_iter : int\n",
|
|
" Máx. iteraciones para el optimizador en estimar_parametros.\n",
|
|
" epsilon : float\n",
|
|
" Pequeño valor para estabilizar el logaritmo.\n",
|
|
" devolver_parametros : bool\n",
|
|
" Si True, devuelve también los parámetros del orden óptimo.\n",
|
|
"\n",
|
|
" Retorna\n",
|
|
" -------\n",
|
|
" orden_optimo : int\n",
|
|
" Orden estimado.\n",
|
|
" varianzas : list[float]\n",
|
|
" Varianzas del residuo para k = 1..orden_max.\n",
|
|
" criterios : list[float]\n",
|
|
" Valores del criterio para k = 1..orden_max.\n",
|
|
" (opcional) parametros_optimos : list[dict]\n",
|
|
" Parámetros {amplitud, frecuencia, fase} del orden óptimo (si devolver_parametros=True).\n",
|
|
" \"\"\"\n",
|
|
" N = len(y)\n",
|
|
" varianzas = []\n",
|
|
" criterios = []\n",
|
|
" mejores_parametros = None\n",
|
|
" parametros_k_opt = None\n",
|
|
"\n",
|
|
" for k in range(1, orden_max + 1):\n",
|
|
" try:\n",
|
|
" # Estimar parámetros para orden k\n",
|
|
" params_k = estimar_parametros(y, num_cosenos=k, max_iter=max_iter)\n",
|
|
"\n",
|
|
" # Reconstruir y calcular residuo\n",
|
|
" _, y_hat = reconstruir_senal(N, params_k)\n",
|
|
" residuo = y - y_hat\n",
|
|
" var_est = float(np.mean(residuo**2))\n",
|
|
"\n",
|
|
" # Guardar métricas\n",
|
|
" varianzas.append(var_est)\n",
|
|
" criterio_k = N * np.log(var_est + epsilon) + C * k\n",
|
|
" criterios.append(criterio_k)\n",
|
|
"\n",
|
|
" # Guardar los parámetros si van siendo los mejores\n",
|
|
" if devolver_parametros:\n",
|
|
" if (parametros_k_opt is None) or (criterio_k < criterios[parametros_k_opt - 1]):\n",
|
|
" parametros_k_opt = k\n",
|
|
" mejores_parametros = params_k\n",
|
|
"\n",
|
|
" except Exception:\n",
|
|
" # Si falla el ajuste para este orden, lo penalizamos fuertemente\n",
|
|
" varianzas.append(np.inf)\n",
|
|
" criterios.append(np.inf)\n",
|
|
"\n",
|
|
" # Selección del orden óptimo\n",
|
|
" orden_optimo = int(np.argmin(criterios) + 1)\n",
|
|
"\n",
|
|
" if devolver_parametros:\n",
|
|
" # Si no se setearon (por fallos), intentar reevaluar los del orden óptimo\n",
|
|
" if mejores_parametros is None or parametros_k_opt != orden_optimo:\n",
|
|
" try:\n",
|
|
" mejores_parametros = estimar_parametros(y, num_cosenos=orden_optimo, max_iter=max_iter)\n",
|
|
" except Exception:\n",
|
|
" mejores_parametros = None\n",
|
|
" return orden_optimo, varianzas, criterios, mejores_parametros\n",
|
|
"\n",
|
|
" return orden_optimo, varianzas, criterios\n",
|
|
"\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 5,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Orden real: 3\n",
|
|
"Orden estimado: 2\n",
|
|
"Varianzas: [2.323137726399167, 1.8383255140487047, 1.837929702038553, 1.8368502482272016, 1.8355685769211567, 1.8287296226777754, inf, inf]\n",
|
|
"Criterios: [np.float64(426.95937009636503), np.float64(315.4275552824338), np.float64(320.8198881079817), np.float64(326.02614158247326), np.float64(331.177142348153), np.float64(334.8107652653657), inf, inf]\n",
|
|
"Parámetros óptimos: [{'amplitud': np.float64(0.07807501614245879), 'frecuencia': np.float64(0.053583028276843), 'fase': np.float64(1.4717432359599454)}, {'amplitud': np.float64(0.9879446102229781), 'frecuencia': np.float64(9.435587905228454), 'fase': np.float64(0.37640582199551625)}]\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"# Señal sintética\n",
|
|
"x, params_true, _ = generar_muestra_coseno(N=500, num_cosenos=3, varianza_ruido=0.05)\n",
|
|
"\n",
|
|
"# Estimo orden con constante C=5.0 (ajustá C a tu gusto)\n",
|
|
"orden_opt, vars_k, crits_k, params_opt = estimador_de_orden(\n",
|
|
" x, orden_max=8, C=5.50, devolver_parametros=True\n",
|
|
")\n",
|
|
"\n",
|
|
"print(\"Orden real:\", len(params_true))\n",
|
|
"print(\"Orden estimado:\", orden_opt)\n",
|
|
"print(\"Varianzas:\", vars_k)\n",
|
|
"print(\"Criterios:\", crits_k)\n",
|
|
"print(\"Parámetros óptimos:\", params_opt)\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 6,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Orden real: 3\n",
|
|
"Tasa de acierto: 0.16666666666666666\n",
|
|
"Órdenes estimados: [4, 6, 7, 5, 3, 6, 3, 1, 5, 3, 5, 2, 4, 1, 4, 1, 6, 1, 4, 3, 6, 5, 4, 3, 8, 4, 6, 7, 4, 8]\n"
|
|
]
|
|
},
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAioAAAHHCAYAAACRAnNyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABQl0lEQVR4nO3deXxM5/4H8M9km+wbCUEWkiBiKYl9CRXSWGqppcq1VSnRWLqJtkiV0LqW1q2tt+iCopRaInbX2sRWW0liS20hSGSRZeb5/eGXU2MSYsw4x/i8X6/zcpZnzvmemTH5zrONSgghQERERKRAFnIHQERERFQaJipERESkWExUiIiISLGYqBAREZFiMVEhIiIixWKiQkRERIrFRIWIiIgUi4kKERERKRYTFXrh5OfnY+rUqdiyZYvcoRjF7du3ERsbi0OHDskdChGR4jBRUbBJkyZBpVI9l2u1bt0arVu3lrZ37doFlUqF1atXP5frP0ylUmHSpEmlHh87dix+/vlnNG7c+LnEM3DgQPj5+Znk3EII9O/fH7t27UL9+vWf+Xx+fn4YOHDgswf2Enj0Pa8EfP3Mw5M+w4o9z8/4FxkTledkyZIlUKlU0mJra4tKlSohIiICX3/9Ne7du2eU61y9ehWTJk3CsWPHjHI+pVm5ciV+++03bN68Ga6urnKH88y+/PJLXLx4EWvXroWNjY3c4Zid06dPY9KkSbh48aLcobzU8vLy8Pbbb6N27dpwcXGBo6Mj6tWrhzlz5qCwsFCv/N27dzF06FB4eHjAwcEBbdq0wZEjR2SInJTASu4AXjaff/45qlatisLCQly/fh27du3C6NGjMXPmTKxfvx5169aVyn766acYN27cU53/6tWriI2NhZ+fH1555ZUyPy4hIeGprmNKeXl5sLLSf2sKIfD3339j8+bN8PHxkSEy47p//z6KioqwadMms0i6lOj06dOIjY1F69at9WrFlPSeN3d5eXk4deoUOnToAD8/P1hYWGD//v0YM2YMDh06hGXLlklltVotOnbsiOPHj+PDDz9E+fLl8e2336J169Y4fPgwAgMDZbyTsintM4wMw2fyOYuMjERoaKi0HRMTgx07dqBTp054/fXXcebMGdjZ2QEArKysTP5mz83Nhb29vaK+zdva2pa4X6VSYezYsc85GtOxtbXFJ5988lSPycnJgYODg4kierko6T1v7tzd3XHw4EGdfe+++y5cXFwwd+5czJw5ExUrVgQArF69Gvv378eqVavQo0cPAECvXr1QvXp1TJw4USepkUNp/we1Wi0KCgpga2tb6mcYGYZNPwrw6quv4rPPPsOlS5fw008/SftLar/cunUrWrRoAVdXVzg6OqJGjRoYP348gAf9Sho2bAgAGDRokNTMtGTJEgAP2uRr166Nw4cPo1WrVrC3t5ceW1p7vUajwfjx41GxYkU4ODjg9ddfR1pamk6Z0trVSzrn/fv3MWnSJFSvXh22trbw8vJC9+7dkZqaKpUpqX336NGjiIyMhLOzMxwdHdG2bVu9D77i5rV9+/Zh7NixUrVxt27dcPPmTb34SvLbb7+hdu3asLW1Re3atbF27doSy2m1WsyePRvBwcGwtbVFhQoVMGzYMNy5c6dM19mxYwdatmwJBwcHuLq6okuXLjhz5oxOmeLX//Tp03jrrbfg5uaGFi1aAHhQu/TFF1+gSpUqsLe3R5s2bXDq1KkSr3X37l2MHj0a3t7eUKvVCAgIwPTp06HVaqUyFy9ehEqlwowZM7Bw4UL4+/tDrVajYcOGSExM1DvnX3/9hR49esDd3R22trYIDQ3F+vXrdcoUFhYiNjYWgYGBsLW1Rbly5dCiRQts3br1ic9PWWIGgBUrViAkJAROTk5wdnZGnTp1MGfOHAAP3g89e/YEALRp00b6/7Br1y4ApffLWrlyJWJjY1G5cmU4OTmhR48eyMzMRH5+PkaPHg1PT084Ojpi0KBByM/P14ln8eLFePXVV+Hp6Qm1Wo1atWph3rx5evf3NK/f+fPn0bNnT7i7u8Pe3h5NmjTBxo0b9cp98803CA4Ohr29Pdzc3BAaGlqmP+r5+fmYOHEiAgICoFar4e3tjY8++kjv3lQqFUaOHCn9H1Gr1QgODkZ8fPwTr1Ga4lquu3fvSvtWr16NChUqoHv37tI+Dw8P9OrVC+vWrdOLqySbN29GWFiY9L5o2LCh3nOxatUqhISEwM7ODuXLl0e/fv1w5coVnTIDBw6Eo6MjUlNT0aFDBzg5OaFv374A/nk+fv75ZwQHB0OtVkvPRUmfYXv37kXDhg1ha2sLf39/LFiwoMTYy/oeSkpKQkREBMqXLw87OztUrVoVgwcPfuJz8yJijYpC/Otf/8L48eORkJCAd955p8Qyp06dQqdOnVC3bl18/vnnUKvVSElJwb59+wAAQUFB+PzzzzFhwgQMHToULVu2BAA0a9ZMOkdGRgYiIyPx5ptvol+/fqhQocJj45oyZQpUKhU+/vhjpKenY/bs2QgPD8exY8ekmp+y0mg06NSpE7Zv344333wTo0aNwr1797B161acPHkS/v7+pd53y5Yt4ezsjI8++gjW1tZYsGABWrdujd27d+t1qn3vvffg5uaGiRMn4uLFi5g9ezZGjhyJX3755bHxJSQk4I033kCtWrUQFxeHjIwMDBo0CFWqVNErO2zYMCxZsgSDBg1CdHQ0Lly4gLlz5+Lo0aPYt28frK2tS73Otm3bEBkZiWrVqmHSpEnIy8vDN998g+bNm+PIkSN6TRQ9e/ZEYGAgpk6dCiEEAGDChAn44osv0KFDB3To0AFHjhxB+/btUVBQoPPY3NxchIWF4cqVKxg2bBh8fHywf/9+xMTE4Nq1a5g9e7ZO+WXLluHevXsYNmwYVCoVvvzyS3Tv3h3nz5+X7unUqVNo3rw5KleujHHjxsHBwQErV65E165d8euvv6Jbt24AHiRacXFxGDJkCBo1aoSsrCwkJSXhyJEjaNeuXanPT1lj3rp1K/r06YO2bdti+vTpAIAzZ85g3759GDVqFFq1aoXo6Gh8/fXXGD9+PIKCggBA+rc0cXFxsLOzw7hx45CSkoJvvvkG1tbWsLCwwJ07dzBp0iQcPHgQS5YsQdWqVTFhwgTpsfPmzUNwcDBef/11WFlZ4ffff8eIESOg1WoRFRUllSvr63fjxg00a9YMubm5iI6ORrly5bB06VK8/vrrWL16tfRcL1q0CNHR0ejRowdGjRqF+/fv488//8ShQ4fw1ltvlXqvWq0Wr7/+Ovbu3YuhQ4ciKCgIJ06cwKxZs3Du3Dn89ttvOuX37t2LNWvWYMSIEXBycsLXX3+NN954A5cvX0a5cuUe+7wCQEFBAbKyspCXl4ekpCTMmDEDvr6+CAgIkMocPXoUDRo0gIWF7vfoRo0aYeHChTh37hzq1KlT6jWWLFmCwYMHIzg4GDExMXB1dcXRo0cRHx8vPRfF/3cbNmyIuLg43LhxA3PmzMG+fftw9OhRnabYoqIiREREoEWLFpgxYwbs7e2lYzt27MDKlSsxcuRIlC9fvtRO9ydOnED79u3h4eGBSZMmoaioCBMnTizx87cs76H09HTpfOPGjYOrqysuXryINWvWPPE1eCEJei4WL14sAIjExMRSy7i4uIj69etL2xMnThQPv0SzZs0SAMTNmzdLPUdiYqIAIBYvXqx3LCwsTAAQ8+fPL/FYWFiYtL1z504BQFSuXFlkZWVJ+1euXCkAiDlz5kj7fH19xYABA554zu+//14AEDNnztQrq9VqpXUAYuLEidJ2165dhY2NjUhNTZX2Xb16VTg5OYlWrVpJ+4qf4/DwcJ3zjRkzRlhaWoq7d+/qXfdhr7zyivDy8tIpl5CQIAAIX19fad///vc/AUD8/PPPOo+Pj48vcX9J1/H09BQZGRnSvuPHjwsLCwvRv39/aV/x69+nTx+dx6enpwsbGxvRsWNHnfscP368AKDzWkyePFk4ODiIc+fO6Zxj3LhxwtLSUly+fFkIIcSFCxcEAFGuXDlx+/Ztqdy6desEAPH7779L+9q2bSvq1Kkj7t+/L+3TarWiWbNmIjAwUNpXr1490bFjx8c+FyUpa8yjRo0Szs7OoqioqNRzrVq1SgAQO3fu1DtW2nu+du3aoqCgQNrfp08foVKpRGRkpM7jmzZtqvO+EEKI3NxcvetERESIatWqSdtP8/qNHj1aABD/+9//pH337t0TVatWFX5+fkKj0QghhOjSpYsIDg4u9XkozY8//igsLCx0zi+EEPPnzxcAxL59+6R9AISNjY1ISUmR9h0/flwAEN98802Zrrd8+XIBQFpCQ0PFn3/+qVPGwcFBDB48WO+xGzduFABEfHx8qee/e/eucHJyEo0bNxZ5eXk6x4qf64KCAuHp6Slq166tU2bDhg0CgJgwYYK0b8CAAQKAGDdunN61AAgLCwtx6tSpEo89+hlma2srLl26JO07ffq0sLS01PmMF6Js76G1a9c+8e+JOWHTj4I4Ojo+dvRPcZa/bt06vSrwslKr1Rg0aFCZy/fv3x9OTk7Sdo8ePeDl5YVNmzY99bV//fVXlC9fHu+9957esdKG6Gk0GiQkJKBr166oVq2atN/LywtvvfUW9u7di6ysLJ3HDB06VOd8LVu2hEajwaVLl0qN7dq1azh27BgGDBgAFxcXaX+7du1Qq1YtnbKrVq2Ci4sL2rVrh1u3bklLSEgIHB0dsXPnzideZ+DAgXB3d5f2161bF+3atSvxeX333Xd1trdt24aCggK89957Ovc5evRovceuWrUKLVu2hJubm06s4eHh0Gg02LNnj0753r17w83NTdourpU7f/48gAdzvuzYsQO9evXCvXv3pPNlZGQgIiICycnJUvW5q6srTp06heTk5FKfj5KUNWZXV1fk5OSUqSnpafTv31+nRqxx48YQQuhVqzdu3BhpaWkoKiqS9j1cy5iZmYlbt24hLCwM58+fR2ZmJoCne/02bdqERo0aSU1+wIPPiaFDh+LixYs4ffo0gAfPxd9//11iM93jrFq1CkFBQahZs6bOc/3qq68CgN57OTw8XKfms27dunB2dpbeH0/Spk0bbN26FatWrcK7774La2tr5OTk6JTJy8uDWq3We2xxv4+8vLxSz79161bcu3cP48aN0+snUvxcJyUlIT09HSNGjNAp07FjR9SsWbPEZrXhw4eXeL2wsDC9z4dHaTQabNmyBV27dtUZBBAUFISIiAi98mV5DxX/LdiwYUOJo6bMDRMVBcnOztZJCh7Vu3dvNG/eHEOGDEGFChXw5ptvYuXKlU+VtFSuXPmpOhE+2sNepVIhICDAoOGeqampqFGjxlN1EL558yZyc3NRo0YNvWNBQUHQarV6fWYeHRFU/If3cf1HipOYkkYUPHrt5ORkZGZmwtPTEx4eHjpLdnY20tPTn3id0u7n1q1beh/cVatWLVOsHh4eOklGcazx8fF6cYaHhwOAXqxPeu5SUlIghMBnn32md86JEyfqnPPzzz/H3bt3Ub16ddSpUwcffvgh/vzzz1Kfm6eNecSIEahevToiIyNRpUoVDB48+Jn6S5T2HBQnrt7e3nr7tVqt9McDAPbt24fw8HCp75GHh4fUD6y43NO8fpcuXSr1vfLwuT7++GM4OjqiUaNGCAwMRFRUlNQk/DjJyck4deqU3nNdvXp1AE9+fwAP3iNl7ZtVoUIFhIeHo0ePHpg3bx46deqEdu3a4fr161IZOzu7Evuh3L9/XzpemuK+brVr1y61zOP+D9asWVPvC42VlVWJzb+A/v/Nkty8eRN5eXll+mwByvYeCgsLwxtvvIHY2FiUL18eXbp0weLFi8vUf+dFxD4qCvH3338jMzNTp632UXZ2dtizZw927tyJjRs3Ij4+Hr/88gteffVVJCQkwNLS8onXedp+JWXxuNqQssRkbKVdU/x//45npdVq4enpiZ9//rnE4x4eHka5TrFnec20Wi3atWuHjz76qMTjxX+Qij3puStOij/44IMSvw0CkN7DrVq1QmpqKtatW4eEhAR89913mDVrFubPn48hQ4Y8c8yenp44duwYtmzZgs2bN2Pz5s1YvHgx+vfvj6VLl5Z6/icp7Tl40nOTmpqKtm3bombNmpg5cya8vb1hY2ODTZs2YdasWQbXgpZFUFAQzp49iw0bNiA+Ph6//vorvv32W0yYMAGxsbGlPk6r1aJOnTqYOXNmiccfTc6M/X+rR48e+OSTT7Bu3ToMGzYMwIPa0mvXrumVLd5XqVIlg65lKLVarddfppixP0/L+h4qnozz4MGD+P3337FlyxYMHjwY//73v3Hw4EE4OjoaNS65MVFRiB9//BEASv3wL2ZhYYG2bduibdu2mDlzJqZOnYpPPvkEO3fuRHh4uNFnOXy02l4IgZSUFJ35Xtzc3HR67Re7dOmSTnONv78/Dh06hMLCwsd2Nn2Yh4cH7O3tcfbsWb1jf/31FywsLPQ+TA3h6+sLQP9+Aehd29/fH9u2bUPz5s2f+oOq+Dql3U/58uWfOPz44Vgffn5v3ryp983W398f2dnZUm3Esyq+nrW1dZnO6e7ujkGDBmHQoEHIzs5Gq1atMGnSpMcmKk8Ts42NDTp37ozOnTtDq9VixIgRWLBgAT777DMEBAQ811k/f//9d+Tn52P9+vU6NQ+PNp88zevn6+tb6nvl4XMBgIODA3r37o3evXujoKAA3bt3x5QpUxATE1PqcFl/f38cP34cbdu2lWWG1OJmnIdrpV555RX873//g1ar1UkQDh06BHt7e73k+mHFzVInT54s9Uvfw/8Hi5u4ip09e1bnOTUGDw8P2NnZlemzpazvoWJNmjRBkyZNMGXKFCxbtgx9+/bFihUrHvv/60XEph8F2LFjByZPnoyqVatKQ99Kcvv2bb19xZO6FVf5Ff+RKylxMMQPP/yg029m9erVuHbtGiIjI6V9/v7+OHjwoM6IhQ0bNug1ybzxxhu4desW5s6dq3ed0r6RWVpaon379li3bp1Oc9ONGzewbNkytGjRAs7OzobensTLywuvvPIKli5dqvOhuXXrVqkfQLFevXpBo9Fg8uTJeucpKip67HP/8HUeLnfy5EkkJCSgQ4cOT4w1PDwc1tbW+Oabb3Set0dH8BTHeuDAgRJ/F+nu3bs6/SvKwtPTE61bt8aCBQtK/Nb78DDwjIwMnWOOjo4ICAh4YvV0WWN+9PwWFhZSAm2q/w+PU1zb8PBrkpmZicWLF+uUe5rXr0OHDvjjjz9w4MABaV9OTg4WLlwIPz8/qX/Eo8+FjY0NatWqBSHEY/sw9OrVC1euXMGiRYv0juXl5ek1Qxrq1q1bJf4f/+677wBAZ26pHj164MaNGzojWG7duoVVq1ahc+fOJfZfKda+fXs4OTkhLi5OaioqVnz90NBQeHp6Yv78+Trvxc2bN+PMmTPo2LGjYTdZCktLS0REROC3337D5cuXpf1nzpzRe4+X9T10584dvefz0b8F5oQ1Ks/Z5s2b8ddff6GoqAg3btzAjh07sHXrVvj6+mL9+vWPnSjo888/x549e9CxY0f4+voiPT0d3377LapUqSJ1tvP394erqyvmz58PJycnODg4oHHjxmVqSy2Ju7s7WrRogUGDBuHGjRuYPXs2AgICdIZQDxkyBKtXr8Zrr72GXr16ITU1FT/99JPecOP+/fvjhx9+wNixY/HHH3+gZcuWyMnJwbZt2zBixAh06dKlxBi++OILaf6YESNGwMrKCgsWLEB+fj6+/PJLg+6rJHFxcejYsSNatGiBwYMH4/bt29LcFNnZ2VK5sLAwDBs2DHFxcTh27Bjat28Pa2trJCcnY9WqVZgzZ440UVVJvvrqK0RGRqJp06Z4++23peHJLi4uZfp9EA8PD3zwwQeIi4tDp06d0KFDBxw9ehSbN29G+fLldcp++OGHWL9+PTp16oSBAwciJCQEOTk5OHHiBFavXo2LFy/qPeZJ/vOf/6BFixaoU6cO3nnnHVSrVg03btzAgQMH8Pfff+P48eMAgFq1aqF169YICQmBu7s7kpKSsHr1aowcOfKx5y9rzEOGDMHt27fx6quvokqVKrh06RK++eYbvPLKK1IfjldeeQWWlpaYPn06MjMzoVarpTkqjK19+/ZSDc+wYcOQnZ2NRYsWwdPTUyepe5rXb9y4cVi+fDkiIyMRHR0Nd3d3LF26FBcuXMCvv/4q1Ti0b98eFStWRPPmzVGhQgWcOXMGc+fORceOHR/b7+1f//oXVq5ciXfffRc7d+5E8+bNodFo8Ndff2HlypXYsmWLThJhqJ9++gnz58+XOsXfu3cPW7ZswdatW9G5c2edmo0ePXqgSZMmGDRoEE6fPi3NTKvRaB7bjAUAzs7OmDVrFoYMGYKGDRtK8w8dP34cubm5WLp0KaytrTF9+nQMGjQIYWFh6NOnjzQ82c/PD2PGjHnm+31UbGws4uPj0bJlS4wYMQJFRUXSZ8vD/bbK+h5aunQpvv32W3Tr1g3+/v64d+8eFi1aBGdn5zJ92XnhPP+BRi+n4qGzxYuNjY2oWLGiaNeunZgzZ47OEOBijw5P3r59u+jSpYuoVKmSsLGxEZUqVRJ9+vTRG8a5bt06UatWLWFlZaUzVDksLKzUIYylDdVcvny5iImJEZ6ensLOzk507NhRZ4hdsX//+9+icuXKQq1Wi+bNm4ukpCS9cwrxYOjdJ598IqpWrSqsra1FxYoVRY8ePXSGHuORoX1CCHHkyBEREREhHB0dhb29vWjTpo3Yv39/ic/xo0P2iu+lpCGqj/r1119FUFCQUKvVolatWmLNmjViwIABesNQhRBi4cKFIiQkRNjZ2QknJydRp04d8dFHH4mrV68+8Trbtm0TzZs3F3Z2dsLZ2Vl07txZnD59WqdM8etf0nB0jUYjYmNjhZeXl7CzsxOtW7cWJ0+eLHGo+L1790RMTIwICAgQNjY2onz58qJZs2ZixowZ0jDc4uHJX331ld61Sno9UlNTRf/+/UXFihWFtbW1qFy5sujUqZNYvXq1VOaLL74QjRo1Eq6ursLOzk7UrFlTTJkyRWfob2nKEvPq1atF+/bthaenp7CxsRE+Pj5i2LBh4tq1azrnWrRokahWrZo0FLT4fVDae37VqlU6jy/tfVXS67N+/XpRt25dYWtrK/z8/MT06dOlYfkXLlyQyj3N65eamip69OghXF1dha2trWjUqJHYsGGDTpkFCxaIVq1aiXLlygm1Wi38/f3Fhx9+KDIzM5/4XBcUFIjp06eL4OBgoVarhZubmwgJCRGxsbE6jwcgoqKi9B5f2vQED0tMTBQ9e/YUPj4+Qq1WCwcHB9GgQQMxc+ZMUVhYqFf+9u3b4u233xblypUT9vb2Iiws7KmG4q5fv140a9ZM+v/VqFEjsXz5cp0yv/zyi6hfv75Qq9XC3d1d9O3bV/z99986ZQYMGCAcHBxKvEZpz0fxsUf/z+zevVuEhIQIGxsbUa1aNTF//ny9z/ji2J/0Hjpy5Ijo06eP9Hx6enqKTp06iaSkpDI/Ry8SlRBG6mFIREREZGTso0JERESKxUSFiIiIFIuJChERESkWExUiIiJSLCYqREREpFhMVIiIiEixXugJ37RaLa5evQonJydZpn8mIiKipyeEwL1791CpUqVSf0up2AudqFy9etUov/NCREREz19aWlqpv05d7IVOVIqnhk5LSzPK770QlVlODlD8K65XrwJP+CFBIiL6R1ZWFry9vR/7Ew/FXuhEpbi5x9nZmYkKPV8P/9y9szMTFSIiA5Sl2wY70xIREZFiMVEhIiIixXqhm36IZGNlBQwY8M86ERGZBD9hiQyhVgNLlsgdBdELTaPRoLCwUO4wyASsra1h+XBfvmfARIWIiJ4rIQSuX7+Ou3fvyh0KmZCrqysqVqz4zPOcMVEhMoQQQG7ug3V7e4ATDhKVWXGS4unpCXt7e07YaWaEEMjNzUV6ejoAwMvL65nOx0SFyBC5uYCj44P17GwOTyYqI41GIyUp5cqVkzscMhE7OzsAQHp6Ojw9PZ+pGYijfoiI6Lkp7pNib28vcyRkasWv8bP2Q2KiQkREzx2be8yfsV5jJipERESkWLImKn5+flCpVHpLVFSUnGERERGZ3JIlS+Dq6ip3GAbbtWsXVCqVyUdvyZqoJCYm4tq1a9KydetWAEDPnj3lDIuIiEhPWloaBg8ejEqVKsHGxga+vr4YNWoUMjIy5A7NrMmaqHh4eKBixYrSsmHDBvj7+yMsLEzOsIiIiHScP38eoaGhSE5OxvLly5GSkoL58+dj+/btaNq0KW7fvl3qYwsKCp5jpE+mtHieRDF9VAoKCvDTTz9h8ODB7GRFymdpCfTo8WAx0uyLRKRcUVFRsLGxQUJCAsLCwuDj44PIyEhs27YNV65cwSeffCKV9fPzw+TJk9G/f384Oztj6NChAB409fj4+MDe3h7dunUrsSZm3bp1aNCgAWxtbVGtWjXExsaiqKhIOq5SqfDdd9+hW7dusLe3R2BgINavX//Y2EuLZ+/evWjZsiXs7Ozg7e2N6Oho5OTkSI/78ccfERoaCicnJ1SsWBFvvfWWNDfKcyUU4pdffhGWlpbiypUrpZa5f/++yMzMlJa0tDQBQGRmZj7HSImIyFB5eXni9OnTIi8vT/dAdnbpy9OUzc0tW9mnkJGRIVQqlZg6dWqJx9955x3h5uYmtFqtEEIIX19f4ezsLGbMmCFSUlJESkqKOHjwoLCwsBDTp08XZ8+eFXPmzBGurq7CxcVFOs+ePXuEs7OzWLJkiUhNTRUJCQnCz89PTJo0SSoDQFSpUkUsW7ZMJCcni+joaOHo6CgyMjJKjb+keFJSUoSDg4OYNWuWOHfunNi3b5+oX7++GDhwoPS4//73v2LTpk0iNTVVHDhwQDRt2lRERkZKx3fu3CkAiDt37pR43VJfayFEZmZmmf9+q/7/xmUXEREBGxsb/P7776WWmTRpEmJjY/X2Z2ZmwtnZ2ZThEb1Q/MZtlDuEJ7o4raPcIZAM7t+/jwsXLqBq1aqwtbX958DjatI7dAA2PvSednD4Z2boR4WFAbt2/bPt4QHcuqVf7in+9B06dAhNmjTB2rVr0bVrV73js2bNwtixY3Hjxg14enrCz88P9evXx9q1a6Uyb731FjIzM7Hxoft48803ER8fL3VGDQ8PR9u2bRETEyOV+emnn/DRRx/h6tWrAB7UqHz66aeYPHkyACAnJweOjo7YvHkzXnvttRLjLymeIUOGwNLSEgsWLJD27d27F2FhYcjJydF9bf5fUlISGjZsiHv37sHR0RG7du1CmzZtcOfOnRI7BZf6WgPIysqCi4tLmf5+K6Lp59KlS9i2bRuGDBny2HIxMTHIzMyUlrS0tOcUIRERveye5nt9aGiozvaZM2fQuHFjnX1NmzbV2T5+/Dg+//xzODo6Sss777yDa9euIfehxKxu3brSuoODA5ydnZ/YJPNoPMePH8eSJUt0rhUREQGtVosLFy4AAA4fPozOnTvDx8cHTk5OUv/Ry5cvl/FZMA5FTKG/ePFieHp6omPHx3/DUqvVUKvVzykqosfIyeEU+kTGlJ1d+rFH+4E97o+yxSPfvy9eNDikYgEBAVCpVDhz5gy6deumd/zMmTNwc3ODh4eHtM/BgM+E7OxsxMbGonv37nrHHq6RsLa21jmmUqmg1Wofe+5H48nOzsawYcMQHR2tV9bHxwc5OTmIiIhAREQEfv75Z3h4eODy5cuIiIh47p1xZU9UtFotFi9ejAEDBsDKSvZwiIhIDk/zh91UZUtRrlw5tGvXDt9++y3GjBkj/Y4N8OAHFn/++Wf079//sQNBgoKCcOjQIZ19Bw8e1Nlu0KABzp49i4CAgGeO+UkaNGiA06dPl3qtEydOICMjA9OmTYO3tzeAB00/cpC96Wfbtm24fPkyBg8eLHcoREREJZo7dy7y8/MRERGBPXv2IC0tDfHx8WjXrh0qV66MKVOmPPbx0dHRiI+Px4wZM5CcnIy5c+ciPj5ep8yECRPwww8/IDY2FqdOncKZM2ewYsUKfPrpp0a/n48//hj79+/HyJEjcezYMSQnJ2PdunUYOXIkgAe1KjY2Nvjmm29w/vx5rF+/XuoX87zJnqi0b98eQghUr15d7lCIiIhKFBgYiKSkJFSrVg29evWCv78/hg4dijZt2uDAgQNwd3d/7OObNGmCRYsWYc6cOahXrx4SEhL0EpCIiAhs2LABCQkJaNiwIZo0aYJZs2bB19fX6PdTt25d7N69G+fOnUPLli1Rv359TJgwAZUqVQLwYJ6zJUuWYNWqVahVqxamTZuGGTNmGD2OslDMqB9DPE2vYSKjUngfFY76IaV63EgQMi9mNeqHiIiIqCRMVIiIiEixOMyGyBCWlg8moSpeJyIik2CiQmQIW1vdmTKJiMgk2PRDRETP3Qs8joPKyFivMRMVIiJ6bopnVc0t7bd6yGwUv8aPzqT7tNj0Q2SInBzA0/PBenq64oYnEymVpaUlXF1dpd+msbe3f+yMrvTiEUIgNzcX6enpcHV1heUz9uNjokJkKH4jJDJIxYoVAeCJP6RHLzZXV1fptX4WTFSIiOi5UqlU8PLygqenJwoLC+UOh0zA2tr6mWtSijFRISIiWVhaWhrtjxmZL3amJSIiIsViokJERESKxUSFiIiIFIt9VIgMYWEBhIX9s05ERCbBRIXIEHZ2wK5dckdBRGT2+FWQiIiIFIuJChERESkWExUiQ+TkAB4eD5acHLmjISIyW+yjQmSoW7fkjoCIyOyxRoWIiIgUi4kKERERKRYTFSIiIlIsJipERESkWExUiIiISLE46ofIEBYWQGjoP+tERGQSTFSIDGFnByQmyh0FEZHZ41dBIiIiUiwmKkRERKRYTFSIDJGbC/j5PVhyc+WOhojIbLGPCpEhhAAuXfpnnYiITII1KkRERKRYTFSIiIhIsZioEBERkWIxUSEiIiLFYqJCREREisVRP0SGUKmAWrX+WSciIpNgokJkCHt74NQpuaMgIjJ7bPohIiIixWKiQkRERIole6Jy5coV9OvXD+XKlYOdnR3q1KmDpKQkucMierzcXCA4+MHCKfSJiExG1j4qd+7cQfPmzdGmTRts3rwZHh4eSE5Ohpubm5xhET2ZEMDp0/+sExGRSciaqEyfPh3e3t5YvHixtK9q1aoyRkRERERKImvTz/r16xEaGoqePXvC09MT9evXx6JFi0otn5+fj6ysLJ2FiIiIzJesicr58+cxb948BAYGYsuWLRg+fDiio6OxdOnSEsvHxcXBxcVFWry9vZ9zxERERPQ8qYSQr4HdxsYGoaGh2L9/v7QvOjoaiYmJOHDggF75/Px85OfnS9tZWVnw9vZGZmYmnJ2dn0vMRACAnBzA0fHBenY24OAgbzyP8Bu3Ue4QnujitI5yh0BEMsnKyoKLi0uZ/n7LWqPi5eWFWsWze/6/oKAgXL58ucTyarUazs7OOgsRERGZL1k70zZv3hxnz57V2Xfu3Dn4+vrKFBFRGalUQPH7lFPoExGZjKyJypgxY9CsWTNMnToVvXr1wh9//IGFCxdi4cKFcoZF9GT29sDFi3JHQURk9mRt+mnYsCHWrl2L5cuXo3bt2pg8eTJmz56Nvn37yhkWERERKYTsP0rYqVMndOrUSe4wiIiISIFkn0Kf6IWUlwc0bPhgycuTOxoiIrMle40K0QtJqwWKf5NKq5U3FiIiM8YaFSIiIlIsJipERESkWExUiIiISLGYqBAREZFiMVEhIiIixeKoHyJDlS8vdwRERGaPiQqRIRwcgJs35Y6CiMjssemHiIiIFIuJChERESkWExUiQ+TlAa1bP1g4hT4RkcmwjwqRIbRaYPfuf9aJiMgkWKNCREREisVEhYiIiBSLiQoREREpFhMVIiIiUiwmKkRERKRYHPVDZCh7e7kjICIye0xUiAzh4ADk5MgdBRGR2WPTDxERESkWExUiIiJSLCYqRIa4fx/o2PHBcv++3NEQEZkt9lEhMoRGA2za9M86ERGZBGtUiIiISLGYqBAREZFiMVEhIiIixWKiQkRERIrFRIWIiIgUi4kKERERKRaHJxMZwsEBEELuKIiIzB5rVIiIiEixmKgQERGRYjFRITLE/ftAz54PFk6hT0RkMkxUiAyh0QCrVz9YOIU+EZHJMFEhIiIixWKiQkRERIrFRIWIiIgUi4kKERERKZasicqkSZOgUql0lpo1a8oZEhERESmI7DPTBgcHY9u2bdK2lZXsIREREZFCyJ4VWFlZoWLFinKHQfR07O2B7Ox/1omIyCRk76OSnJyMSpUqoVq1aujbty8uX75catn8/HxkZWXpLESyUKke/N6Pg8ODdSIiMglZa1QaN26MJUuWoEaNGrh27RpiY2PRsmVLnDx5Ek5OTnrl4+LiEBsbK0OkZC78xm2UO4Qnujito9whkALxvUsvK1lrVCIjI9GzZ0/UrVsXERER2LRpE+7evYuVK1eWWD4mJgaZmZnSkpaW9pwjJnrApqgQMzbOwoyNs2BTVCh3OEREZkv2PioPc3V1RfXq1ZGSklLicbVaDbVa/ZyjItJnqdWgx8ntAIDP2g0HYC1vQEREZkr2PioPy87ORmpqKry8vOQOhYiIiBRA1kTlgw8+wO7du3Hx4kXs378f3bp1g6WlJfr06SNnWERERKQQsjb9/P333+jTpw8yMjLg4eGBFi1a4ODBg/Dw8JAzLCIiIlIIWROVFStWyHl5IiIiUjhF9VEhIiIiehgTFSIiIlIsRQ1PJnpR5Fmr0eC9n6V1IiIyDSYqRIZQqXDb3kXuKIiIzB6bfoiIiEixWKNCZACbokJ8uuM7AMAXrw5BgRVnpiUiMgXWqBAZwFKrQf+jG9H/6EZYajVyh0NEZLaYqBAREZFiMVEhIiIixWKiQkRERIrFRIWIiIgUi4kKERERKRYTFSIiIlIszqNCZID71jZo8e5/pXUiIjINJipEBhAqC/ztUkHuMIiIzB6bfoiIiEixWKNCZABrTSE+2PMjAGBGq3+h0JJT6BMRmQJrVIgMYKXRYNgfazDsjzWw0nAKfSIiU2GiQkRERIrFRIWIiIgUi4kKERERKRYTFSIiIlIsJipERESkWExUiIiISLE4jwqRAe5b26Dd4P9I60REZBpMVIgMIFQWSPbwlTsMIiKzx6YfIiIiUiyDa1RycnKwe/duXL58GQUFBTrHoqOjnzkwIiWz1hQi6sBKAMB/mvbiFPpERCZiUKJy9OhRdOjQAbm5ucjJyYG7uztu3boFe3t7eHp6MlEhs2el0WD0vuUAgAWN3mCiQkRkIgY1/YwZMwadO3fGnTt3YGdnh4MHD+LSpUsICQnBjBkzjB0jERERvaQMSlSOHTuG999/HxYWFrC0tER+fj68vb3x5ZdfYvz48caOkYiIiF5SBiUq1tbWsLB48FBPT09cvnwZAODi4oK0tDTjRUdEREQvNYP6qNSvXx+JiYkIDAxEWFgYJkyYgFu3buHHH39E7dq1jR0jERERvaQMqlGZOnUqvLy8AABTpkyBm5sbhg8fjps3b2LhwoVGDZCIiIheXgbVqISGhkrrnp6eiI+PN1pARERERMU4My2RAfKtrPF6/5nSOhERmUaZE5UGDRpg+/btcHNzQ/369aFSqUote+TIEaMER6RUWgtL/OlVXe4wiIjMXpkTlS5dukCtVgMAunbtaqp4iIiIiCRlTlQmTpxY4jrRy8haU4hBSesBAItDX+fMtEREJmLQqJ/ExEQcOnRIb/+hQ4eQlJRkUCDTpk2DSqXC6NGjDXo80fNkpdFg/K7FGL9rMaw0GrnDISIyWwYlKlFRUSVO7HblyhVERUU99fkSExOxYMEC1K1b15BwiIiIyEwZlKicPn0aDRo00Ntfv359nD59+qnOlZ2djb59+2LRokVwc3MzJBwiIiIyUwYlKmq1Gjdu3NDbf+3aNVhZPd2I56ioKHTs2BHh4eGGhEJERERmzKB5VNq3b4+YmBisW7cOLi4uAIC7d+9i/PjxaNeuXZnPs2LFChw5cgSJiYllKp+fn4/8/HxpOysr6+kCJyIioheKQYnKjBkz0KpVK/j6+qJ+/foAHvyicoUKFfDjjz+W6RxpaWkYNWoUtm7dCltb2zI9Ji4uDrGxsYaEbBC/cRuf27UMdXFaR7lDIJIV/58SmTeDmn4qV66MP//8E19++SVq1aqFkJAQzJkzBydOnIC3t3eZznH48GGkp6ejQYMGsLKygpWVFXbv3o2vv/4aVlZW0JQwkiImJgaZmZnSwl9qJiIiMm8GT6Hv4OCAoUOHGnzhtm3b4sSJEzr7Bg0ahJo1a+Ljjz+GpaWl3mPUarU06RyRnPKtrPFmn6nSOhERmYbBiUpycjJ27tyJ9PR0aLVanWMTJkx44uOdnJxQu3ZtnX0ODg4oV66c3n4ipdFaWOKgD4fTExGZmkGJyqJFizB8+HCUL18eFStW1PndH5VKVaZEhYiIiOhJDEpUvvjiC0yZMgUff/yxUYPZtWuXUc9HZCpWmiL0OR4PAFhe7zUUWfKHyImITMGgT9c7d+6gZ8+exo6F6IVhrSnC5K3zAQCra4czUSEiMhGDRv307NkTCQkJxo6FiIiISIdBXwMDAgLw2Wef4eDBg6hTpw6srXVHPURHRxslOCIiInq5GZSoLFy4EI6Ojti9ezd2796tc0ylUjFRISIiIqMwKFG5cOGCseMgIiIi0mNQH5ViBQUFOHv2LIqKiowVDxEREZHEoEQlNzcXb7/9Nuzt7REcHIzLly8DAN577z1MmzbNqAESERHRy8ugRCUmJgbHjx/Hrl27dH5QMDw8HL/88ovRgiNSqgIrawzqMRGDekxEAafQJyIyGYP6qPz222/45Zdf0KRJE51ZaYODg5Gammq04IiUSmNhiZ3+DeUOg4jI7BlUo3Lz5k14enrq7c/JydFJXIiIiIiehUGJSmhoKDZu3ChtFycn3333HZo2bWqcyIgUzEpThB4ntqHHiW2w0rAzORGRqRjU9DN16lRERkbi9OnTKCoqwpw5c3D69Gns379fb14VInNkrSnCjE2zAQAba7TgFPpERCZiUI1KixYtcOzYMRQVFaFOnTpISEiAp6cnDhw4gJCQEGPHSERERC8pg78G+vv7Y9GiRcaMhYiIiEiHQYlK8bwppfHx8TEoGCIiIqKHGZSo+Pn5PXZ0j0ajMTggIiIiomIGJSpHjx7V2S4sLMTRo0cxc+ZMTJkyxSiBERERERmUqNSrV09vX2hoKCpVqoSvvvoK3bt3f+bAiIiIiIw6prJGjRpITEw05imJFKnAyhojuoyT1omIyDQMSlSysrJ0toUQuHbtGiZNmoTAwECjBEakZBoLS2yq2ULuMIiIzJ5BiYqrq6teZ1ohBLy9vbFixQqjBEZERERkUKKyY8cOnUTFwsICHh4eCAgIgJUVZ+gk82ep1SDi3AEAwJbqTaGxsJQ5IiIi82RQVtG6dWsjh0H0YrEpKsS366YBAILGrEaeDRMVIiJTMGgK/bi4OHz//fd6+7///ntMnz79mYMiIiIiAgxMVBYsWICaNWvq7Q8ODsb8+fOfOSgiIiIiwMBE5fr16/Dy8tLb7+HhgWvXrj1zUERERESAgYmKt7c39u3bp7d/3759qFSp0jMHRURERAQY2Jn2nXfewejRo1FYWIhXX30VALB9+3Z89NFHeP/9940aIBEREb28DEpUPvzwQ2RkZGDEiBEoKCgAANja2uLjjz9GTEyMUQMkIiKil5dBiYpKpcL06dPx2Wef4cyZM7Czs0NgYCDUarWx4yNSpEJLK3zQYbS0TkREpvFMn7DXr1/H7du30apVK6jVaggh9GasJTJHRZZWWF0nXO4wiIjMnkGdaTMyMtC2bVtUr14dHTp0kEb6vP322+yjQkREREZjUKIyZswYWFtb4/Lly7C3t5f29+7dG/Hx8UYLjkipLLUatElNRJvURFhqNXKHQ0Rktgxq+klISMCWLVtQpUoVnf2BgYG4dOmSUQIjUjKbokIsXh0LgFPoExGZkkE1Kjk5OTo1KcVu377NDrVERERkNAYlKi1btsQPP/wgbatUKmi1Wnz55Zdo06aN0YIjIiKil5tBTT9ffvkl2rZti6SkJBQUFOCjjz7CqVOncPv27RJnrCUiIiIyhEE1KrVr18a5c+fQokULdOnSBTk5OejevTuOHj0Kf39/Y8dIREREL6mnrlEpLCzEa6+9hvnz5+OTTz4xRUxEREREAAyoUbG2tsaff/5piliIiIiIdBjU9NOvXz/897//feaLz5s3D3Xr1oWzszOcnZ3RtGlTbN68+ZnPS2RqhZZW+Kzdu/is3bucQp+IyIQM+oQtKirC999/j23btiEkJAQODg46x2fOnFmm81SpUgXTpk1DYGAghBBYunQpunTpgqNHjyI4ONiQ0IieiyJLK/zYoJPcYRARmb2nSlTOnz8PPz8/nDx5Eg0aNAAAnDt3TqfM0/zWT+fOnXW2p0yZgnnz5uHgwYNMVIiIiOjpEpXAwEBcu3YNO3fuBPBgyvyvv/4aFSpUeOZANBoNVq1ahZycHDRt2rTEMvn5+cjPz5e2s7Kynvm6RIaw0GrQ6O9TAIA/qgRDa8GZaYmITOGpEhUhhM725s2bkZOT80wBnDhxAk2bNsX9+/fh6OiItWvXolatWiWWjYuLQ2xs7DNdj8gY1EWFWLF8PABOoU9EZEoGdaYt9mjiYogaNWrg2LFjOHToEIYPH44BAwbg9OnTJZaNiYlBZmamtKSlpT3z9YmIiEi5nqpGRaVS6fVBeZo+KSWxsbFBQEAAACAkJASJiYmYM2cOFixYoFdWrVbzt4SIiIheIk/d9DNw4EApWbh//z7effddvVE/a9asMTggrVar0w+FiIiIXl5PlagMGDBAZ7tfv37PdPGYmBhERkbCx8cH9+7dw7Jly7Br1y5s2bLlmc5LRERE5uGpEpXFixcb9eLp6eno378/rl27BhcXF9StWxdbtmxBu3btjHodIiIiejHJOqWmMWa3JSIiIvPFub+JDFBkaYmprQdJ60REZBpMVIgMUGhpjYWN35A7DCIis/dM86gQERERmRJrVIgMYKHVoPaNVADAyQr+nEKfiMhEmKgQGUBdVIj1P4wFwCn0iYhMiU0/REREpFhMVIiIiEixmKgQERGRYjFRISIiIsViokJERESKxUSFiIiIFIvDk4kMUGRpidnN+0jrRERkGkxUiAxQaGmN2S36yh0GEZHZY9MPERERKRZrVIgMoBJaBNxKAwCklPeGUDHnJyIyBSYqRAawLSzA1u+jABRPoW8rc0REROaJXwOJiIhIsZioEBERkWIxUSEiIiLFYqJCREREisVEhYiIiBSLiQoREREpFocnExmgyNISCxp1l9aJiMg0mKgQGaDQ0hpxbQbLHQYRkdlj0w8REREpFmtUiAygElpUzroJALji7MEp9ImITISJCpEBbAsLsHf+2wA4hT4RkSnxayAREREpFhMVIiIiUiwmKkRERKRYTFSIiIhIsZioEBERkWIxUSEiIiLF4vBkIgNoLCzxQ/2O0joREZkGExUiAxRYWWNC++Fyh0FEZPbY9ENERESKxRoVIkMIAfe8LADAbTtnQKWSOSAiIvPERIXIAHaF+TjyTV8AnEKfiMiU2PRDREREiiVrohIXF4eGDRvCyckJnp6e6Nq1K86ePStnSERERKQgsiYqu3fvRlRUFA4ePIitW7eisLAQ7du3R05OjpxhERERkULI2kclPj5eZ3vJkiXw9PTE4cOH0apVK5miIiIiIqVQVGfazMxMAIC7u3uJx/Pz85Gfny9tZ2VlPZe4iIiISB6KSVS0Wi1Gjx6N5s2bo3bt2iWWiYuLQ2xs7HOOjIiInie/cRvlDuGJLk7raJTzvEz3aijFJCpRUVE4efIk9u7dW2qZmJgYjB07VtrOysqCt7f38wiPSIfGwhKra7eV1omIyDQUkaiMHDkSGzZswJ49e1ClSpVSy6nVaqjV6ucYGVHJCqys8UHHMXKHQURk9mRNVIQQeO+997B27Vrs2rULVatWlTMcIiIiUhhZE5WoqCgsW7YM69atg5OTE65fvw4AcHFxgZ2dnZyhET2eELArfNCxO89azSn0iYhMRNZ5VObNm4fMzEy0bt0aXl5e0vLLL7/IGRbRE9kV5uPMrB44M6uHlLAQEZHxyd70Q0RERFQa/tYPERERKRYTFSIiIlIsJipERESkWExUiIiISLGYqBAREZFiKWJmWqIXjdbCAhtrNJfWiYjINJioEBkg38oGUV1j5A6DiMjs8asgERERKRYTFSIiIlIsJipEBrAruI+L0zvh4vROsCu4L3c4RERmi4kKERERKRYTFSIiIlIsJipERESkWExUiIiISLGYqBAREZFiMVEhIiIixeLMtEQG0FpYYEe1UGmdiIhMg4kKkQHyrWwwuOckucMgIjJ7/CpIREREisVEhYiIiBSLiQqRAewK7uP0zDdweuYbnEKfiMiE2EeFyED2hflyh0BEZPZYo0JERESKxUSFiIiIFIuJChERESkWExUiIiJSLCYqREREpFgc9UNkAK1KhYPetaV1IiIyDSYqRAbIt1bjzbemyR0GEZHZY9MPERERKRYTFSIiIlIsJipEBrAruI/DX7+Fw1+/xSn0iYhMiH1UiAxULi9L7hCIiMwea1SIiIhIsZioEBERkWIxUSEiIiLFYqJCREREisVEhYiIiBSLo36IDKBVqXC8YqC0TkREpiFrjcqePXvQuXNnVKpUCSqVCr/99puc4RCVWb61Gl0GzEKXAbOQb62WOxwiIrMla6KSk5ODevXq4T//+Y+cYRAREZFCydr0ExkZicjISDlDICIiIgV7ofqo5OfnIz8/X9rOyuLMoCQP28L72PbdCABA+JBvcd/aVuaIiIjM0wuVqMTFxSE2NlbuMMyO37iNcofwRBendZQ7BB0qAVTJSpfWiYjINF6o4ckxMTHIzMyUlrS0NLlDIiIiIhN6oWpU1Go11GqOsCAiInpZvFA1KkRERPRykbVGJTs7GykpKdL2hQsXcOzYMbi7u8PHx0fGyIiIiEgJZE1UkpKS0KZNG2l77NixAIABAwZgyZIlMkVFRERESiFrotK6dWsIwSET9OIRKuBcOR9pnYiITOOF6kxLpBT3rW3Rfsi3codBRGT22JmWiIiIFIuJChERESkWExUiA9gW3kfCdyOQ8N0I2BbelzscIiKzxT4qRAZQCaB6xmVpnYiITIM1KkRERKRYTFSIiIhIsZioEBERkWIxUSEiIiLFYqJCREREisVRP0QGECrgb2dPaZ2IiEyDiQqRAe5b26LF8O/lDoOIyOyx6YeIiIgUi4kKERERKRYTFSIDqAvzsW7pGKxbOgbqwny5wyEiMlvso0JkAAshUO96srRORESmwRoVIiIiUiwmKkRERKRYTFSIiIhIsZioEBERkWIxUSEiIiLF4qgfIgNl2DnLHQIRkdljokJkgDwbW4REL5M7DCIis8emHyIiIlIsJipERESkWExUiAygLszHimXjsGLZOE6hT0RkQuyjQmQACyHQJO2ktE5ERKbBGhUiIiJSLCYqREREpFhMVIiIiEixmKgQERGRYjFRISIiIsXiqB8iA+Vaq+UOgYjI7DFRITJAno0tao39Ve4wiIjMHpt+iIiISLGYqBAREZFisemHyADqogLMWzsVADC823jkW9nIHBERkXliokJkAAutFq+eT5LWiYjINNj0Q0RERIqliETlP//5D/z8/GBra4vGjRvjjz/+kDskIiIiUgDZE5VffvkFY8eOxcSJE3HkyBHUq1cPERERSE9Plzs0IiIikpnsicrMmTPxzjvvYNCgQahVqxbmz58Pe3t7fP/993KHRkRERDKTNVEpKCjA4cOHER4eLu2zsLBAeHg4Dhw4IGNkREREpASyjvq5desWNBoNKlSooLO/QoUK+Ouvv/TK5+fnIz8/X9rOzMwEAGRlZZkkPm1+rknOa0zGuPeX5T4B492rpuA+iiPS5OdCK4w38oev6dN5We71ZblPgPeqNKb4G1t8TiHEkwsLGV25ckUAEPv379fZ/+GHH4pGjRrplZ84caIAwIULFy5cuHAxgyUtLe2JuYKsNSrly5eHpaUlbty4obP/xo0bqFixol75mJgYjB07VtrWarW4ffs2ypUrB5VKZfJ4n0VWVha8vb2RlpYGZ2dnucMxqZflXl+W+wRennt9We4T4L2aoxfpPoUQuHfvHipVqvTEsrImKjY2NggJCcH27dvRtWtXAA+Sj+3bt2PkyJF65dVqNdRq3V+sdXV1fQ6RGo+zs7Pi30DG8rLc68tyn8DLc68vy30CvFdz9KLcp4uLS5nKyT4z7dixYzFgwACEhoaiUaNGmD17NnJycjBo0CC5QyMiIiKZyZ6o9O7dGzdv3sSECRNw/fp1vPLKK4iPj9frYEtEREQvH9kTFQAYOXJkiU095kStVmPixIl6TVfm6GW515flPoGX515flvsEeK/myFzvUyVEWcYGERERET1/ss9MS0RERFQaJipERESkWExUiIiISLGYqBAREZFiMVExsT179qBz586oVKkSVCoVfvvtN7lDMom4uDg0bNgQTk5O8PT0RNeuXXH27Fm5wzKJefPmoW7dutKkSk2bNsXmzZvlDsvkpk2bBpVKhdGjR8sditFNmjQJKpVKZ6lZs6bcYZnMlStX0K9fP5QrVw52dnaoU6cOkpKS5A7LqPz8/PReU5VKhaioKLlDMzqNRoPPPvsMVatWhZ2dHfz9/TF58uSy/Y7OC0ARw5PNWU5ODurVq4fBgweje/fucodjMrt370ZUVBQaNmyIoqIijB8/Hu3bt8fp06fh4OAgd3hGVaVKFUybNg2BgYEQQmDp0qXo0qULjh49iuDgYLnDM4nExEQsWLAAdevWlTsUkwkODsa2bdukbSsr8/x4vHPnDpo3b442bdpg8+bN8PDwQHJyMtzc3OQOzagSExOh0Wik7ZMnT6Jdu3bo2bOnjFGZxvTp0zFv3jwsXboUwcHBSEpKwqBBg+Di4oLo6Gi5w3tm5vk/UUEiIyMRGRkpdxgmFx8fr7O9ZMkSeHp64vDhw2jVqpVMUZlG586ddbanTJmCefPm4eDBg2aZqGRnZ6Nv375YtGgRvvjiC7nDMRkrK6sSf2PM3EyfPh3e3t5YvHixtK9q1aoyRmQaHh4eOtvTpk2Dv78/wsLCZIrIdPbv348uXbqgY8eOAB7UJi1fvhx//PGHzJEZB5t+yCQyMzMBAO7u7jJHYloajQYrVqxATk4OmjZtKnc4JhEVFYWOHTsiPDxc7lBMKjk5GZUqVUK1atXQt29fXL58We6QTGL9+vUIDQ1Fz5494enpifr162PRokVyh2VSBQUF+OmnnzB48GDF/4CtIZo1a4bt27fj3LlzAIDjx49j7969ZvMlmTUqZHRarRajR49G8+bNUbt2bbnDMYkTJ06gadOmuH//PhwdHbF27VrUqlVL7rCMbsWKFThy5AgSExPlDsWkGjdujCVLlqBGjRq4du0aYmNj0bJlS5w8eRJOTk5yh2dU58+fx7x58zB27FiMHz8eiYmJiI6Oho2NDQYMGCB3eCbx22+/4e7duxg4cKDcoZjEuHHjkJWVhZo1a8LS0hIajQZTpkxB37595Q7NKJiokNFFRUXh5MmT2Lt3r9yhmEyNGjVw7NgxZGZmYvXq1RgwYAB2795tVslKWloaRo0aha1bt8LW1lbucEzq4W+edevWRePGjeHr64uVK1fi7bffljEy49NqtQgNDcXUqVMBAPXr18fJkycxf/58s01U/vvf/yIyMhKVKlWSOxSTWLlyJX7++WcsW7YMwcHBOHbsGEaPHo1KlSqZxWvKRIWMauTIkdiwYQP27NmDKlWqyB2OydjY2CAgIAAAEBISgsTERMyZMwcLFiyQOTLjOXz4MNLT09GgQQNpn0ajwZ49ezB37lzk5+fD0tJSxghNx9XVFdWrV0dKSorcoRidl5eXXkIdFBSEX3/9VaaITOvSpUvYtm0b1qxZI3coJvPhhx9i3LhxePPNNwEAderUwaVLlxAXF8dEhaiYEALvvfce1q5di127dpll57zH0Wq1yM/PlzsMo2rbti1OnDihs2/QoEGoWbMmPv74Y7NNUoAHHYhTU1Pxr3/9S+5QjK558+Z6UwecO3cOvr6+MkVkWosXL4anp6fU0dQc5ebmwsJCt8uppaUltFqtTBEZFxMVE8vOztb5VnbhwgUcO3YM7u7u8PHxkTEy44qKisKyZcuwbt06ODk54fr16wAAFxcX2NnZyRydccXExCAyMhI+Pj64d+8eli1bhl27dmHLli1yh2ZUTk5Oen2MHBwcUK5cObPre/TBBx+gc+fO8PX1xdWrVzFx4kRYWlqiT58+codmdGPGjEGzZs0wdepU9OrVC3/88QcWLlyIhQsXyh2a0Wm1WixevBgDBgww2+HmwIORiFOmTIGPjw+Cg4Nx9OhRzJw5E4MHD5Y7NOMQZFI7d+4UAPSWAQMGyB2aUZV0jwDE4sWL5Q7N6AYPHix8fX2FjY2N8PDwEG3bthUJCQlyh/VchIWFiVGjRskdhtH17t1beHl5CRsbG1G5cmXRu3dvkZKSIndYJvP777+L2rVrC7VaLWrWrCkWLlwod0gmsWXLFgFAnD17Vu5QTCorK0uMGjVK+Pj4CFtbW1GtWjXxySefiPz8fLlDMwqVEGYydR0RERGZHc6jQkRERIrFRIWIiIgUi4kKERERKRYTFSIiIlIsJipERESkWExUiIiISLGYqBAREZFiMVEhIj1LliyBq6ur3GEYpHXr1hg9evRLd20ic8VEhcgMpKWlYfDgwahUqRJsbGzg6+uLUaNGISMjQ+7QTGbXrl1QqVS4e/euzv41a9Zg8uTJ8gRFREbHRIXoBXf+/HmEhoYiOTkZy5cvR0pKCubPn4/t27ejadOmuH37dqmPLSgoeI6RPh/u7u5wcnKSOwwiMhImKkQvuKioKNjY2CAhIQFhYWHw8fFBZGQktm3bhitXruCTTz6Ryvr5+WHy5Mno378/nJ2dMXToUAAPmnp8fHxgb2+Pbt26lVgTs27dOjRo0AC2traoVq0aYmNjUVRUJB1XqVT47rvv0K1bN9jb2yMwMBDr169/bOz5+fn44IMPULlyZTg4OKBx48bYtWuXdPzSpUvo3Lkz3Nzc4ODggODgYGzatAkXL15EmzZtAABubm5QqVQYOHAgAP3mFz8/P3zxxRfo378/HB0d4evri/Xr1+PmzZvo0qULHB0dUbduXSQlJUmPycjIQJ8+fVC5cmXY29ujTp06WL58uU7sOTk50jm9vLzw73//W+/+7ty5g/79+8PNzQ329vaIjIxEcnLyY58TInqE3D82RESGy8jIECqVSkydOrXE4++8845wc3MTWq1WCCGEr6+vcHZ2FjNmzBApKSkiJSVFHDx4UFhYWIjp06eLs2fPijlz5ghXV1fh4uIinWfPnj3C2dlZLFmyRKSmpoqEhATh5+cnJk2aJJUBIKpUqSKWLVsmkpOTRXR0tHB0dBQZGRmlxj9kyBDRrFkzsWfPHpGSkiK++uoroVarxblz54QQQnTs2FG0a9dO/PnnnyI1NVX8/vvvYvfu3aKoqEj8+uuv0g/OXbt2Tdy9e1cIof/Dib6+vsLd3V3Mnz9fnDt3TgwfPlw4OzuL1157TaxcuVKcPXtWdO3aVQQFBUnP099//y2++uorcfToUZGamiq+/vprYWlpKQ4dOiSdd/jw4cLHx0ds27ZN/Pnnn6JTp07CyclJ59qvv/66CAoKEnv27BHHjh0TERERIiAgQBQUFJTtBSYiwUSF6AV28OBBAUCsXbu2xOMzZ84UAMSNGzeEEA/+aHft2lWnTJ8+fUSHDh109vXu3VsnUWnbtq1eMvTjjz8KLy8vaRuA+PTTT6Xt7OxsAUBs3ry5xNguXbokLC0txZUrV3T2t23bVsTExAghhKhTp45OMvSw4l8mv3Pnjs7+khKVfv36SdvXrl0TAMRnn30m7Ttw4IAAIK5du1bitYR4kDS9//77Qggh7t27J2xsbMTKlSul4xkZGcLOzk669rlz5wQAsW/fPqnMrVu3hJ2dnc7jiOjxrOSqySEi4xFP8SPooaGhOttnzpxBt27ddPY1bdoU8fHx0vbx48exb98+TJkyRdqn0Whw//595Obmwt7eHgBQt25d6biDgwOcnZ2Rnp5eYhwnTpyARqNB9erVdfbn5+ejXLlyAIDo6GgMHz4cCQkJCA8PxxtvvKFzjbJ6+DEVKlQAANSpU0dvX3p6OipWrAiNRoOpU6di5cqVuHLlCgoKCpCfny/dZ2pqKgoKCtC4cWPpHO7u7qhRo4a0febMGVhZWemUKVeuHGrUqIEzZ8489T0QvayYqBC9wAICAqBSqUpMNoAHfyzd3Nzg4eEh7XNwcHjq62RnZyM2Nhbdu3fXO2ZrayutW1tb6xxTqVTQarWlntPS0hKHDx+GpaWlzjFHR0cAwJAhQxAREYGNGzciISEBcXFx+Pe//4333nvvqeJ/OC6VSlXqvuJYv/rqK8yZMwezZ89GnTp14ODggNGjR5tl52MipWNnWqIXWLly5dCuXTt8++23yMvL0zl2/fp1/Pzzz+jdu7f0h7gkQUFBOHTokM6+gwcP6mw3aNAAZ8+eRUBAgN5iYWHYx0j9+vWh0WiQnp6ud86KFStK5by9vfHuu+9izZo1eP/997Fo0SIAgI2NDYAHNTvGtm/fPnTp0gX9+vVDvXr1UK1aNZw7d0467u/vD2tra53n7c6dOzplgoKCUFRUpFMmIyMDZ8+eRa1atYweM5G5YqJC9IKbO3cu8vPzERERgT179iAtLQ3x8fFo164dKleurNNcU5Lo6GjEx8djxowZSE5Oxty5c3WafQBgwoQJ+OGHHxAbG4tTp07hzJkzWLFiBT799FOD465evTr69u2L/v37Y82aNbhw4QL++OMPxMXFYePGjQCA0aNHY8uWLbhw4QKOHDmCnTt3IigoCADg6+sLlUqFDRs24ObNm8jOzjY4lkcFBgZi69at2L9/P86cOYNhw4bhxo0b0nFHR0e8/fbb+PDDD7Fjxw6cPHkSAwcO1EnaAgMD0aVLF7zzzjvYu3cvjh8/jn79+qFy5cro0qWL0WIlMndMVIhecIGBgUhKSkK1atXQq1cv+Pv7Y+jQoWjTpg0OHDgAd3f3xz6+SZMmWLRoEebMmYN69eohISFBLwGJiIjAhg0bkJCQgIYNG6JJkyaYNWsWfH19nyn2xYsXo3///nj//fdRo0YNdO3aFYmJifDx8QHwoLYkKioKQUFBeO2111C9enV8++23AIDKlSsjNjYW48aNQ4UKFTBy5MhniuVhn376KRo0aICIiAi0bt0aFStWRNeuXXXKfPXVV2jZsiU6d+6M8PBwtGjRAiEhIXr3FxISgk6dOqFp06YQQmDTpk16TWREVDqVeJpeeERERETPEWtUiIiISLGYqBAREZFiMVEhIiIixWKiQkRERIrFRIWIiIgUi4kKERERKRYTFSIiIlIsJipERESkWExUiIiISLGYqBAREZFiMVEhIiIixWKiQkRERIr1f6BsHll9ObzzAAAAAElFTkSuQmCC",
|
|
"text/plain": [
|
|
"<Figure size 640x480 with 1 Axes>"
|
|
]
|
|
},
|
|
"metadata": {},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"import matplotlib.pyplot as plt\n",
|
|
"\n",
|
|
"# Parámetros de la simulación\n",
|
|
"iteraciones = 30\n",
|
|
"N = 500\n",
|
|
"orden_real = 3\n",
|
|
"varianza_ruido = 0.001\n",
|
|
"orden_max = 8\n",
|
|
"C = 5.0\n",
|
|
"\n",
|
|
"ordenes_estimados = []\n",
|
|
"aciertos = 0\n",
|
|
"\n",
|
|
"for _ in range(iteraciones):\n",
|
|
" # Generar señal sintética con 'orden_real' cosenos\n",
|
|
" x, _, _ = generar_muestra_coseno(N=N, num_cosenos=orden_real, varianza_ruido=varianza_ruido)\n",
|
|
" \n",
|
|
" # Estimar el orden\n",
|
|
" orden_est, _, _ = estimador_de_orden(x, orden_max=orden_max, C=C)\n",
|
|
" ordenes_estimados.append(orden_est)\n",
|
|
" \n",
|
|
" if orden_est == orden_real:\n",
|
|
" aciertos += 1\n",
|
|
"\n",
|
|
"# Resultados\n",
|
|
"print(\"Orden real:\", orden_real)\n",
|
|
"print(\"Tasa de acierto:\", aciertos / iteraciones)\n",
|
|
"print(\"Órdenes estimados:\", ordenes_estimados)\n",
|
|
"\n",
|
|
"# Histograma de órdenes estimados\n",
|
|
"plt.hist(ordenes_estimados, bins=range(1, orden_max+2), align=\"left\", rwidth=0.7)\n",
|
|
"plt.axvline(orden_real, color=\"red\", linestyle=\"--\", label=\"Orden real\")\n",
|
|
"plt.xlabel(\"Orden estimado\")\n",
|
|
"plt.ylabel(\"Frecuencia\")\n",
|
|
"plt.title(f\"Distribución de órdenes estimados en {iteraciones} corridas\")\n",
|
|
"plt.legend()\n",
|
|
"plt.show()\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"# Segundo intento del detector de orden\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Generador de muestras\n",
|
|
"\n",
|
|
"Esta sección se describe los cambios y mejoras implementadas en el generador de señales cardiorespiratorias sintéticas para simulaciones con radar UWB. La función ha sido adaptada para generar señales más realistas, controlables y útiles para métodos no lineales de análisis espectral.\n",
|
|
"\n",
|
|
"- Se modificaron los intervalos de frecuencia fundamentales para que reflejen rangos fisiológicos realistas para mediciones con radar UWB: la frecuencia cardíaca se estableció entre 0.6 y 1.1 Hz, lo que corresponde aproximadamente a 36 a 66 latidos por minuto, y la frecuencia respiratoria entre 0.15 y 0.35 Hz, equivalente a 9 a 21 respiraciones por minuto. Los centros de estos intervalos se utilizan como semillas para la inicialización de métodos no lineales, siendo 0.85 Hz para la cardíaca y 0.25 Hz para la respiratoria.\n",
|
|
"- Se definió que el número de componentes de cada señal corresponde al fundamental más los armónicos especificados por el usuario. Esto permite simular señales más complejas y con mayor riqueza espectral, manteniendo control sobre la cantidad de información armónica incluida.\n",
|
|
"- Los armónicos no son múltiplos exactos de la frecuencia fundamental. Se introduce un factor de dispersión aleatoria de ±5% para cada armónico, de modo que las frecuencias no sean exactos múltiplos enteros del fundamental, simulando así la inarmonicidad típica de señales fisiológicas reales.\n",
|
|
"- Las amplitudes de los componentes fueron ajustadas para reflejar la relación típica en señales UWB. La señal respiratoria es significativamente mayor que la cardíaca, con amplitudes generadas en el intervalo de 1.0 a 2.0 para respiratoria y de 0.01 a 0.1 para cardíaca, manteniendo la proporción de orden de magnitud observada en mediciones reales.\n",
|
|
"- Se mantiene la generación aleatoria de fases y offsets para cada componente dentro de rangos definidos, lo que asegura variabilidad y realismo en la forma de la señal.\n",
|
|
"- Se añadió la posibilidad de introducir ruido aditivo gaussiano con varianza especificada por el usuario. Además, la función calcula automáticamente la relación señal a ruido (SNR) resultante en decibelios, lo que permite evaluar la calidad de la señal generada en condiciones controladas.\n",
|
|
"- Los valores por defecto para la frecuencia de muestreo y la duración de la señal son 200 Hz y 30 segundos, respectivamente, parámetros adecuados para simulaciones de señales cardiorespiratorias con radar UWB.\n",
|
|
"- Se incluyen diccionarios de salida detallados: uno con los parámetros utilizados (frecuencias, amplitudes, fases, offsets, varianza y SNR), otro con los límites de cada parámetro para posibles optimizaciones, y otro con las semillas de inicialización. Además, se entrega metadata con información de la simulación, como número de componentes, frecuencia de muestreo, duración y factor de dispersión aplicado, facilitando la reutilización y análisis posterior de la señal generada."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 1,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"import numpy as np\n",
|
|
"\n",
|
|
"def generar_muestras_cardiorespiratorias(N_card=3, N_resp=2, fs=200, T=30, sigma2=0.01, epsilon=0.05, seed=None):\n",
|
|
" \"\"\"\n",
|
|
" Generador de señales cardiorespiratorias sintéticas para radar UWB.\n",
|
|
"\n",
|
|
" Args:\n",
|
|
" N_card (int): número de componentes de la señal cardíaca (fundamental + armónicos).\n",
|
|
" N_resp (int): número de componentes de la señal respiratoria (fundamental + armónicos).\n",
|
|
" fs (float): frecuencia de muestreo [Hz].\n",
|
|
" T (float): duración de la señal [s].\n",
|
|
" sigma2 (float): varianza del ruido aditivo gaussiano.\n",
|
|
" epsilon (float): factor de dispersión para armónicos (±5% por defecto).\n",
|
|
" seed (int or None): semilla para reproducibilidad.\n",
|
|
"\n",
|
|
" Returns:\n",
|
|
" signal (np.array): señal sintetizada con ruido.\n",
|
|
" params (dict): diccionario con parámetros de frecuencias, amplitudes, fases.\n",
|
|
" bounds (dict): límites de cada parámetro.\n",
|
|
" seeds (dict): centros de los intervalos para inicialización no lineal.\n",
|
|
" metadata (dict): información de la simulación.\n",
|
|
" \"\"\"\n",
|
|
"\n",
|
|
" if seed is not None:\n",
|
|
" np.random.seed(seed)\n",
|
|
"\n",
|
|
" t = np.arange(0, T, 1/fs)\n",
|
|
"\n",
|
|
" # Intervalos de frecuencias fundamentales\n",
|
|
" f_card_interval = np.array([0.6, 1.1])\n",
|
|
" f_resp_interval = np.array([0.15, 0.35])\n",
|
|
" f_card_center = f_card_interval.mean()\n",
|
|
" f_resp_center = f_resp_interval.mean()\n",
|
|
"\n",
|
|
" # Amplitudes\n",
|
|
" A_card_interval = np.array([0.01, 0.1])\n",
|
|
" A_resp_interval = np.array([1.0, 2.0])\n",
|
|
" A_card_center = A_card_interval.mean()\n",
|
|
" A_resp_center = A_resp_interval.mean()\n",
|
|
"\n",
|
|
" # Fases y offsets\n",
|
|
" phase_card = np.random.uniform(0, 2*np.pi, N_card)\n",
|
|
" phase_resp = np.random.uniform(0, 2*np.pi, N_resp)\n",
|
|
"\n",
|
|
" offset_card = np.random.uniform(-0.1,0.1)\n",
|
|
" offset_resp = np.random.uniform(-0.1,0.1)\n",
|
|
"\n",
|
|
" # Frecuencia fundamental\n",
|
|
" f_card0 = np.random.uniform(*f_card_interval)\n",
|
|
" f_resp0 = np.random.uniform(*f_resp_interval)\n",
|
|
"\n",
|
|
" # Generar armónicos con dispersión\n",
|
|
" f_card = [f_card0]\n",
|
|
" for n in range(2, N_card+1):\n",
|
|
" delta = np.random.uniform(-epsilon, epsilon)\n",
|
|
" f_card.append(n*f_card0*(1+delta))\n",
|
|
"\n",
|
|
" f_resp = [f_resp0]\n",
|
|
" for n in range(2, N_resp+1):\n",
|
|
" delta = np.random.uniform(-epsilon, epsilon)\n",
|
|
" f_resp.append(n*f_resp0*(1+delta))\n",
|
|
"\n",
|
|
" f_card = np.array(f_card)\n",
|
|
" f_resp = np.array(f_resp)\n",
|
|
"\n",
|
|
" # Amplitudes para cada componente\n",
|
|
" A_card = np.random.uniform(A_card_interval[0], A_card_interval[1], N_card)\n",
|
|
" A_resp = np.random.uniform(A_resp_interval[0], A_resp_interval[1], N_resp)\n",
|
|
"\n",
|
|
" # Generar señal\n",
|
|
" signal_card = np.sum([A_card[i]*np.cos(2*np.pi*f_card[i]*t + phase_card[i]) for i in range(N_card)], axis=0) + offset_card\n",
|
|
" signal_resp = np.sum([A_resp[i]*np.cos(2*np.pi*f_resp[i]*t + phase_resp[i]) for i in range(N_resp)], axis=0) + offset_resp\n",
|
|
"\n",
|
|
" signal_clean = signal_card + signal_resp\n",
|
|
"\n",
|
|
" # Ruido y SNR\n",
|
|
" noise = np.random.normal(0, np.sqrt(sigma2), len(t))\n",
|
|
" signal_noisy = signal_clean + noise\n",
|
|
"\n",
|
|
" snr = 10*np.log10(np.var(signal_clean)/sigma2)\n",
|
|
"\n",
|
|
" # Preparar diccionarios de salida\n",
|
|
" params = {\n",
|
|
" 'f_card': f_card,\n",
|
|
" 'f_resp': f_resp,\n",
|
|
" 'A_card': A_card,\n",
|
|
" 'A_resp': A_resp,\n",
|
|
" 'phase_card': phase_card,\n",
|
|
" 'phase_resp': phase_resp,\n",
|
|
" 'offset_card': offset_card,\n",
|
|
" 'offset_resp': offset_resp,\n",
|
|
" 'sigma2': sigma2,\n",
|
|
" 'SNR_dB': snr\n",
|
|
" }\n",
|
|
"\n",
|
|
" bounds = {\n",
|
|
" 'f_card': f_card_interval,\n",
|
|
" 'f_resp': f_resp_interval,\n",
|
|
" 'A_card': A_card_interval,\n",
|
|
" 'A_resp': A_resp_interval,\n",
|
|
" 'phase': [0, 2*np.pi],\n",
|
|
" 'offset': [-0.1, 0.1]\n",
|
|
" }\n",
|
|
"\n",
|
|
" seeds = {\n",
|
|
" 'f_card_center': f_card_center,\n",
|
|
" 'f_resp_center': f_resp_center,\n",
|
|
" 'A_card_center': A_card_center,\n",
|
|
" 'A_resp_center': A_resp_center\n",
|
|
" }\n",
|
|
"\n",
|
|
" metadata = {\n",
|
|
" 'N_card': N_card,\n",
|
|
" 'N_resp': N_resp,\n",
|
|
" 'fs': fs,\n",
|
|
" 'T': T,\n",
|
|
" 'epsilon': epsilon\n",
|
|
" }\n",
|
|
"\n",
|
|
" return signal_noisy, params, bounds, seeds, metadata"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Estimador de orden"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 4,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"import numpy as np\n",
|
|
"from scipy.optimize import minimize\n",
|
|
"import matplotlib.pyplot as plt\n",
|
|
"\n",
|
|
"def estimar_parametros(y, t, N, rangos_freq):\n",
|
|
" \"\"\"\n",
|
|
" Estima parámetros de una suma de N cosenos:\n",
|
|
" y(t) = sum A_k cos(2π f_k t + φ_k)\n",
|
|
"\n",
|
|
" Estrategia:\n",
|
|
" - No lineal: optimiza solo las frecuencias.\n",
|
|
" - Lineal: estima amplitudes y fases para cada conjunto de frecuencias.\n",
|
|
"\n",
|
|
" Parámetros\n",
|
|
" ----------\n",
|
|
" y : array_like\n",
|
|
" Señal observada (vector de tamaño M).\n",
|
|
" t : array_like\n",
|
|
" Vector de tiempos (tamaño M).\n",
|
|
" N : int\n",
|
|
" Número de cosenos a estimar.\n",
|
|
" rangos_freq : list of tuples\n",
|
|
" Lista [(fmin1,fmax1), (fmin2,fmax2), ...] con los rangos válidos\n",
|
|
" para cada frecuencia.\n",
|
|
"\n",
|
|
" Retorna\n",
|
|
" -------\n",
|
|
" params : list of tuples\n",
|
|
" Lista [(A1,f1,phi1), ..., (AN,fN,phiN)].\n",
|
|
" \"\"\"\n",
|
|
"\n",
|
|
" M = len(y)\n",
|
|
"\n",
|
|
" # --- Función de coste: dado un vector de frecuencias, ajusta amplitudes/fases y calcula error ---\n",
|
|
" def coste(frecuencias):\n",
|
|
" # Matriz de diseño con cosenos y senos\n",
|
|
" X = []\n",
|
|
" for f in frecuencias:\n",
|
|
" X.append(np.cos(2*np.pi*f*t))\n",
|
|
" X.append(np.sin(2*np.pi*f*t))\n",
|
|
" X = np.column_stack(X) # M x (2N)\n",
|
|
"\n",
|
|
" # Resolver mínimos cuadrados lineales\n",
|
|
" coef, _, _, _ = np.linalg.lstsq(X, y, rcond=None)\n",
|
|
" y_hat = X @ coef\n",
|
|
" resid = y - y_hat\n",
|
|
" return np.sum(resid**2)\n",
|
|
"\n",
|
|
" # --- Inicialización aleatoria dentro de los rangos ---\n",
|
|
" f0 = np.array([np.random.uniform(low, high) for (low, high) in rangos_freq])\n",
|
|
"\n",
|
|
" # --- Restricciones de cada frecuencia ---\n",
|
|
" bounds = rangos_freq\n",
|
|
"\n",
|
|
" # --- Optimización no lineal en frecuencias ---\n",
|
|
" res = minimize(coste, f0, bounds=bounds, method=\"L-BFGS-B\")\n",
|
|
"\n",
|
|
" f_est = res.x\n",
|
|
"\n",
|
|
" # --- Una vez obtenidas las frecuencias, estimar amplitudes y fases ---\n",
|
|
" X = []\n",
|
|
" for f in f_est:\n",
|
|
" X.append(np.cos(2*np.pi*f*t))\n",
|
|
" X.append(np.sin(2*np.pi*f*t))\n",
|
|
" X = np.column_stack(X)\n",
|
|
"\n",
|
|
" coef, _, _, _ = np.linalg.lstsq(X, y, rcond=None)\n",
|
|
"\n",
|
|
" params = []\n",
|
|
" for k in range(N):\n",
|
|
" a = coef[2*k]\n",
|
|
" b = coef[2*k + 1]\n",
|
|
" A = np.sqrt(a**2 + b**2)\n",
|
|
" phi = np.arctan2(-b, a)\n",
|
|
" params.append((A, f_est[k], phi))\n",
|
|
"\n",
|
|
" return params\n",
|
|
"\n",
|
|
"\n",
|
|
"def estimar_orden(y, t, N_max, C, rangos_freq, eps=1e-8):\n",
|
|
" \"\"\"\n",
|
|
" Estima el orden de la señal (cantidad de cosenos) probando N=1..N_max\n",
|
|
" con el criterio penalizado.\n",
|
|
" \"\"\"\n",
|
|
"\n",
|
|
" M = len(y)\n",
|
|
" criterios = []\n",
|
|
"\n",
|
|
" for N in range(1, N_max+1):\n",
|
|
" # Usamos solo los primeros N rangos\n",
|
|
" params = estimar_parametros(y, t, N, rangos_freq[:N])\n",
|
|
"\n",
|
|
" # Reconstrucción\n",
|
|
" y_hat = np.zeros_like(y)\n",
|
|
" for A, f, phi in params:\n",
|
|
" y_hat += A * np.cos(2*np.pi*f*t + phi)\n",
|
|
"\n",
|
|
" # Varianza del error\n",
|
|
" resid = y - y_hat\n",
|
|
" sigma2 = np.mean(resid**2)\n",
|
|
"\n",
|
|
" # Criterio penalizado (según tu pedido)\n",
|
|
" J = M * np.log(sigma2 + eps) + C * N\n",
|
|
" criterios.append(J)\n",
|
|
"\n",
|
|
" criterios = np.array(criterios)\n",
|
|
" N_hat = np.argmin(criterios) + 1\n",
|
|
"\n",
|
|
" return N_hat, criterios\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 6,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"ename": "IndexError",
|
|
"evalue": "index 10 is out of bounds for axis 0 with size 10",
|
|
"output_type": "error",
|
|
"traceback": [
|
|
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
|
|
"\u001b[1;31mIndexError\u001b[0m Traceback (most recent call last)",
|
|
"Cell \u001b[1;32mIn[6], line 38\u001b[0m\n\u001b[0;32m 33\u001b[0m params_est \u001b[38;5;241m=\u001b[39m estimar_parametros(signal_noisy, t, N_true, rangos_freq)\n\u001b[0;32m 35\u001b[0m \u001b[38;5;66;03m# ============================\u001b[39;00m\n\u001b[0;32m 36\u001b[0m \u001b[38;5;66;03m# 3. Estimar orden\u001b[39;00m\n\u001b[0;32m 37\u001b[0m \u001b[38;5;66;03m# ============================\u001b[39;00m\n\u001b[1;32m---> 38\u001b[0m N_hat, criterios \u001b[38;5;241m=\u001b[39m \u001b[43mestimar_orden\u001b[49m\u001b[43m(\u001b[49m\u001b[43msignal_noisy\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mt\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mN_max\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mN_max\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mC\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mC\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mrangos_freq\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mrangos_freq\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 40\u001b[0m \u001b[38;5;66;03m# ============================\u001b[39;00m\n\u001b[0;32m 41\u001b[0m \u001b[38;5;66;03m# 4. Reconstrucción de la señal\u001b[39;00m\n\u001b[0;32m 42\u001b[0m \u001b[38;5;66;03m# ============================\u001b[39;00m\n\u001b[0;32m 43\u001b[0m y_hat \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mzeros_like(signal_noisy)\n",
|
|
"Cell \u001b[1;32mIn[4], line 91\u001b[0m, in \u001b[0;36mestimar_orden\u001b[1;34m(y, t, N_max, C, rangos_freq, eps)\u001b[0m\n\u001b[0;32m 87\u001b[0m criterios \u001b[38;5;241m=\u001b[39m []\n\u001b[0;32m 89\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m N \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mrange\u001b[39m(\u001b[38;5;241m1\u001b[39m, N_max\u001b[38;5;241m+\u001b[39m\u001b[38;5;241m1\u001b[39m):\n\u001b[0;32m 90\u001b[0m \u001b[38;5;66;03m# Usamos solo los primeros N rangos\u001b[39;00m\n\u001b[1;32m---> 91\u001b[0m params \u001b[38;5;241m=\u001b[39m \u001b[43mestimar_parametros\u001b[49m\u001b[43m(\u001b[49m\u001b[43my\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mt\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mN\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mrangos_freq\u001b[49m\u001b[43m[\u001b[49m\u001b[43m:\u001b[49m\u001b[43mN\u001b[49m\u001b[43m]\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 93\u001b[0m \u001b[38;5;66;03m# Reconstrucción\u001b[39;00m\n\u001b[0;32m 94\u001b[0m y_hat \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mzeros_like(y)\n",
|
|
"Cell \u001b[1;32mIn[4], line 71\u001b[0m, in \u001b[0;36mestimar_parametros\u001b[1;34m(y, t, N, rangos_freq)\u001b[0m\n\u001b[0;32m 69\u001b[0m params \u001b[38;5;241m=\u001b[39m []\n\u001b[0;32m 70\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m k \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mrange\u001b[39m(N):\n\u001b[1;32m---> 71\u001b[0m a \u001b[38;5;241m=\u001b[39m \u001b[43mcoef\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;241;43m2\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mk\u001b[49m\u001b[43m]\u001b[49m\n\u001b[0;32m 72\u001b[0m b \u001b[38;5;241m=\u001b[39m coef[\u001b[38;5;241m2\u001b[39m\u001b[38;5;241m*\u001b[39mk \u001b[38;5;241m+\u001b[39m \u001b[38;5;241m1\u001b[39m]\n\u001b[0;32m 73\u001b[0m A \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39msqrt(a\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m2\u001b[39m \u001b[38;5;241m+\u001b[39m b\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m2\u001b[39m)\n",
|
|
"\u001b[1;31mIndexError\u001b[0m: index 10 is out of bounds for axis 0 with size 10"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"\n",
|
|
"# ============================\n",
|
|
"# Configuración del experimento\n",
|
|
"# ============================\n",
|
|
"fs = 200\n",
|
|
"T = 30\n",
|
|
"sigma2 = 0.01\n",
|
|
"N_card_true = 3\n",
|
|
"N_resp_true = 2\n",
|
|
"C = 10 # constante de penalización\n",
|
|
"N_max = 8 # orden máximo a probar\n",
|
|
"\n",
|
|
"# ============================\n",
|
|
"# 1. Generar señal sintética\n",
|
|
"# ============================\n",
|
|
"signal_noisy, params_true, bounds, seeds, metadata = generar_muestras_cardiorespiratorias(\n",
|
|
" N_card=N_card_true,\n",
|
|
" N_resp=N_resp_true,\n",
|
|
" fs=fs,\n",
|
|
" T=T,\n",
|
|
" sigma2=sigma2,\n",
|
|
" seed=42\n",
|
|
")\n",
|
|
"\n",
|
|
"t = np.arange(0, T, 1/fs)\n",
|
|
"\n",
|
|
"# ============================\n",
|
|
"# 2. Estimar parámetros para orden verdadero\n",
|
|
"# ============================\n",
|
|
"# Aquí se asume que ya tenés `estimar_parametros` definida\n",
|
|
"N_true = N_card_true + N_resp_true\n",
|
|
"rangos_freq = [bounds['f_card']] * N_card_true + [bounds['f_resp']] * N_resp_true\n",
|
|
"\n",
|
|
"params_est = estimar_parametros(signal_noisy, t, N_true, rangos_freq)\n",
|
|
"\n",
|
|
"# ============================\n",
|
|
"# 3. Estimar orden\n",
|
|
"# ============================\n",
|
|
"N_hat, criterios = estimar_orden(signal_noisy, t, N_max=N_max, C=C, rangos_freq=rangos_freq)\n",
|
|
"\n",
|
|
"# ============================\n",
|
|
"# 4. Reconstrucción de la señal\n",
|
|
"# ============================\n",
|
|
"y_hat = np.zeros_like(signal_noisy)\n",
|
|
"for A, f, phi in params_est:\n",
|
|
" y_hat += A * np.cos(2*np.pi*f*t + phi)\n",
|
|
"\n",
|
|
"# ============================\n",
|
|
"# 5. Gráficos\n",
|
|
"# ============================\n",
|
|
"plt.figure(figsize=(12,5))\n",
|
|
"plt.plot(t, signal_noisy, label=\"Señal observada\", alpha=0.7)\n",
|
|
"plt.plot(t, y_hat, label=f\"Reconstrucción (N={N_true})\", linewidth=2)\n",
|
|
"plt.xlabel(\"Tiempo [s]\")\n",
|
|
"plt.ylabel(\"Amplitud\")\n",
|
|
"plt.legend()\n",
|
|
"plt.title(\"Señal sintética vs Reconstrucción\")\n",
|
|
"plt.show()\n",
|
|
"\n",
|
|
"plt.figure(figsize=(6,4))\n",
|
|
"plt.plot(range(1, len(criterios)+1), criterios, marker=\"o\")\n",
|
|
"plt.axvline(N_hat, color=\"r\", linestyle=\"--\", label=f\"Orden estimado = {N_hat}\")\n",
|
|
"plt.xlabel(\"Orden N\")\n",
|
|
"plt.ylabel(\"Criterio penalizado\")\n",
|
|
"plt.title(\"Selección de orden\")\n",
|
|
"plt.legend()\n",
|
|
"plt.show()\n",
|
|
"\n",
|
|
"# ============================\n",
|
|
"# 6. Resultados\n",
|
|
"# ============================\n",
|
|
"print(\"Parámetros verdaderos:\")\n",
|
|
"for i, (A,f,phi) in enumerate(zip(\n",
|
|
" np.concatenate([params_true['A_card'], params_true['A_resp']]),\n",
|
|
" np.concatenate([params_true['f_card'], params_true['f_resp']]),\n",
|
|
" np.concatenate([params_true['phase_card'], params_true['phase_resp']])\n",
|
|
" )):\n",
|
|
" print(f\"Componente {i+1}: A={A:.3f}, f={f:.3f}, phi={phi:.3f}\")\n",
|
|
"\n",
|
|
"print(\"\\nParámetros estimados:\")\n",
|
|
"for i, (A,f,phi) in enumerate(params_est):\n",
|
|
" print(f\"Componente {i+1}: A={A:.3f}, f={f:.3f}, phi={phi:.3f}\")\n",
|
|
"\n",
|
|
"print(f\"\\nOrden verdadero: {N_true}\")\n",
|
|
"print(f\"Orden estimado : {N_hat}\")"
|
|
]
|
|
}
|
|
],
|
|
"metadata": {
|
|
"kernelspec": {
|
|
"display_name": "Python 3",
|
|
"language": "python",
|
|
"name": "python3"
|
|
},
|
|
"language_info": {
|
|
"codemirror_mode": {
|
|
"name": "ipython",
|
|
"version": 3
|
|
},
|
|
"file_extension": ".py",
|
|
"mimetype": "text/x-python",
|
|
"name": "python",
|
|
"nbconvert_exporter": "python",
|
|
"pygments_lexer": "ipython3",
|
|
"version": "3.12.6"
|
|
}
|
|
},
|
|
"nbformat": 4,
|
|
"nbformat_minor": 2
|
|
}
|