ARMA、MA谱估计算法——MATLAB源码

以下是ARMA、MA谱估计算法的MATLAB源码：

function [PSD_ARMA,freq] = ARMA_psd(data,p,q,nfft,fs)
% data: 数据
% p: AR 阶数
% q: MA 阶数
% nfft: FFT 点数
% fs: 采样率
% PSD_ARMA: PSD 估计值
% freq: 频率向量

if nargin < 5
    error('Not enough input arguments!')
end

if ~isvector(data)
    error('Data must be a vector!')
end

if ~isscalar(p) || p < 0 || mod(p,1) ~= 0
    error('AR order must be a non-negative integer!')
end

if ~isscalar(q) || q < 0 || mod(q,1) ~= 0
    error('MA order must be a non-negative integer!')
end

if ~isscalar(nfft) || nfft < length(data)
    error('FFT length must be greater than or equal to the length of data!')
end

if ~isscalar(fs) || fs <= 0
    error('Sampling frequency must be a positive number!')
end

% 计算自相关系数
[acf,lags] = xcorr(data,p+q,'biased');
acf = acf(lags >= 0);

% 构造 Yule-Walker 方程
R = toeplitz(acf(1:p));
r = acf(2:p+1);

% 解 Yule-Walker 方程得到 AR 系数
ar = -inv(R)*r;

% 构造 MA 模型
ma = [1;zeros(q-1,1);-ar(p:-1:1)];

% 计算 PSD
[H,freq] = freqz(1,ma,nfft,fs);
PSD_ARMA = abs(H).^2;

% 画图
plot(freq,10*log10(PSD_ARMA))
xlabel('Frequency (Hz)')
ylabel('PSD (dB/Hz)')
title('ARMA PSD Estimate')
grid on

end

注意，这里的ARMA模型是指自回归滑动平均模型，不是自回归移动平均模型（ARMA）模型。如果需要ARMA模型的PSD估计，可以使用MATLAB自带的`pwelch`函数。