Prony's method and matrix pencil method performance on determining the complex natural resonance frequencies of a linear system

==Resumen==

La caracterización de un sistema físico es un asunto de suma importancia para abordar algunos problemas de la física aplicada y de la ingeniería. Las frecuencias complejas de resonancia de un sistema, las cuales están incluidas en su respuesta impulsiva son parte característica de dicho sistema y ayudan a describirlo. Son pocos los trabajos escritos en idioma inglés que muestran una comparación entre los métodos discretos que extraen las frecuencias complejas de resonancia a través del estudio de la respuesta impulsiva del sistema que se esté estudiando. Mucho menos común es encontrar este tipo de estudios comparativos escritos en idioma castellano. Dada esta situación, en este trabajo de investigación se describen, prueban y comparan algunos de los métodos numéricos discretos más importantes para estimar las frecuencias naturales complejas de resonancia de un sistema, a través del estudio de su respuesta impulsiva en diferentes escenarios de simulación. De acuerdo con los resultados obtenidos, el método del haz de matrices con filtro SVD es el método menos sensible al ruido, mientras que el método Prony y sus diferentes versiones son los más eficientes desde el punto de vista computacional. Se discuten los escenarios de simulación que podrían ser más adecuados para cada método.

'''Palabras clave''': Resonancias naturales complejas, método del haz de matrices, método de Prony, descomposición en valores singulares, método de mínimos cuadrados totales

==Abstract==

Characterization of a physical system is an important issue to approach some applied physics and engineering problems. The complex natural resonance frequencies of the system which are included in its impulsive response are characteristic of such system and are part of its description. Few works written in English language show a comparison among discrete methods that extract natural complex frequencies from a system impulsive response. Much less common is to find works written in Spanish language about this important research topic. Given this situation, important discrete numeric methods to estimate the complex natural resonance frequencies of a system through its impulsive response are described, tested and compared in different simulation scenarios in this document. According to the obtained results, the matrix pencil method with a SVD filter is the less sensitive method to the noise, while the Prony method and its different versions are the fastest ones. Scenarios that could be more suitable for each method are discussed.

'''Keywords''': Complex Natural Resonances (CNRs), Matrix Pencil (MP) method, Prony's method, Singular Value Decomposition (SVD), Total Least Square (TLS) method

==1. Introducción==

El estudio de un sistema físico a partir de su respuesta impulsiva es un asunto relevante y vigente en el mundo de la física aplicada y la ingeniería. La respuesta impulsiva de un sistema lineal invariante en el tiempo contiene las denominadas frecuencias naturales complejas de resonancia (FCNR), distintivas de dicho sistema, a partir de las cuales el sistema puede ser identificado sin necesidad de conocer ''a priori'' su constitución interna.

