用Burg法估計AR模型的引數。

藉助如圖所示的格型預測誤差濾波器，伯格法通過求出前向預測誤差和後向預測誤差的平均功率來選取最佳的反射係數k，使誤差的平均功率取得最小值，進而通過反饋求出模型係數和噪聲方差，該演算法可直接通過分析觀測資料得到需要的模型引數，而不用取求解計算量較大自相關函式，

burg法的實現原理

其實AR模型用來估計功率譜是很好用的，比傳統的Welch演算法好很多，矩形性很好，用來估計訊號的矩形等效頻寬效果很好

用Burg法估計AR模型的引數的流程圖

matlab程式碼為：CSDN中沒有matlab專屬的，為了好看就這樣貼程式碼了，哼~~~~~~~~~~~~~

function [psdviaBurg, f, p] = Burg(x, Fs, varargin)
%MYBURG      根據burg演算法實現的AR模型功率譜計算
% psdviaBurg 根據burg演算法求出的功率譜值
% f          頻率軸引數
% p          模型階次
% x          輸出訊號
% Fs         取樣率
% varargin   若為數值型，則為AR模型階次
%            若為字串，則為定階準則，AR模型階次由程式確定
%
% 解析輸入引數內容
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, 20e5, Fs);
psdviaBurg = E(end)*abs(h).^2./Fs;
psdviaBurg=psdviaBurg/abs(max(psdviaBurg));
psdviaBurg=(10*log10(abs(psdviaBurg)));
 

還缺個函式，就是根據給定的資料來計算AR模型的引數和誤差，並最小化誤差，這個包含了AR模型的核心的演算法

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 = [a; 0] + k(m)*[0; conj(flipud(a))];
    
    % 當前階次的誤差功率
    E(m + 1) = (1 - conj(k(m))*k(m))*E(m);
end