Prony谱分析研究
这是最近研究的一个很有意思的东西~
之前发了一篇笔记式的AR谱分析算法的东西，当时说了，AR模型是谱分析里面最最最基础的模型算法，它的思想就是通过一个采样得到的信号x[n]，算出一系列AR参数ai来，使得下式重构出来的信号与原来的信号误差最小：
\(x_{re}[n]=-\sum_{i=1}^{p}a_ix[n-i]\)
或者说计算出一系列AR参数ai来，使得\(x[n]=-\sum_{i=1}^{p}a_ix[n-i]\)成立！
其中ai的个数，也就是上式中的p，就是AR模型的阶数，之前那篇博文就列举和对比了几种计算ai的算法及仿真。
然后就是今天的主角，Prony模型，如果要阅读下去，首先要对上面的那几行要有所理解，因为Prony是基于AR模型来求解的！！
首先什么是Prony模型呢？Prony模型就是说，我们把原来的信号分解成下式：
\(x(t)=\sum A_ie^{(\alpha_i+j\omega_i)t}\)
说明几点，首先就是说，上面的式子表示，Prony模型构建思想就是把原始信号x(t)分解成一系列的不同幅度，不同衰减系数，不同频率的信号的叠加。就好像傅里叶变换是把原始信号分解成一系列不同频率幅度的正弦信号的和一样，只是我们这里多了个衰减系数αi，其他参数也解释一下，ωi就是每个子信号的频率，Ai是信号的幅度，必须一提的是，Ai可以是个复数，什么意思？就是说每个子信号还可以有一定的初始相位，懂？
哦，对，另外一点就是，离散信号的话，就是把上式的t变成nΔt就行了嘛。
然后Prony模型要干的事自然就是把这个振幅相位频率衰减全部给求解出来啦~没错，Prony模型就是根据采样得到的信号x[n]把每个子信号上述三个参数给解出来！！
但是怎么解呢？首先为了方便讨论，我们先把上面的分解的式子做点处理，只是为了书写方便而已，就是定义一个b和z：（注：下面公式里面把初始相位θ抽取出来了，所以这里的A就是实数了！）【决定了，以后该是乖乖用Latex敲公式再截图了！！唯一比较烦的就是还要去背景的白色。。】
接下来推导之前，先回到我们的AR模型：
AR模型解出来的a下脚标是从1开始到p的，所以我们可以额外定义：
\(a_0=1\)
那么就可以得到下面的式子成立：
然后再将我们的Prony模型代进去，就可以得到下面几步的推导成立(别怕，都很简单的推导来的。。)：
为什么要这么变换？因为我们需要先把z给解出来，所以把b移到一边，然后有为了让式子对一切n成立，所以把n也移到一边了，所以最后我们得到了上面最后的式子如果想对一切n成立，那么必然有第二个求和符号Σ里面的东西恒为0，也就是：
所以，我们只要先算出原本信号x[n]的AR模型参数ai，i=1,2,3…p，p为模型阶数，然后再定义a0=1，解上面的多项式方程，解出来的p个根就是我们的zi了，再看看我们上面是怎么定义zi的？我们只要知道我们信号的采样频率，就可以知道采样周期Δt，那么就可以把p个衰减系数αi和频率ωi都给解出来了！！
然后，我们再看看怎么把幅度解出来，呵呵，zi已知的情况下，不就是解下面这个方程了么？
我们解出所有的参数后就可以算出原本信号的功率谱了：
额，不知道我说清楚了没？
下面依旧老规矩，还是仿真：
首先我们先构建一个信号，它由三个衰减正弦信号组成：
第一个信号幅度是4，衰减系数是4，频率是4；
第二个幅度是2，衰减系数是3，频率是31.5；
第三个幅度是3，衰减系数是1，频率15；
mathematica code
然后写个算Prony参数的函数，不过这之前需要一个计算AR模型参数的函数，按照我之前那篇的研究结果，我比较喜欢Marple算法：
mathematica code
【这里容我扯个题外话，之前那篇AR模型博文的Marple里面，对于C矩阵，我是采用Table建表，然后还用了Σ求和函数，但是随着模型阶数p一增大，Mathematica的计算速度慢到无法忍受，后面换成现在这个矩阵点乘来代替，提高了80+倍的运算速度~然后后面的代码我就不说了，我都是做过一些优化的。】
然后就是计算Prony模型参数了，下面这个函数最后返回四个量，a, z, A, ValidLength，分别是Marple算法下AR模型的a参数，包含衰减系数和频率信息的z，包含幅度信息和初始相位信息的A，还有一个ValidLength；这个我要解释一下ValidLength，就是说，我们上面的仿真里面，使用了3个衰减正弦波来构成一个信号，但是现实中，我们并不可能事前获知信号里面有多少个衰减正弦分量，而根据上面的推导，模型阶数p决定了算出来的正弦衰减分量的个数，所以我们把p取大一些，像上面代码所示，p取的是信号长度的0.4倍，然后我们找出算出来的p个分量里面，|A|最大的ValidLength个分量，就是这个意思。
mathematica code
然后我们计算出了各个衰减正弦分量，再用这些分量叠加，得到我们模型重构的信号，可以发现基本完美重构了这个原始信号：
mathematica code
虽然你可能觉得这不是废话么？上面的推导不就是说明这么件事情的么？但是，亲，其实不是的，我们的模型阶数远远大于原本用来生成这个信号的衰减正弦波的数量(本仿真中是3)，那么也就是说，算出来的分量里面，会有很多别的分量混在里面，最可怕的是什么呢？别的分量里面有衰减系数α是正的分量，懂什么意思么？衰减系数为正，正激励震荡，也就是说随着n的增大，这个分量的值会越来越大(而衰减的正好相反，会逐步变小)以至于到不可忽略的地步，所以如果计算AR模型的算法不好的话，这里重构可能不会和原来的波形重叠的！！
我们可以写个函数，以非常直观的方式看看我们计算出来的各个分量，然后每个分量按|A|的大小排序，再用特殊颜色重点标出前ValidLength个|A|最大的分量，函数如下：
mathematica code
然后运行结果如下，由于p很大【数据长度100的0.4倍，也就是40】，所以给大家看前面一部分就是了：
仔细看看上表，可以发现前6个的|A|远远大于后面的，而且正好表征着上面的三个正弦衰减信号，对比一下，可以发现都是对的，|A|的值有一点点误差，但是基本可以接受。
然后就是后面，可以发现我之前说的那种“正”衰减出现了，但是由于它的|A|极小，负15次方数量级了，在n从1取到100这个范围内，基本可以无视，但是如果信号长度超长，达到十万，那么你可以想象一下这些正的分量的印象有多可怕，就算幅度很小，也会造成不可估量的误差！！！
然后我们再来看看功率谱，写个函数算功率谱，改成分贝单位Db，然后标出最大的三个频率的位置：
mathematica code
有人可能要说了，假设我现在知道信号是由3个正弦衰减信号构成的，那么我们取p=6会怎么样？【之所以要取两倍，看上面那个表就知道，我们原本信号写的是exp的叠加，而不是sin或者cos，所以要两个正弦，频率一正一负合成一个exp】
下面是结果，首先是重构的信号：
看到么？这个时候就不会重叠了，原因就是我们用了一个6阶的AR模型来描述这个信号，而6阶根本描述不来，所以才会出现这种状况！！
回头看看我们原本说的AR模型的思想，x[n]=-∑aix[n-i]，也就是说，如果6阶的AR模型足够精确的话，那么我们只要知道x[1]到x[6]，以及6阶AR模型的6个a，这么12个数我们就可以重构出这个信号来了！！但事实是，不行，所以就会出现上面的情况了。【从这个角度思考，可以考虑用这个做数据压缩，虽然基本肯定是有损的！！】
然后再看看各个Prony模型参数：
频率还算“勉强”算的不错，但是衰减系数和幅度什么的就完全不行了，原因就是上面所说的。我们也可以认为这个算法，阶数对频率不是很敏感【还是说频率对阶数不敏感？？傻傻分不清楚了。。】。
然后我们再看看功率谱:
其实可以预料到的，因为频率基本估算正确了，但是幅度就肯定是错的了。所以如果你想用谱特征来识别一些东西的话，然后你的识别算法又对幅度很敏感，那么模型阶数一定要选高一些！！
【完】
本文内容遵从CC版权协议，转载请注明出自http://www.kylen314.com
你好，我没有Mathematica的基础，请问能贴出MATLAB的代码吗？多谢了~
先前，我也尝试过用Python写过Prony’s method，苦于数学基础较为薄弱，一般信号处理专著上的推导公式涉及到SVD-TLS之类的数学，没能吃透。
额，不好意思，Matlab代码暂无，但是个人觉得并不难实现，只要获得了AR模型的参数，那么本文前面的推导编程是很容易实现的（其实就是解两个方程）。python如果你用的是numpy或者Scipy这个的话，那么写起来难度应该和Matlab差不多。。不太清楚你们需要这个的需求，如果只是为了实现算法，而不太追求计算时间的话，其实可以直接暴力解决大型线性方程，而不用去实现教科书那些递归解法。
后面求解拟合参数那部分看得懂，就是前面求解AR模型的Marple LS算法没太看明白（特别是优化后的代码），请问哪里能找到这个算法的推导呢，我找了半天都没找到。
如果要推导的话需要找Marple的那篇论文，教科书上好像没怎么见过。。如果是算法流程的话随便找篇中文期刊论文就可以了。。或者你看之前那篇AR的博文http://www.kylen314.com/archives/1530中的代码【In[513]那张图】，虽然是Mma写的，但是可读性还是很好地；简单解释一下，代码里面Table是构造矩阵的意思，CP是二维矩阵，b是一维矩阵，Inverse是求逆，上角标T是倒置；