prepping for pr

examples/fourier.ipynb

@@ -49,7 +49,7 @@
     "\n",
     "### FS - Fourier Series\n",
     "$$\n",
-    "X(k) = \\int_{0}^{P}\n",
+    "X(k) = \\int_{0}^{P} x(n) e^{-j \\omega_k n}\n",
     "$$"
    ]
   }

examples/frequency_estimators.ipynb

@@ -0,0 +1,57 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import os, sys\n",
+    "from typing import NoReturn\n",
+    "import logging\n",
+    "\n",
+    "import numpy as np\n",
+    "from numpy.typing import ArrayLike\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "import scipy as sp\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dir = os.path.abspath('../ssp')\n",
+    "dir = os.path.dirname(dir)\n",
+    "sys.path.append(dir)\n",
+    "\n",
+    "from ssp import spectrum\n",
+    "\n"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": ".env.py.312",
+   "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.1"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

ssp/spectrum.m

@@ -5,4 +5,4 @@
   sigma = min(diag(d);
   index = find(diag(d) == sigma);
   vmin = v(:, index);
-  end;
+  end;x

ssp/spectrum.py

@@ -190,25 +190,25 @@ def music(x: ArrayLike, p: int, M: int) -> ArrayLike:
 
     Page 430, Figure 8.19
     """
-    _x = np.array(x).ravel().reshape(-1, 1)
+    _x = np.array(x, ndmin=1)
     if p + 1 > M or len(_x) < M:
         raise ValueError("Size of signal covariance matrix is inappropriate.")
 
     R = covar(x, M)
     d, v = np.linalg.eig(R)
-    ddiag = np.diag(d)
-    i = np.argsort(ddiag)
-    # y = ddiag[i]
+    y = np.sort(d)
+    i = np.argsort(d)
     nfft = max(len(_x) + 1, 1024)
     Px = np.zeros(nfft, dtype=complex)
 
     for j in range(M - p):
         Px = Px + np.abs(np.fft.fft(v[:, i[j]], nfft))
 
+    Ax = Px.copy()
     Px = -20 * np.log10(Px)
     logger.debug(f'{Px.shape=}')
 
-    return Px
+    return Px, Ax
 
 
 def ev(x: ArrayLike, p: int, M: int) -> np.ndarray:

tests/test_spectrum.py

@@ -15,7 +15,7 @@
 
 def test_phd() -> None:
     # two complex exponentials in white noise have the following autocorrelation sequence
-    rx = np.array([6, 1.92705 + 4.58522j, -3.42705 + 3.49541j], dtype=complex)
+    rx = np.array([6, 1.92705 + 4.58522j, -3.42705 + 3.49541j], dtype=complex) / 6
     logger.warning(f'{rx=}')