利用burg算法，用matlab写一个估计函数功率谱密度的程序

以下是用MATLAB实现Burg算法的程序，用于估计信号功率谱密度：

function [Pxx, w] = burg_psd(x, order, Nfft)
% x: 输入信号
% order: Burg算法的阶数
% Nfft: FFT的点数

% 求自相关系数
[r, lags] = xcorr(x, 'biased');
r = r(lags >= 0);
lags = lags(lags >= 0);

% 初始化Burg算法参数
a = zeros(order + 1, order + 1);
e = r(1);
for i = 1:order
    % 计算反射系数
    k = -a(i, 1:i)' * r(i+1:-1:2) / e;
    % 更新预测误差
    e = e * (1 - k' * k);
    % 更新AR系数
    a(i+1, 1:i+1) = [a(i, 1:i) + k(end:-1:1)', k(end:-1:1)];
end

% 计算功率谱密度
w = linspace(0, pi, Nfft);
Pxx = abs(polyval(a(end, :), exp(-1i*w))).^2 / e;

使用示例：

% 生成一个信号
fs = 1000;
t = 0:1/fs:1-1/fs;
x = cos(2*pi*50*t) + randn(size(t));

% 计算功率谱密度
[Pxx, w] = burg_psd(x, 10, 1024);

% 绘制功率谱密度图
plot(w/pi*fs/2, 10*log10(Pxx))
xlabel('Frequency (Hz)')
ylabel('Power/frequency (dB/Hz)')
title('Burg PSD Estimate')

其中，x为输入信号，order为Burg算法的阶数，Nfft为FFT的点数。函数的输出为功率谱密度Pxx和对应的频率向量w。