diff --git a/Simulaciones/Detector_Orden.ipynb b/Simulaciones/Detector_Orden.ipynb index 7adefe9..a2374e7 100644 --- a/Simulaciones/Detector_Orden.ipynb +++ b/Simulaciones/Detector_Orden.ipynb @@ -368,7 +368,7 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 1, "metadata": {}, "outputs": [], "source": [ @@ -493,6 +493,237 @@ "\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": {