Entre las técnicas más utilizadas para estimar las FCNR de un sistema dado (en el dominio del tiempo) se reconocen dos grandes métodos, a saber: el método de Prony y el método de haz de matrices, los cuales son algoritmos con un fundamento teórico que no está basado en el análisis de Fourier. Vale la pena mencionar que en el dominio de la frecuencia existen otros métodos, ver por ejemplo <span id='citeF-1'></span>[[#cite-1|[1]]], los cuales, sin embargo, no serán motivo de estudio en este trabajo. La versión original del método de Prony fue propuesta por Gaspard Riche de Prony en 1795 <span id='citeF-2'></span>[[#cite-2|[2]]-<span id='citeF-3'></span>[[#cite-3|3]]] y la del método del haz de matrices por Y. Hua y T.K. Sarkar en 1989 <span id='citeF-4'></span>[[#cite-4|[4]]]. A partir de entonces y hasta la fecha de hoy, se han planteado nuevas versiones de ambos métodos, algunas de las cuales han demostrado ser más robustas en presencia de ruido <span id='citeF-5'></span>[[#cite-5|[5-9]]].

Si bien el método de Prony fue concebido originalmente para tratar de modelar el fenómeno de expansión de los gases y el método del haz de matrices recibió al principio una atención especial en el campo de la tecnología de radar para la identificación de blancos no amigables en un contexto militar, en el proceso de revisión bibliográfica realizada durante el desarrollo de este trabajo se encontró una gran utilidad práctica de las distintas versiones de ambos métodos, para aplicaciones tan disímiles como las siguientes: determinación de las frecuencias naturales de oscilación de una red eléctrica <span id='citeF-10'></span>[[#cite-10|[10]]], detección de fallas y monitoreo de armónicos de la red de suministro de potencia eléctrica <span id='citeF-11'></span><span id='citeF-12'></span><span id='citeF-13'></span>[[#cite-11|[11-13]]], monitoreo de armónicos de convertidores electrónicos de sistemas de potencia eléctrica <span id='citeF-22'></span>[[#cite-14|[14]]], estudio de los sistemas de control <span id='citeF-15'></span>[[#cite-15|[15]]], procesamiento de señales de radar <span id='citeF-16'></span><span id='citeF-17'></span><span id='citeF-18'></span><span id='citeF-19'></span><span id='citeF-20'></span><span id='citeF-21'></span>[[#cite-16|[16-21]]], detección de armas de mano escondidas en la ropa <span id='citeF-22'>[[#cite-22|[22]]-</span><span id='citeF-23'></span>[[#cite-23|23]]], diseño y análisis de metamateriales <span id='citeF-24'>[[#cite-24|[24]]-</span><span id='citeF-25'></span>[[#cite-25|25]]], diseño de líneas de transmisión y análisis de antenas de microondas <span id='citeF-26'></span>[[#cite-26|[26]]-<span id='citeF-27'></span>[[#cite-27|27]]], monitoreo y diagnóstico biomédico <span id='citeF-28'></span>[[#cite-28|[28-32]]], procesamiento de imágenes de resonancia magnética nuclear <span id='citeF-33'></span>[[#cite-33|[33-35]]], microscopía de superresolución <span id='citeF-36'></span>[[#cite-36|[36]]], medición no invasiva de glucosa en sangre mediante microondas <span id='citeF-37'></span>[[#cite-37|[37]]], así como análisis y procesamiento de señales provenientes de sistemas de audio<span id='citeF-38'></span>[[#cite-38|[38]]], telescopios astronómicos <span id='citeF-39'></span>[[#cite-39|[39]]], detectores de ondas gravitacionales <span id='citeF-40'></span>[[#cite-40|[40]]], sonar <span id='citeF-41'></span>[[#cite-41|[41-43]]], sismógrafos <span id='citeF-44'></span>[[#cite-44|[44]]] y medidores de fuerzas de Coriolis <span id='citeF-45'></span>[[#cite-45|[45]]].

Vale la pena mencionar las aplicaciones que realizan el proceso de detección e identificación de un objeto en calidad de dispersor electromagnético a partir de su respuesta impulsiva <span id='citeF-16'></span>[[#cite-16|[16-25]]], ya que los métodos numéricos estudiados en este trabajo son los que comunmente son utilizados para dichos procesos de detección e identificación. Para estas aplicaciones, el conjunto de elementos físicos conformados por el objeto dispersor, la señal impulsiva que lo excita y su respectiva respuesta es lo que de ahora en adelante llamaremos sistema lineal. Visto el dispersor como un sistema, su respuesta impulsiva tiene como rasgo distintivo estar compuesta por un conjunto de FCNR dadas por <math>s_n=-\sigma _n+\jmath \omega _n</math>, las cuales permiten, entre otras cosas, identificar unívocamente al propio objeto y en conjunto con cierto conocimiento previo, caracterizarlo desde un punto de vista físico. Las FCNR propias de un objeto dispersor electromagnético se correlacionan directamente con su geometría y constitución material y se distribuyen sobre el plano complejo formando “ramas” simétricas respecto del eje real con cierto orden y separación que es también característico del objeto. Los objetos “altamente” resonantes, como por ejemplo un dipolo metálico, poseen un cojunto de FCNR que se alojan en ramas muy próximas al eje imaginario, alejándose de éste a medida que aumenta la frecuencia. Las FCNR de los objetos “suavemente” resonantes, en cambio, como por ejemplo una esfera metálica, se alojan en ramas que se separan mucho más notablemente del eje imaginario <span id='citeF-46'></span>[[#cite-46|[46]]]. La extracción de las componentes imaginarias <math>\omega _n</math> de estas FCNR no puede realizarse con precisión usando el análisis de Fourier por que al seccionar la función del sistema <math>H(s)</math> con el plano <math>\sigma =0</math> para obtener a <math>H(\omega )</math>, las <math>\omega _n</math> pueden quedar “escondidas” o ''aliased'', lo cual se observa en mayor medida para sistemas suavemente resonantes en comparación con los altamente resonantes <span id='citeF-29'></span>[[#cite-29|[29]]].

También vale la pena agregar que por lo general el conjunto de valores de <math display="inline">\omega _n</math> no son múltiplos de una frecuencia armónica natural fundamental como suele indicar la teoría de Fourier para señales periódicas, sino que más bien están distribuidas de manera arbitraria en el espectro, donde su valor es característico del sistema que se esté estudiando.

Aún dado el gran auge del uso del método de Prony, del método del haz de matrices y de sus variantes, no abundan las publicaciones arbitradas acerca del desempeño comparativo de dichos métodos en idioma castellano. Durante el proceso de revisión bibliográfica, en idioma castellano solo se consiguió el contexto histórico en el cual se crea el método de Prony, su descripción y su comparación con el método de la descomposición armónica de Pisarenko <span id='citeF-47'></span>[[#cite-47|[47]]]. También se encontraron en idioma castellano algunas publicaciones en las que se aplica el método de Prony para modelar una falla en una central termonuclear <span id='citeF-48'></span>[[#cite-48|[48]]], para la determinación de las frecuencias naturales de oscilación de una red eléctrica de potencia <span id='citeF-49'></span>[[#cite-49|[49]]-<span id='citeF-50'></span>[[#cite-49|50]]], para el procesamiento de imágenes de alta resolución <span id='citeF-51'></span>[[#cite-51|[51]]] y para la creación de un software para la predicción del valor de acciones financieras <span id='citeF-52'></span>[[#cite-52|[52]]].

Incluso en idioma inglés, son escasos los estudios comparativos publicados en revistas arbitradas acerca del desempeño de cada una de las variantes del método de Prony y del método del haz de matrices en diferentes escenarios de simulación. Entre los estudios comparativos encontrados durante el trabajo de revisión pueden citarse <span id='citeF-30'></span>[[#cite-30|[30]],<span id='citeF-40'></span>[[#cite-40|40]],<span id='citeF-53'></span>[[#cite-53|53-55]]], en los cuales se examinan solamente alguna de las versiones de los métodos de Prony y del haz de matrices en algunos escenarios de simulación particulares. Entre estos estudios, uno de los más detallados que se encontró fue el descrito en <span id='citeF-30'></span>[[#cite-30|[30]]], el cual muestra el desempeño de varios métodos generalmente utilizados en sistemas de análisis de señales bio-médicas, pero probando solo una de las versiones del método de Prony y del método del haz de matrices, sin tomar en cuenta el tiempo de cálculo de cada algoritmo, sin variar la frecuencia de muestreo y además variando los niveles de ruido en un rango no tan amplio como en el trabajo de investigación realizado en nuestro artículo. Otro trabajo de relevancia, de alta rigurosidad y enfoque estadístico, en el cual se muestra el desempeño de varias variantes del método de Prony y una variante del haz de matrices para diferentes números de muestras, sin variar la frecuencia de muestreo y tampoco el nivel de ruido presente en las señales que se estudian, se reporta en <span id='citeF-40'></span>[[#cite-40|[40]]].

Vale la pena agregar que también durante el proceso de revisión bibliográfica se encontró un importante estudio analítico no comparativo del desempeño particular del método de Prony <span id='citeF-56'></span>[[#cite-56|[56]]], en el cual se estudia la precisión teórica de dicho método en presencia de ruido.

Con base en el contexto anterior, con la intención de realizar un pequeño aporte al conocimiento y documentarlo en idioma castellano nos planteamos describir, analizar y comparar mediante simulación el desempeño de cinco métodos derivados de las teorías originales de Prony y del haz de matrices, a saber: el método de Prony clásico, el método de Prony con filtrado a través de la descomposición de valores singulares (DVS), el método de Prony usando mínimos cuadrados totales, el método del haz de matrices sin filtro y el método del haz de matrices con filtrado usando DVS. Para ello ponemos a prueba cada método en dos escenarios distintos de simulación, uno en el que se usan distintos niveles de ruido añadido y otro en el que se prueban distintas frecuencias de muestreo, ambos para una ventana de observación constante. En cada caso se registran los tiempos de cómputo de cada método. Luego, a partir de este trabajo comparativo, se determinan los escenarios más favorables para cada método, con los fines de proveer información valiosa que puede servir para desarrollar aplicaciones en los campos de la física y la ingeniería en general.

==2. Breve descripción de los cinco métodos estudiados.==

En este apartado describiremos brevemente las prescripciones de los métodos de Prony clásico, de haz de matrices y sus derivados. Específicamente revisaremos los fundamentos matemáticos de los métodos de Prony, de Prony con filtrado de valores singulares, de Prony con mínimos cuadrados, de haz de matrices y de haz de matrices con filtrado de valores singulares. Todos estos métodos se implementan en el dominio del tiempo y permiten estimar las frecuencias naturales complejas de resonancia de un sistema lineal invariante en el tiempo (LIT). Las frecuencias naturales complejas de un sistema LIT coinciden con los polos de la respuesta <math display="inline">H(s)</math> del sistema.

La respuesta impulsiva <math>h(t)</math> de un sistema lineal está compuesta por una combinación lineal de <math>N'</math> cosenos amortiguados de la forma

<span id="eq-1"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> h(t)=\sum _{n=1}^{N'} |r_n|e^{-\sigma _nt}\cos (\omega _nt+\phi _n)=\sum _{n=1}^{N'}\Re \left\{r_ne^{s_nt} \right\} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}

que a su vez forma un par transformado con la función <math>H(s)</math> del sistema:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> h(t)=\sum _{n=1}^{N'} \Re \left\{r_ne^{s_nt} \right\}\Leftrightarrow </math> <math>\sum _{n=1}^{N'}\frac{1}{2}\left(\frac{r_n}{s-s_n}+\frac{r_n^\ast }{s-s_n^\ast }\right)=H(s) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

donde <math>N'=N/2</math>, siendo <math>N</math> el número par de frecuencias naturales complejas de resonancia del sistema, <math display="inline">r_n</math> y <math display="inline">r_n^\ast </math>, con <math display="inline">r_n=|r_n|e^{\jmath \phi _n}</math>, son residuos complejos, y  <math display="inline">s_n</math> y <math display="inline">s_n^\ast </math> son los polos complejos del sistema también llamados frecuencias complejas naturales de resonancia, con <math>s_n=-\sigma _n+\jmath \omega _n</math> y <math>\sigma _n>0</math>, donde el símbolo <math display="inline">\ast </math> denota conjugación.

Los métodos de Prony y de haz de matrices en todas sus versiones, son métodos de estimación mediante ajuste numérico que parten de un número finito, dígase <math display="inline">M</math>, de muestras de <math display="inline">h(t)</math> tal cual se muestra en la <span id='img-1'></span>[[#img-1|Figura 1]].

<div id='img-1'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"| [[Image:Draft_Barroso_243065238-sistema2.png|360px|Sistema de reconstrucción de la respuesta impulsiva de un sistema lineal.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding:10px;"| '''Figura 1'''. Sistema de reconstrucción de la respuesta impulsiva de un sistema lineal
|}


En la <span id='img-1'></span>[[#img-1|Figura 1]], el sistema lineal <math>H(s)</math> es estimulado con un impulso <math>\delta (t)</math> ideal. Como resultado, el sistema responde produciendo a la salida una respuesta impulsiva de la forma dada por la ecuación <span id='eq-1'></span>[[#eq-1|(1)]]. Luego, la respuesta <math>h(t)</math> del sistema es muestreada a una tasa <math>T_s</math> uniforme: <math>h_{m}=h(mT_s)</math>. Las muestras <math>h_{m}</math>, contaminadas con ruido <math>n_m</math> en parte proveniente del propio sistema y en parte de tipo numérico, dan lugar a la señal <math>y_m</math>:

<span id="eq-3"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> y_m = y(mT_s) = \displaystyle \sum _{n=1}^{N} r_n e^{s_n mT_s } + n_m = </math> <math>  \displaystyle \sum _{n=1}^{N} r_n {z_n}^{m} + n_m </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

con <math display="inline">m=0,1,2,\ldots ,M-1</math>, y donde <math display="inline">z_n=e^{s_n T_s}</math> son los polos de <math display="inline">H(z)</math>, y a su vez <math display="inline">H(z)</math> es la representación discreta del sistema <math display="inline">H(s)</math>.

El conjunto de muestras <math>\{ y_m\} </math> alimenta al algoritmo de estimación o método de ajuste objeto de estudio en este trabajo, el cual permite estimar las frecuencias complejas naturales de resonancia <math>\hat{s}_n</math> a partir de las cuales finalmente se construye una estimación de la respuesta impulsiva que se denominará de aquí en adelante <math>\hat{h}(t)</math>.

===2.1 Método de Prony clásico===

Gaspard Riche de Prony creó el método que hoy lleva su nombre en la última parte del siglo XVIII <span id='citeF-2'></span>[[#cite-2|[2]]] para tratar de modelar el fenómeno de expansión de los gases. Prony asumió originalmente que el número <math display="inline">N</math> de frecuencias complejas naturales era conocido y la señal <math display="inline">y(t)</math> bajo estudio tenía un ruido añadido despreciable, de modo que la condición <math display="inline">h(t) \approx y(t)</math> se cumplía.

El primer paso del algoritmo de Prony consiste en organizar las muestras <math display="inline">h_m</math> obtenidas a partir del muestreo uniforme de la señal <math display="inline">y(t)</math>, con <math display="inline">m=0,1,2,\ldots , M-1</math>, en un vector columna <math display="inline">\mathbf {h}</math>:

<span id="eq-4"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> {\mathbf {h}}=\left(h_0,h_1,\ldots ,h_{M-1}\right)^{\mathit{H}} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}

donde el superíndice <math display="inline">\mathit{H}</math> denota conjugación transpuesta. Luego, tomando en cuenta que en ausencia de ruido <math display="inline">y_m=h_m=h(mT_s)</math>, el valor de cada muestra puede expandirse a partir de la ecuación <span id='eq-3'></span>[[#eq-3|(3)]] de la forma

<span id="eq-5"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \begin{array}{rlcl}h_0~=&r_1&+...+&r_N\\ h_1~=&r_1e^{s_1T_s} &+...+&r_Ne^{{(s_N)}T_s}\\ \vdots ~= & \\ h_m~=&r_1e^{s_1mT_s} &+...+&r_Ne^{{(s_N)}mT_s }\\ \vdots ~= & \\ h_{M-1}~=&r_1e^{s_1 (M-1)T_s}&+...+&r_Ne^{(s_N)(M-1)T_s } \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

o de forma más compacta

<span id="eq-6"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \begin{array}{rclcccl}h_0 & = & r_1 & + & \ldots & + & r_N\\ h_1 & = & r_1z_1 & + & \ldots & + & r_Nz_N\\ \vdots & = & \\ h_m & = & r_1z_1^m & + & \ldots & + & r_Nz_N^m\\ \vdots & = & \\ h_{M-1} & = & r_1z_1^{(M-1)} & + & \ldots & + & r_Nz_N^{(M-1)} \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|}

donde <math display="inline">\{ z_n\} </math>, con <math display="inline">n=1,2,\ldots ,N</math>, y <math display="inline">z_n=e^{s_nT_s}</math> son los polos de <math display="inline">H(z)</math>. Por otra parte, <math display="inline">H(z)</math> es la representación discreta del sistema original <math display="inline">H(s)</math>.

El sistema de ecuaciones dado por <span id='eq-6'></span>[[#eq-6|(6)]] contiene <math display="inline">2N</math> incógnitas: <math display="inline">N</math> residuos complejos desconocidos y <math display="inline">N</math> polos <math display="inline">\{ z_n\} </math> desconocidos. Por esa razón se requiere un número <math display="inline">M</math> de muestras mayor o igual a <math display="inline">2N</math> para resolver dicho conjunto de ecuaciones para <math display="inline">\{ r_n\} </math> y <math display="inline">\{ z_n\} </math>.

Luego, los valores de <math display="inline">\{ r_n\} </math> y <math display="inline">\{ z_n\} </math> se estiman tal que, mediante la técnica de mínimos cuadrados, con la <math display="inline">\hat{h}(t)</math> obtenida, se minimice el error <math display="inline">\xi </math> dado por

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\xi =\sqrt{\sum _{m=0}^{M-1}\left(\hat{h}_m-h_m\right)^2}</math>
|}
|}

Posteriormente, el método de Prony prosigue con los tres pasos siguientes:

<ol>

<li>La solución de un sistema lineal algebraico de ecuaciones de por lo menos orden <math display="inline">N</math>, para los coeficientes <math display="inline">\{ a_n\} </math> del sistema <math display="inline">H(z)</math>:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

H(z)=\frac{b_0z^K+b_1z^{K-1}+\ldots +b_{K-1}z+b_K}{a_0z^N+a_1z^{N-1}+\ldots + a_{N-1}z+a_N}=\frac{N(z)}{D(z)} </math>
|}
|}</li>

donde se asume que <math display="inline">a_0=1</math> y <math display="inline">K<N</math>. Observe que los ceros del polinomio en el denominador son los polos o frecuencias naturales de <math display="inline">H(z)</math>.
<li>El cálculo de las <math display="inline">N</math> raíces <math display="inline">\{ z_n\} </math> del polinomio característico <math display="inline">D(z)</math>:

<span id="eq-7"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

 \begin{array}{rl}

D(z)=a_0z^N+a_1z^{N-1}+\ldots +a_N = & 0\\ = a_0(z-z_N)(z-z_{N-1})\times \ldots \times (z-z_1)= & 0\\ = a_0\Pi _{n=1}^N(z-z_n)= & 0 \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}</li>
<li>La solución de un sistema de <math display="inline">N</math> ecuaciones algebraicas para los residuos <math display="inline">\{ r_n\} </math>. </li>

</ol>

Estos pasos se explican de manera detallada a continuación.

Para comenzar la estimación de coeficientes <math display="inline">\{ a_n\} </math>, el sistema dado por la ecuación <span id='eq-6'></span>[[#eq-6|(6)]] se multiplica sucesivamente por <math display="inline">(a_N,a_{N-1},\ldots ,a_1,a_0)^H</math> como se indica a continuación:

<span id="eq-8"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \begin{array}{crlll}{a_N}~\times & (h_0 & = r_1 & +\ldots+&r_N)\\ {\vdots }~~~~\times & (h_1 & = r_1z_1 & +\ldots+&r_Nz_N)\\ {a_0}~~\times & \vdots & & & {\vdots } \\ & (h_m & = r_1z_1^m &+\ldots+&r_Nz_N^m)\\ & \vdots & = & & {\vdots } \\ & (h_{M-1} & = r_1z_1^{(M-1)} &+\ldots+&r_Nz_N^{(M-1)}) \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}

entonces, también se puede realizar el proceso iniciándolo en la fila siguiente

<span id="eq-9"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \begin{array}{crlll}&(h_0 & = r_1 & +\ldots+&r_N)\\ {a_N}~\times & (h_1 & = r_1z_1 & +\ldots +&r_Nz_N)\\ {\vdots }~~~~\times & \vdots & & & {\vdots } \\ {a_0}~~\times & (h_m & = r_1z_1^m & +\ldots+&r_Nz_N^m)\\ & \vdots & = & & {\vdots }\\ & (h_{M-1} & = r_1z_1^{(M-1)} & +\ldots+&r_Nz_N^{(M-1)}) \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|}

y así sucesivamente, hasta la última fila como se observa a continuación

<span id="eq-10"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \begin{array}{ccrll}&&(h_0 & = r_1 & +...+r_N)\\ & & (h_1 & = r_1z_1 & +...+r_Nz_N)\\ & & \vdots & = & \\ {a_N} & \times & (h_m & = r_1z_1^m & +...+r_Nz_N^m)\\ {\vdots } & \times & \vdots & = & \\ {a_0} & \times & (h_{M-1} & = r_1z_1^{(M-1)} & +...+r_Nz_N^{(M-1)}) \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|}

Si nos centramos ahora en el sistema que resulta de completar la operación de multiplicación de la ecuación <span id='eq-8'></span>[[#eq-8|(8)]] y se suman los términos en cada lado de las <math display="inline">N+1</math> primeras ecuaciones, se obtiene:

<span id="eq-11"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>a_Nh_0+\ldots +a_0h_N =</math> <math>a_N(r_1 + \ldots +r_N) + \ldots +a_0(r_1z_1^N+\ldots+r_Nz_N^N) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}

la cual puede ser reordenada como

<span id="eq-12"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> a_Nh_0+\ldots +a_0h_N = </math> <math> (a_Nr_1 + \ldots +a_0 r_1 z_1^N)+</math>
|-
| style="text-align: center;" | <math>(a_Nr_2 + \ldots +a_0 r_2 z_2^N)+</math> <math>\ldots +(a_Nr_N+\ldots +a_0r_Nz_N^N) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}

Luego de notar que los términos del lado derecho de la ecuación <span id='eq-12'></span>[[#eq-12|(12)]] son todas versiones del polinomio característico <span id='eq-7'></span>[[#eq-7|(7)]], puede escribirse que

<span id="eq-13"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> a_Nh_0+\ldots +a_0h_N = 0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
|}

Luego, esta ecuación puede generalizarse mediante la repetición de un proceso similar aplicado a todos los sistemas lineales que van desde <span id='eq-8'></span>[[#eq-8|(8)]] hasta <span id='eq-10'></span>[[#eq-10|(10)]], para obtener:

<span id="eq-14"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>a_Nh_m+a_{N-1}h_{m+1}+\ldots+a_0h_{m+N}=0 </math>
|-
| style="text-align: center;" | <math> \mbox{para } \quad m=0,1,\ldots , M-N-1 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
|}

por eso, un sistema de <math display="inline">M-N</math> ecuaciones lineales con <math display="inline">N-1</math> incógnitas puede escribirse:

<span id="eq-15"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \begin{array}{rl}a_Nh_0+\ldots+a_0h_{N}= & 0\\ a_Nh_1+\ldots+a_0h_{N+1}= & 0\\ \vdots & \\ a_Nh_{M-N-1}+\ldots +a_0h_{M-1}= & 0 \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
|}

Luego, recordando que <math display="inline">a_0=1</math>, el sistema de ecuaciones <span id='eq-15'></span>[[#eq-15|(15)]] puede ser escrito nuevamente, con la forma

<span id="eq-16"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \begin{array}{rl}a_Nh_0+\ldots+a_1h_{N-1}~= &-h_{N}\\ a_Nh_1+\ldots+a_1h_{N}~= &-h_{N+1}\\ \vdots & \\ a_Nh_{M-N-1}+\ldots+a_1h_{M-2}~= &-h_{M-1} \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
|}

Y éste a su vez, puede ser reescrito en la forma matricial siguiente

<span id="eq-17"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \underbrace{\begin{pmatrix} h_0 & h_1 & \ldots & h_{N-1}\\ h_1 & h_2 & \ldots & h_{N}\\ \vdots & \vdots & \vdots & \vdots \\ h_{M-N-1} & h_{M-N-2} &\ldots & h_{M-2} \end{pmatrix}}_{\mathbf{H}} \underbrace{\begin{pmatrix} a_N\\ a_{N-1}\\ \vdots \\ a_1 \end{pmatrix}}_{\mathbf{a}}= </math> <math> \underbrace{\begin{pmatrix} -h_{N}\\ -h_{N+1}\\ \vdots \\ -h_{M-1} \end{pmatrix}}_{\mathbf{b}} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
|}

donde <math display="inline">\mathbf {H}</math> es una matriz de Hankel de orden <math display="inline">(M-N)\times N</math>, <math display="inline">\mathbf {b}</math> es el vector <math display="inline">(M-N)\times 1</math> de valores conocidos, ambos construidos a partir de las muestras <math display="inline">h_m</math> de <math display="inline">h(t)</math>, mientras que <math display="inline">\mathbf {a}</math> es el vector incógnita con los coeficientes <math display="inline">a_n</math>, con <math display="inline">n=1,2,\ldots , N</math> del polinomio característico de <math display="inline">H(z)</math>.

Como la ecuación <span id='eq-17'></span>[[#eq-17|(17)]] representa un sistema de ecuaciones linealmente independientes con más ecuaciones que incógnitas no hay manera de encontrar una solución exacta para el vector <math display="inline">\mathbf {a}</math>. La mejor estimación del vector <math display="inline">\mathbf {a}</math> se obtiene, sin embargo, mediante la técnica de aproximación por mínimos cuadrados <span id='citeF-57'></span>[[#cite-57|[57]]] de la forma:

<span id="eq-18"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \mathbf{a}=\mathbf{H}^+\mathbf{b} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
|}

donde <math display="inline">\mathbf{H}^+=({\mathbf{H}^{H}}\mathbf{H})^{-1}\mathbf{H}^H</math> es la inversa de Moore-Penrose, o pseudo-inversa, de la matriz <math display="inline">\mathbf{H}</math>.

Una vez que los valores de <math display="inline">\{ a_n\} </math> son obtenidos, las frecuencias naturales (discretas) <math display="inline">\{ z_n\} </math> del sistema <math display="inline">H(z)</math> podrían estimarse encontrando las raíces del polinomio característico dado por <span id='eq-7'></span>[[#eq-7|(7)]].

A partir de <math display="inline">\{ z_n\} </math> las frecuencias complejas naturales de resonancia <math display="inline">\{ s_n\} </math> pueden estimarse, ya que <math display="inline">z_n=e^{s_nT_s}=e^{-\sigma _nT_s}e^{\pm \jmath \omega _n T_s}</math>, y además se cumple que <math display="inline">|z_n|=e^{-\sigma _nT_s}</math> y que <math display="inline">\angle z_n=\pm \omega _nT_s</math>, entonces

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>s_n = \frac{\ln{(z_n)}}{T_s} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
|}

Finalmente, de ser necesario, el sistema de ecuaciones de la ecuación <span id='eq-6'></span>[[#eq-6|(6)]] puede ser utilizado para estimar los <math display="inline">N</math> residuos <math display="inline">\{ r_n\} </math> y construir la señal estimada <math display="inline">\hat{h}(t)</math>:

<span id="eq-20"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \underbrace{\begin{pmatrix} 1 & 1 & \ldots & 1\\ z_1 & z_2 & \ldots & z_N\\ \vdots & \vdots & \vdots & \vdots \\ z_1^m & z_2^m & \vdots & z_N^m\\ \vdots & \vdots & \vdots & \vdots \\ z_1^{M-1} & z_2^{M-1} & \ldots & z_N^{M-1}\\ \end{pmatrix}}_{\mathbf{Z}} \underbrace{\begin{pmatrix} r_1\\ r_2\\ \vdots \\ r_N \end{pmatrix}}_{\mathbf{r}}= \underbrace{\begin{pmatrix} h_0\\ h_1\\ \vdots \\ h_m\\ \vdots \\ h_{M-1} \end{pmatrix}}_{\mathbf{h}} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
|}

Y de aquí

<span id="eq-21"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \mathbf{r}=\mathbf{Z}^+\mathbf{h} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
|}

donde <math display="inline">\mathbf{Z}^+</math> es la matriz inversa de Moore-Penrose de <math display="inline">\mathbf{Z}</math>, que a su vez es una matriz de Vandermonde transpuesta.

===2.2 Método de Prony clásico utilizando filtrado con descomposición en valores singulares (DVS)===

Cuando el número exacto de frecuencias naturales complejas de <math display="inline">h_m</math> se desconoce y además el ruido <math display="inline">n_m</math> no se puede despreciar, el método clásico de Prony en general falla. Para remediar este inconveniente se puede proceder, suponiendo un valor máximo de polos <math display="inline">N_\mathrm{MAX}</math> inicial, a filtrar el ruido mediante una descomposición en valores singulares de la matriz de Hankel del sistema antes de aplicar el método de Prony clásico <span id='citeF-5'></span>[[#cite-5|[5]]]. Una versión similar a esta estrategia se conoce con el nombre de método de Kumaresan-Tufts <span id='citeF-6'></span>[[#cite-6|[6]],<span id='citeF-40'></span>[[#cite-40|40]]].

En este método, así como en el método clásico de Prony, se procede primero a la creación de una matriz de Hankel como <math display="inline">\mathbf {H}</math> en la ecuación <span id='eq-17'></span>[[#eq-17|(17)]], pero que esta vez designaremos como <math display="inline">\mathbf {Y}</math> debido a la presencia del ruido, la cual tiene la siguiente estructura

<span id="eq-22"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\mathbf {Y}}={\begin{pmatrix} y_{0} &y_{1} &...&y_{(N_\mathrm{MAX}-1)}\\ y_{1} &y_{2} &... &...\\ y_{2} &y_{3} &... &...\\ . & . &... &...\\ . & . &... &y_{M-3}\\ y_{(M-N_\mathrm{MAX}-1)} &. &... &y_{M-2}\\ \end{pmatrix}} </math>
|}
|}

de dimensiones <math display="inline">(M-N_\mathrm{MAX}) \times N_\mathrm{MAX}</math> tal que <math display="inline">N_\mathrm{MAX} > N</math>.

Luego, la matriz <math display="inline">\mathbf {Y}</math> es factorizada mediante DVS de la forma

<span id="eq-22"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \mathbf{Y}=\mathbf{U}\mathbf{\Sigma }{\mathbf{V}}^{\mathit{H}} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
|}

donde <math display="inline">\mathbf{U}</math> es una matriz unitaria de dimensiones <math display="inline">(M-N_\mathrm{MAX})\times (M-N_\mathrm{MAX})</math> cuyas columnas son los autovectores de <math display="inline">\mathbf{Y}\mathbf{Y}^{\mathit{H}}</math>, <math display="inline">\mathbf{V}</math> es una matriz de dimensiones <math display="inline">(N_\mathrm{MAX}+1) \times (N_\mathrm{MAX}+1)</math> cuyas columnas son los autovectores de <math display="inline">\mathbf{Y}^{\mathit{H}}\mathbf{Y}</math> y <math display="inline">\mathbf{\Sigma }</math> es una matriz diagonal de dimensiones <math display="inline">(M-N_\mathrm{MAX}) \times (N_\mathrm{MAX}+1)</math> la cual contiene los valores singulares de <math display="inline">\mathbf{Y}</math> (raíz cuadrada de los autovalores de <math display="inline">{\mathbf{Y}}^{\mathit{H}}{\mathbf{Y}}</math> ) ordenados en forma descendente tal que <math display="inline">\sigma _1 > \sigma _2 > ... > \sigma _{N_\mathrm{MAX}+1}</math> <span id='citeF-9'></span>[[#cite-9|[9]],<span id='citeF-58'></span>[[#cite-58|58]]].

Debido a la presencia de ruido, el número de valores singulares no nulos en <math display="inline">\mathbf{\Sigma }</math> se incrementa más allá del orden <math display="inline">N</math> del sistema, el cual no se conoce a priori. Los valores singulares originales de la matriz “limpia” <math display="inline">\mathbf{H}</math> son perturbados por el ruido en <math display="inline">\mathbf{Y}</math> y además otros valores singulares adicionales aparecen en el proceso. Dada esta situación, los errores que se deben al ruido pueden suprimirse por medio de la eliminación de los valores singulares espurios de la matriz <math display="inline">\mathbf{\Sigma } </math> <span id='citeF-6'></span>[[#cite-6|[6]],<span id='citeF-40'></span>[[#cite-40|40]]].

Con base en este principio y tomando en cuenta que los valores singulares asociados al ruido poseen menor potencia, se establece una tolerancia <math display="inline">T</math> relativa al máximo valor singular <math display="inline">\sigma _{1}</math> a partir de un umbral <math display="inline">\sigma _T</math> de acuerdo a la siguiente ecuación

<span id="eq-23"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \frac{\sigma _T}{\sigma _{1}} \approx 10^{-p} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
|}

donde <math display="inline">p</math> es el número de dígitos significativos en los datos debido a la resolución esperada, <math display="inline">T</math> es la tolerancia relativa admitida, y <math display="inline">\sigma _T</math> es el umbral de filtrado. Es natural esperar que el número de valores singulares <math display="inline">\sigma _n</math> que cumplen con la condición <math display="inline">\sigma _n > \sigma _T</math> coincida con <math display="inline">N</math> en condiciones ideales. Por otra parte, los valores singulares que cumplen con la condición <math display="inline">\sigma _n < \sigma _T</math> se consideran espurios. Un ejemplo de aplicación de este procedimiento para determinar el número verdadero de polos del sistema analizado se muestra en la <span id='img-2'></span>[[#img-2|Figura 2]], donde ya se ha calculado el valor de <math display="inline">\sigma _T</math>, umbral con el cual se determina el número de polos del sistema, que en este caso particular mostrado en la figura resulta ser <math display="inline">N=8</math>.

<div id='img-2'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 65%;"
|-
|style="padding:10px;"| [[Image:Draft_Barroso_243065238-NoisyPencil2.png|510px|Caso típico de los valores singulares de la diagonal de la matriz Σ en función de la n-ésima posición en la diagonal de dicha matriz, para una respuesta transitoria y(t) que proviene de un sistema lineal de ocho polos.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding:10px;"| '''Figura 2'''. Caso típico de los valores singulares de la diagonal de la matriz <math>\mathbf {\Sigma }</math> en función de la <math>n-</math>ésima posición en la diagonal de dicha matriz, para una respuesta transitoria <math>y(t)</math> que proviene de un sistema lineal de ocho polos
|}


Seguidamente, utilizando este criterio los valores singulares de la matriz <math display="inline">\mathbf {\Sigma }</math> con una tolerancia relativa inferior a <math display="inline">T</math> se reemplazan con ceros, para generar una matriz <math display="inline">\mathbf {\Sigma ^\prime }</math>, la cual podría interpretarse como una matriz <math display="inline">\mathbf {\Sigma }</math> de bajo nivel de ruido o “limpia” <span id='citeF-5'></span>[[#cite-5|[5]]]. Con esta nueva matriz <math display="inline">\mathbf {\Sigma ^\prime }</math> se crea luego una matriz de Hankel <math display="inline">\mathbf {Y^\prime }</math> de bajo nivel de ruido o “limpia” a partir de la expresión

<span id="eq-24"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf {Y^\prime } ={\mathbf {U}}{{\mathbf {\Sigma ^\prime }}}{{\mathbf {V}}^H} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
|}

Usando la matriz <math display="inline">\mathbf{Y^{\prime}}</math> en lugar de <math display="inline">\mathbf{H}</math> en la ecuación <span id='eq-18'></span>[[#eq-18|(18)]] se obtiene el vector <math display="inline">\mathbf{a}</math> con los coeficientes del polinomio característico de <math display="inline">H(z)</math> de la forma:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{a}=\mathbf {Y^\prime }^+\mathbf{b}</math>
|}
|}

a partir de lo cual se siguen todos los pasos subsiguientes del método clásico de Prony.

===2.3 Método de Prony usando mínimos cuadrados totales===

El uso del método de mínimos cuadrados totales en conjunto con el método de Prony, descrito en la presente sección, presupone que el número <math display="inline">N</math> de frecuencias naturales se conoce ''a priori'' y que la señal bajo análisis está contaminada de ruido. En este método se aplica una descomposición DVS a la matriz <math display="inline">{\mathbf {Y}}</math> ampliada de la forma <math display="inline">[\mathbf {Y},-\mathbf {b}]</math> <span id='citeF-7'></span>[[#cite-7|[7]],<span id='citeF-8'></span>[[#cite-8|8]]].

La modificación sustancial planteada surge cuando no se propone resolver directamente el sistema de ecuaciones <math display="inline">\mathbf {Y}\mathbf {a} = \mathbf {b}</math>, sino el sistema <math display="inline">[{\mathbf {Y}},{-\mathbf {b}}][{\mathbf {a}}^{\mathit{H}} ~~-1]^{^{\mathit{H}}}=0</math>.

Como la matriz aumentada <math display="inline">[\mathbf{Y},-\mathbf{b}]</math> contiene muestras contaminadas con ruido, se la expande mediante una descomposición en valores singulares de la forma <span id='citeF-59'></span>[[#cite-59|[59]]]:

<span id="eq-25"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>[\mathbf {Y},-\mathbf {b}] = \mathbf{U}\mathbf{\Sigma }\mathbf{V}^{\mathit{H}} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (25)
|}

Luego, comprobamos que la matriz <math display="inline">\mathbf {V}</math> de <span id='eq-25'></span>[[#eq-25|(25)]] puede representarse como

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf {V}=\begin{pmatrix}\mathbf {V_{pp}} &~\mathbf {v_{pq}} \\ \mathbf {v_{qp}} &~v_{qq} \\ \end{pmatrix} </math>
|}
|}

donde <math display="inline">\mathbf {v_{pq}}</math> es el último vector columna de <math display="inline">\mathbf {V}</math> sin su último elemento, <math display="inline">\mathbf {v_{qp}}</math> es el último vector fila de <math display="inline">\mathbf {V}</math> sin su último elemento, y <math display="inline">v_{qq}</math> es el elemento faltante de ambos vectores. Y el vector <math display="inline">-v_{qq}^{-1}{\mathbf {v_{pq}}}</math> contiene los coeficientes del polinomio característico de Prony dado en la ecuación <span id='eq-7'></span>[[#eq-7|(7)]] <span id='citeF-8'></span>[[#cite-8|[8]]]:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\mathbf {a}}\approx -v_{qq}^{-1}{\mathbf {v_{pq}}}</math>
|}
|}

A partir de aquí se aplican los pasos subsiguientes del método de Prony clásico.

===2.4 Método de haz de matrices en ausencia de ruido===

El método de haz de matrices o de haz de funciones generalizado (GPOF) es un método de un solo paso de cálculo, el cual permite estimar los polos discretos <math display="inline">\{ z_n\} </math> del sistema <math display="inline">H(z)</math> como los autovalores generalizados de un haz lineal de matrices de la forma:

<span id="eq-26"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \mathbf{H_2}-\lambda \mathbf{H_1} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (26)
|}

donde las matrices <math display="inline">\mathbf{H_1}</math> y <math display="inline">\mathbf{H_2}</math> se derivan de la matriz aumentada de Hankel <math display="inline">[\mathbf{H},-\mathbf {b}]</math>, de dimensiones <math display="inline">(M-N)\times (N+1)</math>

<span id="eq-27"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> {[\mathbf {H},-\mathbf {b}]}=\begin{pmatrix}h_{0} &h_{1}&..&h_{N}\\ h_{1} &h_{2}&..& h_{N+1}\\ . & . &..& .\\ . & . &..& .\\ .    & . &..& .\\ h_{M-N-1}& &..&h_{M-1}\\ \end{pmatrix}_{(M-N)\times ({N+1})} </math>
|}
|}

Las matrices <math display="inline">\mathbf{H_1}</math> y <math display="inline">\mathbf{H_2}</math> se construyen a partir de <math display="inline">[\mathbf{H},-\mathbf {b}]</math> eliminando de esta la segunda y la primera columna, respectivamente:

<span id="eq-27"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\mathbf {H_1}}= \begin{pmatrix}h_{0} &h_{1}& ...&h_{N-1}\\ h_{1} &h_{2}& ... &...\\ h_{2} &h_{3}& ... &...\\ . & . & ... &...\\ . & . & ... &h_{M-3}\\ h_{M-N-1}& . & ... &h_{M-2}\\ \end{pmatrix}_{(M-N) \times N} </math>
|}
|}

<span id="eq-27"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\mathbf {H_2}}= \begin{pmatrix}h_{1} &h_{2} &...&h_{N}\\ h_{2} &h_{3} &... &...\\ h_{3} &h_{4} &... &...\\ . & . &... &...\\ .    & . &... &h_{M-2}\\ h_{M-N}& . &... &h_{M-1}\\ \end{pmatrix}_{(M-N) \times N} </math>
|}
|}

Fácilmente se comprueba que las matrices <math display="inline">\mathbf{H_1}</math> y <math display="inline">\mathbf{H_2}</math> pueden ser factorizadas de la forma <span id='citeF-4'></span>[[#cite-4|[4]]]:

<span id="eq-27"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \mathbf{H_1}= \mathbf{Z_1RZ_2} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (27)
|}

<span id="eq-28"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \mathbf{H_2}= \mathbf{Z_1RZ_0Z_2} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (28)
|}

donde

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{Z_1}=\begin{pmatrix}1 & 1 & \ldots &1\\ z_1 & z_2 & \ldots & z_N\\ \vdots & \vdots & \ldots & \vdots \\ z_1^{M-N-1} & z_2^{M-N} & \ldots & z_N^{M-1} \end{pmatrix}_{(M-N)\times N} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (29)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{R}=\begin{pmatrix}r_1 & 0 & \ldots & 0\\ 0 & r_2 & \ldots & 0\\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \vdots & r_N \end{pmatrix}_{N\times N} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (30)
|}

<span id="eq-31"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{Z_0}=\begin{pmatrix}z_1 & 0 & \ldots & 0\\ 0 & z_2 & \ldots & 0\\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \vdots & z_N \end{pmatrix}_{N\times N} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (31)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{Z_2}= \begin{pmatrix}1 & z_1 & \ldots & z_1^{N-1}\\ 1 & z_2 & \ldots & z_2^{N-1}\\ \vdots & \vdots & \ldots & \vdots \\    1 & z_N & \ldots & z_N^{N-1} \end{pmatrix}_{N\times N} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
|}

Ahora bien, combinando apropiadamente las ecuaciones <span id='eq-27'></span>[[#eq-27|(27)]] y <span id='eq-28'></span>[[#eq-28|(28)]] se obtiene:

<span id="eq-33"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{H_1^+}\mathbf{H_2} = \mathbf{Z_2^+}\mathbf{Z_0}\mathbf{Z_2} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
|}

donde <math display="inline">\mathbf{H_1^+}</math> es la pseudoinversa de Moore-Penrose de <math display="inline">\mathbf{H_1}</math>:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{H_1^+}=\left({\mathbf{H_1^\mathit{H}}}\mathbf{H_1}\right)^{-1}\mathbf{H_1^\mathit{H}} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
|}

De la inspección de la ecuación <span id='eq-33'></span>[[#eq-33|(33)]], puede concluirse que los autovalores de <math display="inline">\mathbf{H_1^+}\mathbf{H_2}</math> son los autovalores de la matriz diagonal <math display="inline">\mathbf{Z_0}</math> dado que <math display="inline">\mathbf{Z_2^+}\mathbf{Z_0}\mathbf{Z_2}</math> es una forma de diagonalización de la matriz <math display="inline">\mathbf{H_1^+}\mathbf{H_2}</math>. Por tanto, los autovalores de la matriz <math display="inline">\mathbf{H_1^+}\mathbf{H_2}</math> son iguales a los elementos de la matriz <math display="inline">\mathbf{Z_0}</math> que de acuerdo a la ecuación <span id='eq-31'></span>[[#eq-31|(31)]] son las raíces <math display="inline">z_n</math> que se desea estimar.

Una vez obtenidas las raíces <math display="inline">z_n</math>, las frecuencias naturales complejas <math display="inline">s_n</math> pueden determinarse como <math display="inline">s_n =\ln (z_n)/\Delta t</math>.

===2.5 Método de haz de matrices en presencia de ruido===

La estrategia que describiremos a continuación también se conoce como método de haz de matrices con mínimos cuadrados totales <span id='citeF-58'></span>[[#cite-58|[58]]], y se utiliza cuando las muestras están contaminadas de ruido y además el número <math display="inline">N</math> de frecuencias naturales complejas del sistema bajo estudio se desconoce. Con este método los polos discretos <math display="inline">z_n</math> del sistema <math display="inline">H(z)</math> se estiman también en un solo paso como solución del problema generalizado de autovalores del haz lineal de matrices de la forma:

<span id="eq-35"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \mathbf{Y_2}-\lambda \mathbf{Y_1} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (35)
|}

donde las matrices <math display="inline">\mathbf{Y_1}</math> y <math display="inline">\mathbf{Y_2}</math> se derivan de la matriz de Hankel <math display="inline">\mathbf{Y}</math> de dimensiones <math display="inline">(M-L)\times (L+1)</math>, construida a partir de las muestras <math display="inline">y_m=h_m+n_m</math> (<span id='img-1'></span>[[#img-1|Figura 1]]):

<span id="eq-36"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\mathbf {Y}}= \begin{pmatrix}y_{0} & y_{1}&.. &y_{L}\\ y_{1} & y_{2}&.. & y_{(L+1)}\\ . & . &.. & .\\ . & . &.. & .\\ y_{M-L-2} & . &.. & y_{M-2}\\ y_{M-L-1}& &.. &y_{M-1}\\ \end{pmatrix}_{(M-L) \times ({L+1})} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (36)
|}

donde <math display="inline">L</math> se conoce como parámetro del haz, y se utiliza para eliminar algunos de los efectos del ruido presente en los datos. Los valores recomendados de <math display="inline">L</math> se ubican en el rango entre <math display="inline">M/3</math> y <math display="inline">M/2</math>. Para estos valores de <math display="inline">L</math> la varianza de los valores estimados de <math display="inline">z_n</math> debida al ruido ha resultado mínima <span id='citeF-9'></span>[[#cite-9|[9]]].

La matriz <math display="inline">\mathbf{Y_1}</math> se construye por medio de la eliminación de la última columna de <math display="inline">\mathbf{Y}</math>, mientras que la matriz <math display="inline">\mathbf{Y_2}</math> se construye mediante la eliminación de la primera columna de <math display="inline">\mathbf{Y}</math>:

<span id="eq-37"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\mathbf {Y_1}}= \begin{pmatrix}y_{0} &y_{1}& ...&y_{L-1}\\ y_{1} &y_{2}& ... &y_{L}\\ y_{2} &y_{3}& ... &y_{L+1}\\ . & . & ... &...\\ y_{M-L-2} & . & ... &y_{M-3}\\ y_{M-L-1}& . & ... &y_{M-2}\\ \end{pmatrix}_{(M-L) \times L} </math>
|}
|}

<span id="eq-37"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\mathbf {Y_2}}= \begin{pmatrix}y_{1} &y_{2} &...&y_{L}\\ y_{2} &y_{3} &... &y_{L+1}\\ y_{3} &y_{4} &... &y_{L+2}\\ . & . &... &...\\ y_{M-L-1} & . &... &y_{M-2}\\ y_{M-L}& . &... &y_{M-1}\\ \end{pmatrix}_{(M-L) \times L} </math>
|}
|}

Debido a la presencia de ruido en las muestras, <math display="inline">L-N</math> polos matemáticos falsos aparecen, acompañando las <math display="inline">N</math> frecuencias naturales complejas propias del sistema. Una descomposición DVS de <math display="inline">\mathbf {Y}</math> permite discriminar las frecuencias naturales de los polos espurios. En efecto:

<span id="eq-37"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \mathbf{Y}=\mathbf{U}\mathbf{\Sigma } \mathbf{V^\mathit{H}} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (37)
|}

donde vale la pena recordar que <math display="inline">\mathbf{\Sigma }</math> es una matriz diagonal de dimensiones <math display="inline">(M-L) \times (L+1)</math> la cual contiene los valores singulares de <math display="inline">\mathbf{Y}</math> ordenados de forma descendente tal que <math display="inline">\sigma _1 > \sigma _2 > ... > \sigma _{L+1}</math> <span id='citeF-9'></span>[[#cite-9|[9]],<span id='citeF-58'></span>[[#cite-58|58]]].

Definiendo una tolerancia <math display="inline">T</math> tal como se hizo en la Sección 2.1 (ecuación <span id='eq-23'></span>[[#eq-23|(23)]]), se construye la matriz diagonal cuadrada <math display="inline">\mathbf {\Sigma }_\mathit{T}</math> de dimensiones <math display="inline">T\times T</math> con los <math display="inline">T</math> valores singulares más significativos de la matriz <math display="inline">\mathbf{Y}</math>, y luego las matrices rectangulares truncadas <math display="inline">{\mathbf{U}_\mathit{T}}</math> y <math display="inline">{\mathbf{V}_\mathit{T}}</math>, por medio de la supresión de todas las columnas más allá de la <math display="inline">T</math>-ésima columna de <math display="inline">\mathbf{U}</math> y <math display="inline">\mathbf{V}</math>, respectivamente.

Por medio de la eliminación de la última fila de la matriz <math display="inline">\mathbf{V}_\mathit{T}</math> se crea la matriz <math display="inline">\mathbf{V_1}</math> y mediante la eliminación de la primera fila de <math display="inline">\mathbf{V}_\mathit{T}</math> se crea la matriz <math display="inline">\mathbf{V_2}</math>. Las matrices obtenidas <math display="inline">\mathbf{V_1}</math> y <math display="inline">\mathbf{V_2}</math> cumplen con la relación

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \mathbf{Y_1}=\,\mathbf {U}_\mathit{T}\mathbf {\Sigma }_\mathit{T}{\mathbf {{V_1}}}^\mathit{H}</math>
|-
| style="text-align: center;" | <math> \mathbf{Y_2}=\,\mathbf {U}_\mathit{T}\mathbf {\Sigma }_\mathit{T}{\mathbf {{V_2}}}^\mathit{H}</math>
|}
|}

Luego, tomando en cuenta que el problema de autovalores generalizado <math display="inline">(L\times T)\Rightarrow (T\times T) </math>:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{Y_2}-\lambda \mathbf{Y_1}\Rightarrow \mathbf{Y_1^+}\mathbf{Y_2}-\lambda \mathbf{I} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (38)
|}

es equivalente al problema ordinario de autovalores dado por

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{V_2^\mathit{H}}-\lambda \mathbf{V_1^\mathit{H}}\Rightarrow \mathbf{V_1^{{\mathit{H}}^+}}\mathbf{V_2}^\mathit{H}-\lambda \mathbf{I} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (39)
|}

los polos <math display="inline">z_n</math> del sistema <math display="inline">H(z)</math> se estiman como los autovalores de la matriz <math display="inline">{\mathbf{V_1^{{\mathit{H}}^+}}}{\mathbf{V_2}}^{\mathit{H}}</math> <span id='citeF-9'></span>[[#cite-9|[9]]].

==3. Resultados==

Con la intención de comparar su desempeño, a continuación se presentan los resultados obtenidos a partir de la programación en MATLAB de los métodos descritos en la Sección 2 en diferentes escenarios de simulación. Como señal de prueba se usó una compuesta de la suma de cuatro cosenos amortiguados de la forma:

<span id="eq-40"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> y(t) = r_1e^{-\sigma _1 t}\cos (\omega _1 t)+r_2e^{-\sigma _2 t}\cos (\omega _2 t)+\ldots </math>
|-
| style="text-align: center;" | <math>\dots + r_3e^{-\sigma _3 t}\cos (\omega _3 t)+ r_4e^{-\sigma _4 t}\cos (\omega _4 t)+n(t) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (40)
|}

donde tanto la potencia de ruido blanco Gausiano <math display="inline">n(t)</math> así como los valores de las frecuencias naturales complejas <math display="inline">s_n=-\sigma _n\pm \jmath \omega _n</math>, con <math display="inline">\sigma _n>0</math>, y sus residuos <math display="inline">r_n</math> asociados fueron definidos ''ad hoc'' para cada experimento.

===3.1 Comparación del desempeño de cada algoritmo en presencia de diferentes niveles de ruido añadido===

Este primer experimento fue realizado con la idea de examinar la precisión de cada uno de los métodos numéricos presentados en la Sección 2 en la estimación de las frecuencias naturales complejas asociadas al sistema lineal bajo estudio en presencia de diferentes niveles de ruido blanco Gausiano añadido. Para realizar este experimento, se seleccionó una señal <math display="inline">y(t)</math> de la forma dada por la ecuación <span id='eq-40'></span>[[#eq-40|(40)]] con <math display="inline">r_1=10</math>, <math display="inline">r_2=5</math>, <math display="inline">r_3=r_4=0</math>, <math display="inline">\sigma _1=1.1</math>, <math display="inline">\omega _1=2\pi </math>, <math display="inline">\sigma _2=1.4</math> y <math display="inline">\omega _2=4\pi </math>. Nótese que dado que <math display="inline">r_3=r_4=0</math>, cada algoritmo, de ser exitoso, solo debe estimar las frecuencias naturales complejas <math display="inline">s_1</math>, <math display="inline">s_2</math>, <math display="inline">s_5</math> y <math display="inline">s_6</math>.

En este experimento se utilizó un formato numérico de doble precisión para la estimación de los valores de <math display="inline">s_1</math>, <math display="inline">s_2</math>, <math display="inline">s_5</math> y <math display="inline">s_6</math>. La señal <math display="inline">y(t)</math> fue muestreada desde <math display="inline">t_0=0</math> hasta <math display="inline">t_f=M t_s</math>, con <math display="inline">M=</math> 80, con una frecuencia de muestreo de <math display="inline">f_s=11</math>, correspondiente a un tiempo de muestreo de <math display="inline">t_s = 1/f_s = 1/11</math>.

Para los métodos de haz de matrices con filtro DVS, y de Prony con filtro DVS, se utilizó un número máximo de componentes de la matriz diagonal <math display="inline">{\mathbf {Z}_0}</math> de <math display="inline">N_{\mathrm{MAX}}=50</math>. El factor de truncamiento <math display="inline">T</math> se fijó en <math display="inline">T=4</math>, mientras que el parámetro del haz de matrices <math display="inline">L</math> se fijó en <math display="inline">L=[5M/12]=33</math> donde la utilización de los corchetes, en este caso, simbolizan la función parte entera.

Por otra parte, para los métodos de Prony clásico y con mínimos cuadrados totales, así como para el de haz de matrices sin filtro DVS, se fijó el número de frecuencias naturales complejas a calcular en <math display="inline">N=T=4</math>.

Una vez definidos los valores de entrada, cada algoritmo fue corrido en presencia de ruido un número de <math display="inline">K=200</math> veces para 20 valores distintos de la relación señal a ruido, desde SNR=10 hasta SNR=153. Con los valores estimados de <math display="inline">\hat{s}_1</math> (<math display="inline">\hat{\sigma }_1</math> y <math display="inline">\hat{\omega }_1</math>) se estimaron, a su vez, las varianzas normalizadas <math display="inline">\mathrm{Var}(\hat{\sigma }_1) / \sigma _1^2</math> y <math display="inline">\mathrm{Var}(\hat{\omega }_1) / \omega _1^2</math> para cada valor de SNR a partir de las ecuaciones:

<span id="eq-41"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathrm{Var}(\hat{\sigma }_1)=\dfrac{\displaystyle \sum _{k=1}^{K}(\hat{\sigma }_{1,k}-\bar{\sigma }_1)^2}{K-1} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (41)
|}

<span id="eq-42"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathrm{Var}(\hat{\omega }_1)=\dfrac{\displaystyle \sum _{k=1}^{K}(\hat{\omega }_{1,k}-\bar{\omega }_1)^2}{K-1} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (42)
|}

donde <math display="inline">\hat{\sigma }_{1,k}</math> y <math display="inline">\hat{\omega }_{1,k}</math> son los valores estimados de <math display="inline">\sigma _{1}</math> y <math display="inline">\omega _{1}</math> en la <math display="inline">k</math>-ésima corrida realizada, mientras que <math display="inline">\bar{\sigma }_1</math> y <math display="inline">\bar{\omega }_1</math> son los respectivos promedios aritméticos de <math display="inline">\hat{\sigma }_1</math> y <math display="inline">\hat{\omega }_1</math>, estimados en <math display="inline">K=200</math> realizaciones.

Los resultados obtenidos se presentan en las <span id='img-3'></span>[[#img-3|Figuras 3]] y <span id='img-4'></span>[[#img-4|4]]. En tales figuras se han trazado las curvas de <math display="inline">\left(\mathrm{Var}(\hat{\sigma }_1) / \sigma _1^2\right)^{-1}</math> y <math display="inline">\left(\mathrm{Var}(\hat{\omega }_1) / \omega _1^2\right)^{-1}</math>, respectivamente, en ambos casos como una función de la relación señal a ruido (SNR), para cada uno de los métodos analizados, de manera similar a la realizada en <span id='citeF-4'></span>[[#cite-4|[4]]].

<div id='img-3'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 55%;"
|-
|style="padding:10px;"| [[Image:Draft_Barroso_243065238-ProniBasicoTLSyPencilSNRiterativounpoloresnLinea2.png|400px|Valor de [Var(̂σ₁) / σ₁²]⁻¹ en decibelios versus la relación señal a ruido (SNR) en decibelios, para cada uno de los métodos analizados.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding:10px;" | '''Figura 3'''. Valor de <math>\left[\mathrm{Var}(\hat{\sigma }_1) / \sigma _1^2\right]^{-1}</math> en decibelios versus la relación señal a ruido (SNR) en decibelios, para cada uno de los métodos estudiados
|}
<div id='img-4'></div>


{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 55%;"
|-
|style="padding:10px;"|[[Image:Draft_Barroso_243065238-ProniBasicoTLSyPencilSNRiterativounpoloimsnLinea2.png|376px|Valor de [Var(̂ω₁) / ω₁²]⁻¹ en decibelios versus la relación señal a ruido (SNR) en decibelios, para cada uno de los métodos estudiados.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding:10px;"| '''Figura 4'''. Valor de <math>\left[\mathrm{Var}(\hat{\omega }_1) / \omega _1^2\right]^{-1}</math> en decibelios versus la relación señal a ruido (SNR) en decibelios, para cada uno de los métodos estudiados
|}

===3.2 Comparación del desempeño de cada algoritmo para distintas frecuencias de muestreo===

Este segundo experimento fue realizado con la idea de examinar la precisión de cada uno de los métodos numéricos presentados en la Sección 2 cuando se estiman las frecuencias naturales complejas de un sistema en presencia de un nivel de ruido añadido fijo, al utilizar distintas frecuencias de muestreo <math display="inline">f_s</math>. Para realizar este experimento, se utilizó una señal <math display="inline">y(t)</math> tal que <math display="inline">r_1=10</math>, <math display="inline">r_2=r_3=r_4=0</math>, con <math display="inline">\sigma _1=1.1</math> y <math display="inline">\omega _1=2\pi </math>. En este caso cada algoritmo debe estimar las frecuencias complejas <math display="inline">s_1</math> y <math display="inline">s_1^\ast </math> del sistema lineal. Para realizar este experimento, se utlizó un ruido del tipo Gausiano blanco <math display="inline">n(t)</math> con una potencia relativa apropiada para proveer una relación señal a ruido resultante de SNR=60 dB.

Para todas las simulaciones realizadas en este experimento se utilizó un formato numérico de doble precisión. La señal <math display="inline">y(t)</math> fue muestreada en una ventana de observación que va desde <math display="inline">t_0=0</math> hasta <math display="inline">t_f=4</math>, para un rango de frecuencias de muestreo (<math display="inline">f_s</math>) que va desde <math display="inline">{f_s}_{MIN}=2.5</math> hasta <math display="inline">{f_s}_\mathrm{MAX}=27</math>. Para lograr esta condición hubo que cambiar para cada simulación el número de muestras <math display="inline">M</math>, tal que se cumpliese <math display="inline">M=[{f_s}(t_f-t_0)]</math> con la idea de mantener la ventana de observación constante para cada valor de <math display="inline">f_s</math> utilizado en cada simulación.

Para las simulaciones realizadas en este experimento con el método de Prony con filtro DVS y para el haz de matrices con filtro DVS, se fijó el número máximo de componentes <math display="inline">N_\mathrm{MAX}</math> de la matriz diagonal <math display="inline">{\mathbf {Z}_0}</math> en <math display="inline">N_\mathrm{MAX}=20</math>, mientras que el parámetro del haz de matrices <math display="inline">L</math> se fijó en <math display="inline">L=[5M/12]=33</math>, para un factor de truncamiento <math display="inline">T=2</math> pues son dos las frecuencias complejas naturales a calcular.

Por otra parte, para las simulaciones realizadas con el método de Prony clásico, con el método de Prony con mínimos cuadrados totales y con el método del haz de matrices sin filtro DVS, se utilizó un número de frecuencias naturales complejas <math display="inline">N=2</math>.

Luego, para cincuenta valores de frecuencia de muestreo <math display="inline">f_s</math> equiespaciadas, desde <math display="inline">f_s=2.5</math> hasta <math display="inline">f_s=27.5</math> se muestra el inverso de las varianzas normalizadas <math display="inline">\mathrm{Var}(\hat{\sigma }_1) / \sigma _1^2</math> y <math display="inline">\mathrm{Var}(\hat{\omega }_1) / \omega _1^2</math>, donde <math display="inline">\mathrm{Var}(\hat{\sigma }_1)</math> y <math display="inline">\mathrm{Var}(\hat{\omega }_1)</math> están definidas por las ecuaciones (41) y (42) para <math display="inline">K=40</math> realizaciones. Como resultado, en las <span id='img-5'></span>[[#img-5|Figuras 5]] y <span id='img-6'></span>[[#img-6|6]] se muestran las gráficas de los valores <math display="inline">\left(\mathrm{Var}(\hat{\sigma }_1) / \sigma _1^2\right)^{-1}</math> y <math display="inline">(\mathrm{Var}\left(\hat{\omega }_1) / \omega _1^2\right)^{-1}</math> en función de la frecuencia <math display="inline">f_s</math> de muestreo, para cada algoritmo bajo estudio.

<div id='img-5'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 50%;"
|-
|style="padding:10px;"| [[Image:Draft_Barroso_243065238-ProniBasicoTLSyPencilfsiterativounpoloresnLinea2.png|400px|Valor de [Var(̂σ₁) / σ₁²]⁻¹ en decibelios ''versus'' la frecuencia de muestreo fₛ, para cada uno de los métodos estudiados.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding:10px;"| '''Figura 5'''. Valor de <math>\left[\mathrm{Var}(\hat{\sigma }_1) / \sigma _1^2\right]^{-1}</math> en decibelios ''versus'' la frecuencia de muestreo <math>f_s</math>, para cada uno de los métodos estudiados
|}

<div id='img-6'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 50%;"
|-
|style="padding:10px;"| [[Image:Draft_Barroso_243065238-ProniBasicoTLSyPencilfsiterativounpoloimsnLinea2.png|400px|Valor de [Var(̂ω₁) / ω₁²]⁻¹ en decibelios ''versus'' la frecuencia de muestreo fₛ, para cada uno de los métodos estudiados.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding:10px;"| '''Figura 6'''. Valor de <math>\left[\mathrm{Var}(\hat{\omega }_1) / \omega _1^2\right]^{-1}</math> en decibelios ''versus'' la frecuencia de muestreo <math>f_s</math>, para cada uno de los métodos estudiados
|}

===3.3 Comparación de los tiempos de cómputo que cada algoritmo necesita para determinar las frecuencias naturales complejas de un sistema dado===

Este experimento fue realizado con el objetivo de comparar el tiempo de cómputo de cada algoritmo para estimar las ocho frecuencias naturales complejas que forman parte de la señal <math display="inline">y(t)</math> dada por la ecuación <span id='eq-40'></span>[[#eq-40|(40)]] variando la ventana temporal de observación de la señal y manteniendo constante la tasa de muestreo.

En este experimento la señal <math display="inline">y(t)</math> fue construida artificialmente usando los siguientes valores: <math display="inline">r_1=10</math>, <math display="inline">r_2=7</math>, <math display="inline">r_3=3</math>, <math display="inline">r_4=1</math>, <math display="inline">\sigma _1=1.1</math>, <math display="inline">\sigma _2=1.4</math>, <math display="inline">\sigma _3=2</math>, <math display="inline">\sigma _4=3</math>, <math display="inline">\omega _1=2\pi </math>, <math display="inline">\omega _2=4\pi </math>, <math display="inline">\omega _3=6\pi </math> y <math display="inline">\omega _4=10\pi </math> usando un formato numérico de doble precisión para todos los cómputos.

Las muestras de <math display="inline">y(t)</math> fueron tomadas durante cuatro ventanas temporales de duración distinta, desde <math display="inline">t_0=0</math> hasta <math display="inline">t_f=M t_s</math>, con <math display="inline">M=\{ 100,200,300,400\} </math>, usando en todos los casos la misma frecuencia de muestreo de <math display="inline">f_s=11</math>, correspondiente a un tiempo de muestreo de <math display="inline">t_s = 1/f_s = 1/11</math>.

Para las simulaciones realizadas en este experimento con el haz de matrices con filtro DVS y con el método de Prony con filtro DVS, se utilizó un número máximo de componentes de la matriz diagonal <math display="inline">{\mathbf {Z}_0}</math> de <math display="inline">N_\mathrm{MAX}=20</math>. El parámetro del haz de matrices se fijó para cada caso en <math display="inline">L=[5M/12]</math> y el factor de truncamiento se fijó en <math display="inline">T=8</math>.

Por otra parte, para las estimaciones realizadas con los métodos de Prony clásico, de Prony con mínimos cuadrados totales y de haz de matrices sin filtro DVS, se fijó un número de frecuencias naturales complejas a estimar de <math display="inline">N=T=8</math>.

En la <span id='table-1'></span>[[#table-1|Tabla 1]] se muestra el tiempo de cómputo en milisegundos empleado por cada algoritmo durante la estimación de la frecuencias naturales complejas <math display="inline">s_n</math> según el número <math display="inline">M</math> de muestras tomadas a <math display="inline">y(t)</math>.

<div class="center" style="font-size: 75%;">'''Tabla 1'''. Tiempo de cómputo en milisegundos empleado por cada algoritmo durante la estimación de la frecuencias <math>s_n</math> según el número <math>M</math> de muestras tomadas a <math>y(t)</math></div>

<div id='tab-1'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%;width:auto;" 
|-style="text-align:center"
!   style="text-align:left"|  Algoritmo !!  <math display="inline">M=100</math>   !!   <math display="inline">M=200</math>  !! <math display="inline"> M=300</math> !! <math display="inline"> M=400</math> 
|- style="text-align:center"
|  style="text-align:left"| Prony clásico      
|  0.00029 
|  0.00033 
|  0.00039 
|  0.00045
|- style="text-align:center"
|  style="text-align:left"| Prony con DVS 
|  0.00165 
|  0.00218 
|  0.00490 
|  0.00427 
|- style="text-align:center"
|  style="text-align:left"| Prony con  mínimos cuadrados  totales
|  0.00167 
|  0.00434 
|  0.01150 
|  0.02474
|- style="text-align:center"
| style="text-align:left"| Haz de matrices  
|  0.00616 
|  0.04337 
|  0.13920 
|  0.27570
|- style="text-align:center"
| style="text-align:left"| Haz de matrices con DVS
|  0.00460
|  0.02660
|  0.07272 
|  0.17300
|}

==4. Análisis de resultados==

Del análisis comparativo de los resultados observados de todas las simulaciones realizadas durante el desarrollo de este trabajo, puede concluirse que, en presencia de bajos niveles de ruido, todos los métodos estudiados son capaces de estimar las frecuencias naturales de resonancia del sistema bajo estudio. Sin embargo, del análisis comparativo de los resultados expuestos en la Sección 3.1, en las <span id='img-3'></span>[[#img-3|Figuras 3]] y  <span id='img-4'></span>[[#img-4|4]], puede concluirse que el método convencional más preciso y robusto cuando la señal <math display="inline">y(t)</math> está contaminada con ruido es el método de haz de matrices convencional con descomposición DVS. Este resultado es consistente con los estudios reportados en <span id='citeF-4'></span>[[#cite-4|[4]]], en los cuales se observa que la precisión de los resultados provenientes de este método está muy cerca del límite de Cramer-Rao, el cual es el límite teórico de máxima precisión posible para este tipo de algoritmos estudiados.

Del estudio realizado en la misma Sección 3.1 llama la atención el buen desempeño que tiene el método de Prony cuando utiliza como insumo una matriz <math display="inline">\mathbf {Y}</math> que se ha limpiado mediante una descomposición DVS. También es importante resaltar la significativa diferencia de tolerancia al ruido que tiene el método de haz de Matrices con y sin filtrado DVS, pues al no utilizar la descomposición DVS, el método de haz de Matrices pierde capacidad de tolerar al ruido presente en las muestras.

De los resultados mostrados en la Sección 3.2, particularmente en las <span id='img-5'></span>[[#img-5|Figuras 5]] y  <span id='img-6'></span>[[#img-6|6]], pudo constatarse que en la medida en que la frecuencia de muestreo aumenta dentro de una misma ventana de observación, la precisión que brinda cada una de las versiones de los métodos de Prony y de haz de matrices aumenta, con la excepción del método de Prony clásico que más bien parece mejorar su desempeño en la medida en que la frecuencia de muestreo <math display="inline">f_s</math> se aproxima a <math display="inline">6\omega _1/(2\pi )</math> (<span id='img-6'></span>[[#img-6|Figura 6]]). De hecho, de acuerdo a la información presentada en estas gráficas, el método de Prony clásico puede brindar una mejor precisión que varios de los métodos utilizados si es pequeña la cantidad de muestras que se dispone dentro de una ventana de observación.

Otro aspecto interesante que puede observarse comparando los resultados de las <span id='img-3'></span>[[#img-3|Figuras 3]] y <span id='img-4'></span>[[#img-4|4]], y también comparando los resultados de las <span id='img-5'></span>[[#img-5|Figuras 5]] y <span id='img-6'></span>[[#img-6|6]], es que, en general, la precisión con la que los métodos analizados calculan el factor de amortiguamiento <math display="inline">\sigma _1</math> es menor que la con la que calculan la frecuencia angular de resonancia <math display="inline">\omega _1</math>. También llama la atención el desempeño del método de Prony clásico cuando, en presencia de poco ruido añadido, casi iguala la precisión del método de haz de matrices con SVD para el cálculo de <math display="inline">\omega _1</math>, mientras que el mismo método de Prony clásico es el menos preciso para el cálculo de <math display="inline">\sigma _1</math> bajo las mismas condiciones. La mayoría de estas observaciones son consistentes con las realizadas en el estudio comparativo presentado en <span id='citeF-54'></span>[[#cite-54|[54]]].

Del estudio realizado en la Sección 3.3, pudo constatarse a través de los resultados expuestos en la <span id='table-1'></span>[[#table-1|Tabla 1]], que las distintas versiones basadas en el método de Prony son más rápidas en el cálculo que los métodos que utilizan la teoría del haz de Matrices. Esta diferencia en la velocidad de cálculo se incrementa en la medida en que el número de muestras utilizadas se hace más grande. De acuerdo al análisis realizado al código del programa que se elaboró, este fenómeno se debe a que las rutinas de MATLAB que calculan las raíces del polinomio característico de Prony son más rápidas que las rutinas que calculan los autovalores de las matrices dadas para el haz de matrices. Estos resultados son congruentes con el estudio presentado en <span id='citeF-38'></span>[[#cite-38|[38]]], en el cual se menciona la dificultad de utilizar el método del haz de matrices en la medida en que se incrementa el número de muestras <math display="inline">M</math> tomadas a la señal <math display="inline">y(t)</math>.

==5. Conclusiones y recomendaciones==

La respuesta impulsiva de un sistema lineal tiene como rasgo distintivo un conjunto de frecuencias naturales complejas de resonancia las cuales le son propias. Durante el desarrollo de este trabajo, se estudiaron y probaron en distintos escenarios de simulación, algunos de los más importantes algoritmos de cálculo numérico que existen para estimar tales frecuencias complejas naturales del sistema. Entre estos algoritmos, se probó el método de Prony clásico, el método de Prony con filtrado usando descomposición en valores singulares, el método de Prony usando mínimos cuadrados totales, el método de haz de matrices sin filtro y el método de haz de matrices con filtrado mediante descomposición en valores singulares. También durante el desarrollo de este trabajo, se realizó una revisión bibliográfica que evidenció la importancia que tienen estos métodos en diferentes aplicaciones de la física y de la ingeniería.

Se evaluó y comparó mediante simulación la precisión obtenida en la estimación de las frecuencias naturales complejas <math display="inline">s_n</math> de cada uno de los algoritmos en presencia de ruido. Pudo concluirse que el método de haz de matrices con descomposición y filtrado DVS es el más robusto en presencia ruido de los métodos convencionales estudiados, tal y como sugiere la bibliografía revisada. De hecho, es difícil encontrar un algoritmo que supere este desempeño pues su precisión está muy cerca del límite teórico de precisión de Cramer-Rao <span id='citeF-4'></span>[[#cite-4|[4]]] para algoritmos de este tipo. Por tanto, el método del haz de matrices puede utilizarse en aplicaciones en que el tiempo de cálculo no sea una variable crítica a tomar en cuenta, debido a su excelente precisión en ambientes de medición ruidosos.

El notable desempeño en presencia de ruido del método del haz de matrices con descomposición DVS se debe, justamente, a la descomposición y filtrado DVS, la cual permite discriminar y eliminar los valores singulares asociados al ruido presente en la matriz <math display="inline">\mathbf {Y}</math>. De hecho, el método pierde su robustez cuando no se utiliza la descomposición DVS. En este mismo sentido, cuando se utiliza filtrado a partir de una descomposición DVS en el método de Prony clásico en presencia de ruido se obtiene un desempeño comparable al del método haz de matrices.

La frecuencia de muestreo y el número de muestras utilizadas para crear la matriz <math display="inline">\mathbf {Y}</math> son variables importantes cuando se quiere aumentar la precisión de la estimación de cada frecuencia natural compleja con ruido presente. A excepción del método de Prony clásico, los métodos estudiados en los escenarios dados aumentan su precisión cuando la frecuencia de muestreo aumenta.

Por otra parte, el método clásico de Prony y sus variantes presentaron mayor velocidad de cálculo que los métodos basados en un haz de matrices, debido a que las rutinas utilizadas para el cálculo de las raíces del polinomio característico de Prony son más rápidas que las que estiman los autovalores de las matrices definidas para el método de haz de matrices. Esto podría indicar que si la relación señal a ruido fuese lo suficientemente alta y si se tiene una aplicación que necesite realizar el cálculo de <math display="inline">s_n</math> tan rápido como sea posible, el método de Prony y sus variantes pudiesen ser más apropiados que el método del Haz de Matrices y sus variantes.

==Agradecimientos==

Proyecto financiado por el Concurso de Proyectos Regulares de Investigación, año 2019, código LPR19-03, Universidad Tecnológica Metropolitana, Santiago de Chile.

==Referencias==

<div class="auto" style="text-align: left;width: auto; margin-left: auto; margin-right: auto;font-size: 85%;">

<div id="cite-1"></div>
[[#citeF-1|[1]]] Lee W., Sarkar T.K., Moon H., Salazar-Palma M. Computation of the natural poles of an object in the frequency domain using the Cauchy method. IEEE Antennas and Wireless Propagation Letters, 11:1137-1140, 2012.

<div id="cite-2"></div>
[[#citeF-2|[2]]] Prony R. Essai éxperimental et analytique: Sur les lois de la dilatabilité de fluides élastique et sur celles de la force expansive de la vapeur de l'alkool à différentes températures. Journal de l'École Polytechnique, 1:24-76, 1795.

<div id="cite-3"></div>
[[#citeF-3|[3]]] Easley D. Prony's method: Determining the number of exponential modes and the optimal sample period. MSc. Thesis, Oklahoma State University, United States of America, 1982.

<div id="cite-4"></div>
[[#citeF-4|[4]]] Hua Y., Sarkar T.K. Generalized pencil-of-function method for extracting poles of an EM system from its transient response. IEEE Transactions on Antennas and Propagation, 37(2):229-234, 1989.

<div id="cite-5"></div>
[[#citeF-5|[5]]] Younan N.H., Taylor C.D. On using the SVD-Prony method to extract poles of an EM system from its transient response. Electromagnetics, 11(2):223-233, 1991.

<div id="cite-6"></div>
[[#citeF-6|[6]]] Porat B., Friedlander B. On the accuracy of the Kumaresan-Tufts method for estimating complex damped exponentials. IEEE Transactions on Acoustics, Speech, and Signal Processing, 35(2):231-235, 1987.

<div id="cite-7"></div>
[[#citeF-7|[7]]] Rahman M.D., Yu K.B. Total least square approach for frequency estimation using linear prediction. IEEE Transactions on Acoustics Speech and Signal Processing, 35(10):1440-1454, 1987.

<div id="cite-8"></div>
[[#citeF-8|[8]]] Gavin H.P. Total least squares. Duke University, Fall 2017. http://people.duke.edu/~hpgavin/SystemID/CourseNotes/TotalLeastSquares.pdf.

<div id="cite-9"></div>
[[#citeF-9|[9]]] Sarkar T.K., Pereira O. Using the matrix pencil method to estimate the parameters of a sum of complex exponentials. IEEE Antennas and Propagation Magazine, 37(1):48-55, 1995.

<div id="cite-10"></div>
[[#citeF-10|[10]]] Crow M.L., Singh A. The matrix pencil for power system modal extraction. IEEE Transaction on Power Systems, 20(1):501-502, 2005.

<div id="cite-11"></div>
[[#citeF-11|[11]]] Xu M.M., Xiao L.Y., Wang H.F. A Prony-based method of locating short-circuit fault in DC distribution system. 2nd IET Renewable Power Generation Conference (RPG 2013), 1-4, September 2013.

<div id="cite-12"></div>
[[#citeF-12|[12]]] Yang L., Jiao, Z., Kang, X., Wang, X. A novel matrix pencil method for real-time power frequency phasor estimation under power system transients. 12th IET International Conference on Developments in Power System Protection (DPSP 2014), 12-32, 2014.

<div id="cite-13"></div>
[[#citeF-13|[13]]] Tawfik M.M., Morcos M.M. On the use of Prony method to locate faults in loop systems by utilizing modal parameters of fault current. IEEE Transactions on Power Delivery, 20(1):532-534, 2005.

<div id="cite-14"></div>
[[#citeF-14|[14]]] Costa F.F., Fernandes D.A., De Almeida L.A.L., Naidu S.R. Prony's method versus FFT for analyzing power converters signals. 2005 European Conference on Power Electronics and Applications, 11-14 September 2005.

<div id="cite-15"></div>
[[#citeF-15|[15]]] Chen J., Fu P., Méndez-Barrios C., Niculescu S., Zhang, H. Stability and instability intervals of polynomially dependent systems: An matrix pencil analysis. 2017 American Control Conference (ACC),  24-26 May, 2017.

<div id="cite-16"></div>
[[#citeF-16|[16]]] Tesche F. On the analysis of scattering and antenna problems using the singularity expansion technique. IEEE Transactions on Antennas and Propagation, 21(1):53-62, 1973.

<div id="cite-17"></div>
[[#citeF-17|[17]]] Hua Y., Sarkar T.K. A discussion of E-pulse method and Prony's method for radar target resonance retrieval from scattered field. IEEE Transactions on Antennas and Propagation, 37(7):944-946, 1989.

<div id="cite-18"></div>
[[#citeF-18|[18]]] Sabio V. Target recognition in ultra-wideband sar imagery. Army Research Laboratory TAD-A283 462, 1994.

<div id="cite-19"></div>
[[#citeF-19|[19]]] Andreev M.V., Drobakhin O.O. Feature of Prony's method application for natural frequencies estimation from the frequency response. 2016 8th International Conference on Ultrawideband and Ultrashort Impulse Signals (UWBUSIS), 5-11 September, 2016.

<div id="cite-20"></div>
[[#citeF-20|[20]]] Yoshizawa A., Uchida S. A discriminant-based RMSE improvement technique for classical Prony method in small array radars. 2020 17th European Radar Conference (EuRAD), pp. 298-301, 10-15 January 2021.

<div id="cite-21"></div>
[[#citeF-21|[21]]] Drissi K.E.K., Poljak D. The matrix pencil method applied to smart monitoring and radar. 17th International Conference on Computational Methods and Experimental Measurements, 13-24, May 2015. 

<div id="cite-22"></div>
[[#citeF-22|[22]]] Harmer S.W., Andrews D.A., Rezgui N.D., Bowring N.J. Detection of handguns by their complex natural resonant frequencies. IET Microwaves, Antennas Propagation, 4(9):1182-1990, 2010.

<div id="cite-23"></div>
[[#citeF-23|[23]]] Harmer S.W., Cole S.D., Bowring N. J., Rezgui N. D., Andrews D. On body concealed weapon detection using a phased antenna array. Progress In Electromagnetics Research, 124:187-210, 2012.

<div id="cite-24"></div>
[[#citeF-24|[24]]] Li M., Li S., Yu Y., Ni X., Chen R. Design of random and sparse metalens with matrix pencil method. Optics Express, 26:24702-24711, 2018.

<div id="cite-25"></div>
[[#citeF-25|[25]]] Kokkinos T., Adve R.S., Sarris C.D. Spectral analysis of negative refractive index metamaterials utilizing signal processing techniques and time-domain simulations. IEEE/ACES International Conference on Wireless Communications and Applied Computational Electromagnetics, pp. 409-412, 3-7 April, 2005.

<div id="cite-26"></div>
[[#citeF-26|[26]]] Singh A. Analysis of leaky wave antennas using the matrix pencil method. MSc. Thesis. Conconcordia University (Canada), 2014.

<div id="cite-27"></div>
[[#citeF-27|[27]]] Sabett K.F., Katehi L.P.B., Sarabandi K. Wavelet-based CAD modeling of microstrip discontinuities using least square Prony's method. 1997 IEEE MTT-S International Microwave Symposium Digest, 1799-1802, 8-13 June, 1997.

<div id="cite-28"></div>
[[#citeF-28|[28]]] Avramenko D.V., Andrejev V.G. Spectral estimation of the photoplethysmographic signal by the double sided Prony method. 2019 8th Mediterranean Conference on Embedded Computing (MECO), 1-4, 10-14 June, 2019.

<div id="cite-29"></div>
[[#citeF-29|[29]]] Bhuiyan M., Malyarenko E.V., Pantea M.A., Seviaryn F.M., Maev R.G. Advantages and limitations of using matrix pencil method for the modal analysis of medical percussion signals. IEEE Transactions on Biomedical Engineering, 60(2):417-426, 2013.

<div id="cite-30"></div>
[[#citeF-30|[30]]] Hossein Q., Hojjat A. A comparative study of signal processing methods for structural health monitoring. Journal of Vibroengineering, 18(4):2186-2204, 2016.

<div id="cite-31"></div>
[[#citeF-31|[31]]] Fernández A., de Santiago L., López Guillén E., Rodríguez Ascariz J. M., Miguel Jiménez J. M., Boquete L. Coding Prony’s method in MATLAB and applying it to biomedical signal filtering. BMC Bioinformatics, 19(451): 2186-2204, 2018.

<div id="cite-32"></div>
[[#citeF-32|[32]]] Pushkareva A. V., Markuleva M. V. The research of noise stability of the Prony's method while processing of cardiological time series. 2019 International Multi-Conference on Engineering, Computer and Information Sciences (SIBIRCON), 0476–0479, 2019.

<div id="cite-33"></div>
[[#citeF-33|[33]]] Fricke S. N., Seymour J.D., Battistel M. D., Freedberg D. I., Eads C. D., Augustine M. P. Data processing in NMR relaxometry using the matrix pencil. Journal of Magnetic Resonance, 313: 106704, 2020.

<div id="cite-34"></div>
[[#citeF-34|[34]]] Bauman G., Bieri O. Matrix pencil decomposition of time-resolved proton MRI for robust and improved assessment of pulmonary ventilation and perfusion. Magnetic resonance in medicine, 77(1): 336-342, 2016.

<div id="cite-35"></div>
[[#citeF-35|[35]]] Wang Y., Lee H., Apte D. Quantitative NMR spectroscopy by matrix pencil methods. International Journal of Imaging Systems and Technology, 4(3): 201-206, 1992.

<div id="cite-36"></div>
[[#citeF-36|[36]]] Ehler M., Kunis S., Peter T., Richter C. A randomized multivariate matrix pencil method for superresolution microscopy. Electronic Transactions on Numerical Analysis, 51:63-74, 2019.

<div id="cite-37"></div>
[[#citeF-37|[37]]] Li Q., Xiao X., Kikkawa T. Noninvasive blood glucose level detection based on matrix pencil method and artificial neural network. Journal of Electrical Engineering & Technology, 16(4): 2183-2190, 2021.

<div id="cite-38"></div>
[[#citeF-38|[38]]] Laroche J. The use of the matrix pencil method for the spectrum analysis of musical signals. The Journal of the Acoustical Society of America, 94: 1958-1965, 1993.

<div id="cite-39"></div>
[[#citeF-39|[39]]] Avramenko D. V., Andrejev V. G. Spectral analysis of light reflections from cosmic objects by the modified Prony's method. 2018 7th Mediterranean Conference on Embedded Computing (MECO), 1-4, 2018.

<div id="cite-40"></div>
[[#citeF-40|[40]]] Berti E., Cardoso V., González J. A., Sperhake U. Mining information from binary black hole mergers: A comparison of estimation methods for complex exponentials in noise. Physical Review D, 75(12): 124017, 2007.

<div id="cite-41"></div>
[[#citeF-41|[41]]] Adve R. Target identification using matrix pencil, Final report. University of Toronto (Canada), 2017.

<div id="cite-42"></div>
[[#citeF-42|[42]]] Chaparro-Arce D., Gallego A., Albarracin-Vargas F., Gutierrez, C., Vega F., Pedraza C. Matrix pencil method applied to the compression of audio data in naval operations. 2020 IEEE International Conference on Computational Electromagnetics (ICCEM), 254-256, 2020.

<div id="cite-43"></div>
[[#citeF-43|[43]]] Chaparro-Arce D., Gutierrez S., Gallego A., Gutierrez C., Vega F., Pedraza C. Extraction of complex natural resonances of ships acoustic signals using matrix pencil method. 2020 IEEE Colombian Conference on Communications and Computing (COLCOM), 1-4, 2020.

<div id="cite-44"></div>
[[#citeF-44|[44]]] Persichkin A., Shpilevoy A. About the method of estimating the parameters of seismic signals. IKBFU’s Vestnik Theoretical and Experimental Physics, 10: 122-125, 2015.

<div id="cite-45"></div>
[[#citeF-45|[45]]] Li M., Henry M. Signal processing methods for coriolis mass flow metering in two-phase flow conditions. 2016 IEEE International Conference on Industrial Technology (ICIT), 690-695, 2016.

<div id="cite-46"></div>
[[#citeF-46|[46]]] Zaky Y. Y., Fortino N., Dauvignac J. -Y., Seyfert F., Olivi M., Baratchart L. Comparison of SEM methods for poles estimation from scattered field by canonical. 2020 IEEE Radar Conference (RadarConf20), 1-6, 2020.

<div id="cite-47"></div>
[[#citeF-47|[47]]] Hidalgo P. Métodos de Prony y Pisarenko para análisis espectral. BSc. Tesis, Escuela Politécnica Nacional (Ecuador), 1985.

<div id="cite-48"></div>
[[#citeF-48|[48]]] Castillo R., Ortiz J., Calleros G. Análisis de estabilidad de los eventos ocurridos en Laguna Verde. Memorias del XVI Congreso Anual de la SNM / XXIII Reunión Anual de la SMSR, 1-13, 2005.

<div id="cite-49"></div>
[[#citeF-49|[49]]] Issouribehere P. E., Barbero J. C., Issouribehere F., Rodríguez J. Análisis de oscilaciones subsincrónicas derivadas de fallas en sistemas de 500 KV. Experiencias de aplicación del Método de Prony. XIV Encuentro Regional Iberoamericano de la Cigré, 2011.

<div id="cite-50"></div>
[[#citeF-50|[50]]] Issouribehere P. E., Barbero J.C., Issouribehere F. Desarrollo de una herramienta computacional para la detección de modos de oscilación en sistemas de potencia basada en el análisis de Prony. XV Encuentro Regional Iberoamericano de la Cigré, 2013.

<div id="cite-51"></div>
[[#citeF-51|[51]]] Quiroz C. Aplicaciones del modelo autoregresivo y del algoritmo de Prony. MSc. Tesis, Instituto Politécnico Nacional (Mexico), 2002.

<div id="cite-52"></div>
[[#citeF-52|[52]]] Arizaga-Jasso A. Análisis de predicción de tendencia en acciones financieras empleando el método de Prony. MSc. Thesis, Instituto Tecnológico y de Estudios Superiores de Occidente (Mexico), 2020.

<div id="cite-53"></div>
[[#citeF-53|[53]]] Almunif A., Fan L., Miao Z. A tutorial on data-driven eigenvalue identification: Prony analysis, matrix pencil, and eigensystem realization algorithm. International Transactions on Electrical Energy Systems, 30(4): 1-17, 2020.

<div id="cite-54"></div>
[[#citeF-54|[54]]] Sarrazin F. Sharaiha A., Pouliguen P., Chauveau J., Collardey S., Potier P. Comparison between matrix pencil and Prony methods applied on noisy antenna responses. 2011 Loughborough Antennas and Propagation Conference, 1-4, 2011.

<div id="cite-55"></div>
[[#citeF-55|[55]]] Sarrazin F., Chauveau J., Pouliguen P., Potier P., Sharaiha A. Accuracy of singularity expansion method in time and frequency domains to characterize antennas in presence of noise. IEEE Transactions on Antennas and Propagation, 62(3): 1261-1269, 2014.

<div id="cite-56"></div>
[[#citeF-56|[56]]] Lee J. H., Cho Y. S., Jeong J. H. Performance analysis of the LS Prony method for estimating parameters of damped sinusoids. IET Radar, Sonar and Navigation, 13(11): 1918–1933, 2019.

<div id="cite-57"></div>
[[#citeF-57|[57]]] Leon S. J. Linear Algebra with its applications. 9th Edition. Pearson, 2015.

<div id="cite-58"></div>
[[#citeF-58|[58]]] Rao S. M., Vechinski D. A. Finite Conducting Bodies: TDIE Solution. Rao S. M. (Ed.), Time Domain Electromagnetics. Academic Press, 1999.

<div id="cite-59"></div>
[[#citeF-59|[59]]] Brunton S. L., Kutz J. N. Data-driven science and engineering: Machine learning, dynamical systems and control. Cambridge University Press, 2019.