前阵子自学了一下谱分析的东西，嘛~也就一点皮毛，现在把笔记放上来，Mathematica版本的。。
其实说白了AR就是谱分析里面一个弱到不能再弱的算法了，但是没办法，后面很多新方法，都是以此为基础的，而且。。俺们老板92年发表的一个谱分析的算法也是以此为基础的。。囧。。
谱分析，就是根据一段信号有限时间范围内的采样得到的样本，这个样本既包含噪声，也包含信号，谱分析目的就是为了把原始信号的功率谱给计算出来！
哦，对了，千万千万别问我问什么不用福利叶变换。。。。。
我们这里实验用的是一个衰减的信号，采样频率50Hz，50个点，模型阶数为6，其实这些参数怎么样都好。
(*清除变量*)
Clear["Global`*"];
Subscript[f, s] = 50;(*采样频率*)
M = 50;(*点数*)
tRange = Range[0, (M - 1)/Subscript[f, s], 1/Subscript[f, s]];
X = N[Table[
E^(-4 t) Cos[8 Pi t] + E^(-3 t) Cos[6 Pi t], {t, tRange}]];
ListLinePlot[X, PlotRange -> All]
p = 6;(*AR模型阶数*)
AR模型的基础就是把对原始信号x[n]建立了一下模型：
\(x[n]=-\sum_{i=1}^{p}x[n-i]\alpha_i+w[n]\)
其中w[n]是噪声，而\(\alpha\)参数就是我们整个模型所需要求解的参数，参数确定，模型就确定了；
p是阶数，看上面的mathematica代码也知道。。
然后这个模型的求解有很多种，这里仿真4中，其中前三种可以从数学上证明等效，也就是结果是一样的，计算出来。
第一种是最小二乘法，简单易懂，就是我们采样信号为x，根据AR模型重构出来(所谓的递推)的信号是x_bar，那么就是要误差平方和最小，推导我不管了，太乱了，这里直接贴mathematica代码，如果懂代码的，几眼就能看懂是什么意思。
(*最小二乘*)
yy[i_] := 0;
yy[i_ /; 0 < i <= M] := X[[i]];
y = Transpose[{Table[yy[i], {i, 1, p + M}]}];
Y = Table[yy[i - j], {i, 1, p + M}, {j, 1, p}];
Subscript[a, ML] =
Flatten[-(Transpose[y].Y).Inverse[Transpose[Transpose[Y].Y]]];
然后就是第二种方法，这种方法的有名之处在于，他可以衍生出一种递归解法。。
这种算法本身计算量比较大，因为要计算很多自相关系数。。。欧，对了，这个r函数在下面的那个算法中有定义。。。就是计算自相关系数rx[n]的
(*Yule-Walker解法*)
Subscript[r, 1] = Table[r[n], {n, 1, p}];
Subscript[R, 1] = Table[r[i - j], {i, 1, p}, {j, 1, p}];
Subscript[a, YW] = LinearSolve[Subscript[R, 1], -Subscript[r, 1]];
下一种就是上面那个算法的递归解法，计算速度和内存消耗会小很多，一般随便一本谱分析的书都会有下面这个递推的推导，由于书写推导涉及到大量的矩阵操作，作为不方便写数学公式的博客，就只能请各位感兴趣的话，自己去看书吧~
(*Levinson-Durbin递归模型*)
(*自相关系数*)
Subscript[r, x][x_, n_] := (1/Length[x])*
Total[(#1*RotateLeft[#1, n] & )[
Join[Table[0, {Length[x]}], x]]];
a[1, 1] := k[1];
r[n_] := Subscript[r, x][X, n];
k[1] := -(r[1]/r[0]);
\[Sigma][1] := (1 - k[1]^2)*r[0];
DD[n_] := r[n + 1] + Sum[a[n, i]*r[n + 1 - i], {i, 1, n}];
a[n_, n_] := a[n, n] = k[n];
k[n_] := -(DD[n - 1]/\[Sigma][n - 1]);
\[Sigma][n_] := \[Sigma][n] = (1 - k[n]^2)*\[Sigma][n - 1];
a[m_, i_] := a[m, i] = a[m - 1, i] + k[m]*a[m - 1, m - i];
Subscript[a, LD] = Table[a[p, i], {i, 1, p}];
然后之前说过，上面三种算法算出来的\(\alpha\)都是一样的，数学上可以证明。。
下面这个Marple算法是一种叫做Burg算法的改进，他的计算复杂度很大，但是，谱分析结果很好，就是大概这么一回事儿~后面有证据。。
(*Marple算法*)
x[n_] := X[[n + 1]];
c[i_, j_] :=
Sum[x[n - i]*x[n - j], {n, p, M - 1}] +
Sum[x[n - p + i]*x[n - p + j],
{n, p, M - 1}];
CP = Table, {i, 1, p}, {j, 1, p}];
b = Table, {i, 1, p}];
Subscript[a, Marple] = -b . Transpose[Inverse[CP]];
Grid[Prepend[
Table[{Subscript[\[Alpha], i], Subscript[a, ML][[i]],
Subscript[a, LD][[i]],
Subscript[a, YW][[i]], Subscript[a, Marple][[i]]}, {i, 1,
p}],
{\[Alpha], "最小二乘法", "L-D递归解法", "L-W解法", "Marple算法"}],
Background -> {None, {Lighter[Yellow, 0.9], {White,
Lighter[Blend[{Blue, Green}], 0.8]}}},
Dividers -> {{Darker[Gray, 0.8], Darker[Gray, 0.8],
Darker[Gray, 0.8], Darker[Gray, 0.6],
Darker[Gray, 0.6]}, {Darker[Gray, 0.6],
Darker[Gray, 0.6], {}, Darker[Gray, 0.6]}},
Alignment -> {{Center, Center, Center}},
Frame -> Darker[Gray, 0.6], ItemStyle -> 14,
Spacings -> {Automatic, 0.8}]
既然Marple算出来的参数和前面三种相差那么大，那么重构出原来的信号的效果如何呢？看下面的重构图，可以明显发现，前面三个算法，搞了这么多种，还是被Marple爆了十条街。。。大概。。。
Subscript[x, ReconML][n_] := Piecewise[{{X[[n]], n < = p},
{-Sum[X[[n - i]]*Subscript[a, ML][[i]], {i, 1, p}], n > p}}];
Subscript[x, ReconLD][n_] := Piecewise[{{X[[n]], n < = p},
{-Sum[X[[n - i]]*Subscript[a, LD][[i]], {i, 1, p}], n > p}}];
Subscript[x, ReconYW][n_] := Piecewise[{{X[[n]], n < = p},
{-Sum[X[[n - i]]*Subscript[a, YW][[i]], {i, 1, p}], n > p}}];
Subscript[x, ReconMarple][n_] := Piecewise[{{X[[n]], n < = p},
{-Sum[X[[n - i]]*Subscript[a, Marple][[i]], {i, 1, p}],
n > p}}];
(ListLinePlot[{X, #1, X - #1},
PlotLegends -> {"原始信号", "重构信号", "误差曲线"},
PlotStyle -> {Red, Dashed, Black}, PlotRange -> All,
PlotLabel -> "最小二乘法"] & )[
Table[Subscript[x, ReconML][i], {i, 1, M}]]
(ListLinePlot[{X, #1, X - #1},
PlotLegends -> {"原始信号", "重构信号", "误差曲线"},
PlotStyle -> {Red, Dashed, Black}, PlotRange -> All,
PlotLabel -> "Marple"] & )[
Table[Subscript[x, ReconMarple][i], {i, 1, M}]]
最后再来简单的看一下普分析出来的结果。。。嘛~结果表明，虽然有两种频率在里面，但是由于采样时间比较短的缘故，所以最后只检测出一种频率，一个检测出4Hz那个，一个检测出3Hz那个，嘛~这种初等的算法就是这样的啦~
(*谱估计*)
Function[{f0, x, xmin, xmax},
Plot[f0, {x, xmin, xmax}, PlotLabel -> "Marple算法",
Epilog ->
Function[
Bias, (({PointSize[Large], Point[#1], Arrow[{#1, #1 + Bias}],
Text[Style[#1, Bold, Red],
Flatten[#1 + Bias], {-1, 0}]} & ) /@
{{x /. Last[#1], First[#1]}} & )[
Maximize[{f0, x >= xmin, x < = xmax}, {x}]]][
{Subscript[f, s]/20, -0.01}], PlotRange -> All,
PlotLabel -> "功率谱"]][
\[Sigma][p]/
Abs[1 + Sum[
Subscript[a, Marple][[i]]*
Exp[-((I*\[Omega]*2*Pi*i)/Subscript[f, s])], {i, 1, p}]]^
2, \[Omega], 0, Subscript[f, s]/2]
Function[{f0, x, xmin, xmax},
Plot[f0, {x, xmin, xmax}, PlotLabel -> "最小二乘算法",
Epilog ->
Function[
Bias, (({PointSize[Large], Point[#1], Arrow[{#1, #1 + Bias}],
Text[Style[#1, Bold, Red],
Flatten[#1 + Bias], {-1, 0}]} & ) /@
{{x /. Last[#1], First[#1]}} & )[
Maximize[{f0, x >= xmin, x < = xmax}, {x}]]][
{Subscript[f, s]/20, -0.01}], PlotRange -> All,
PlotLabel -> "功率谱"]][
\[Sigma][p]/
Abs[1 + Sum[
Subscript[a, ML][[i]]*
Exp[-((I*\[Omega]*2*Pi*i)/Subscript[f, s])], {i, 1, p}]]^2,
\[Omega], 0, Subscript[f, s]/2]
【完】
本文内容遵从CC版权协议，转载请注明出自http://www.kylen314.com
