Yule-Walker AR法的功率谱估计：pyulear函数

（1）Pxx=pyulear（x，order）：用Yule-Walker AR法对离散时间信号x进行功率谱估计。输入参数order为 AR模型的阶数。如果x为实信号，则返回结果为“单边”功率谱；如果x为复信号，则返回结果为“双边”功率谱。

（2）Pxx=pyulear（x，order，NFFT）：参数NFFT用来指定FFT运算所采用的点数

如果x为实信号、NFFT为偶数，则Pxx的长度为（NFFT/2+1）；

如果x为实信号、NFFT为奇数，则Pxx的长度为（NFFT+l）/2；

如果x为复信号，则Pxx的长度为NFFT；

NFFT的默认值为256。

（3）［Pxx，w］=pyulear（…）：输出参数w为和估计PSD的位置一一对应的归一化角频率答派，单位为rad/sample，其范围如下：

如果x为实信猜埋号，则w的范围为［0，pi］；

如果x为复信号，则w的范围为［0，2*pi］。

（4）［Pxx，f］=pyulear（…，Fs）：同时返回和估计PSD的位置一一对应的线性频率f，单位为Hz，参数Fs为采样频率，单位也是Hz。当Fs为空矩阵“［］”时，则使用默认值1 Hz。输出参数f的范围如下：

如果x为实信号，则f的清兆贺范围为［0，Fs/2］；

如果x为复信号，则f的范围为［0，Fs］。

（5）［Pxx，w］=pyulear（…，‘twosided’）：在［0，2*pi］区间上进行功率谱的“双边”估计；twosided也可以由onesided代替。

（6）［Pxx，f］=pyulear（…，Fs，‘twosided’）：在［0，Fs］区间上进行功率谱的“双边”估计；twosided也可以由onesided代替。

（7）pyulear（…）：没有输出参数，在当前图形窗口里绘制出PSD估计结果图，坐标分别为dB和归一化频率。

［例4-4］用Yule-Walker AR 法进行PSD 估计，结果如图4-14 所示。

Fs=500；t=0：1/Fs：1；

x=sin（2*pi*60*t）+4*sin（2*pi*110*t）+sin（2*pi*210*t）+randn（size（t））；

pyulear（x，20，［］，Fs）。

图4-14 Yule-Walker AR法功率谱估计结果图

1、将下面函数复制保存为myBurg.m文件。就可以计算。预测误差E 、系数a 、误差功率psdviaBurg

function [psdviaBurg, f, p,a,E] = myBurg(x, Fs, varargin)

%MYBURG 根据burg算法实现的AR模型功率谱计算

% psdviaBurg 根据burg算法求出的功率谱值

% f 频率轴参数

% p 模型阶次

% a AR模型参数

%E AR模型误差

% x 输出信号

% Fs 采样率

% varargin 若为数值型，则为AR模型阶次

%若为字符串，则为定阶准则，AR模型阶次由程序确定

%

% $Author: lskyp

% $Date: 2010.6.26

% 解析输入参数内容

error(nargchk(3, 3, nargin))% 该函数的输入必须为三个个

if strcmp(class(varargin{1}), 'double')

p = varargin{1}

elseif ischar(varargin{1})

criterion = varargin{1}

else

error('参数2必须为数值型或者字符串')

end

x = x(:)

N = length(x)

% 模型参数求解

if exist('p', 'var') % p变量是否存在，存在则不需要定阶，直接使用p阶

[a, E] = computeARpara(x, p)

else % p不存在，需要定阶，定阶准则即criterion

p = ceil(N/3)% 阶次一如凯春般不超过信号长度的1/3

% 计算1到p阶的误差

[a, E] = computeARpara(x, p)

% 根据误差求解目标函数最小值

kc = 1:p + 1

switch criterion

case 'FPE'

goalF = E.*(N + (kc + 1))./(N - (kc + 1))

case 'AIC'

goalF = N.*log(E) + 2.*kc

end

[minF, p] = min(goalF)% p就是目标函数最小的位置，也即定阶准则给出的阶次

% 使用p阶重新求解AR模型参数

[a, E] = computeARpara(x, p)

end

% 计算功率谱密度

[h, f] = freqz(1, a, [], Fs)

psdviaBurg = E(end)*abs(h).^2./Fs

function [a, E] = computeARpara(x, p)

% 根据信号序列x和阶次p计算AR模型参数和误差

N = length(x)

% 初始值

ef = x% 前向预测误差

eb = x% 后向预测误差

a = 1% 初始模型渣耐参数

E = x'*x/N% 初始误差

k = zeros(1, p)% 为反射系数预分配空间，提高循环速度

E = [E k]% 为误差预分配空间，提高速度

for m = 1:p

% 根据burg算法步骤，首先计算m阶的反射系数

efm = ef(2:end)% 前一阶次的前向预测误差

ebm = eb(1:end - 1)% 前一阶次的后向预测误差

num = -2.*ebm'*efm % 反射系数的分子项

den = efm'*efm + ebm'*ebm% 反射系数的分母项

k(m) = num./den% 当前阶次的反射系数

% 更新前后向预测误差

ef = efm + k(m)*ebm

eb = ebm + conj(k(m))*efm

% 更新模型系数a

a = [a0] + k(m)*[0conj(flipud(a))]

% 当前阶次的孙并误差功率

E(m + 1) = (1 - conj(k(m))*k(m))*E(m)

end

2、例如：

参考论坛中的例子。

randn('state', 1)

x = randn(100, 1)

y = filter(1, [1 1/2 1/3 1/4 1/5], x)

pburg(y, 4, [], 1000)

[psd, f, p,a,E] = myBurg(y, 1000, 'FPE')

figure

a,E

plot(f, 10*log10(psd))

结果：图还是一样，就补贴出来了。

a =

1.0000

0.3116

0.3647

0.2086

0.2088

0.0425

E =

1.02240.98590.91670.89890.86440.8629

---