spectrum.m

arma模型的功率谱估计代码。还包括系统的极零点分析等完整软件包。

function [w,P] = spectrum(a, b, N)
% 计算传输函数H(z) = B/A的频谱, 在[0,pi)之间取N个点
alength = length(a);
blength = length(b);
f = 0: 1/N: 1;
f = f(1:N);
w = f*pi;
for n=1:N
    ar = 0;
    ma = 0;
    for m=1:alength
        ar = ar + a(m)*exp(-j*w(n)*(m-1));
    end
    for m=1:blength
        ma = ma + b(m)*exp(-j*w(n)*(m-1));
    end
    P(n) = (abs(ma/ar))^2;
    %P(n) = (abs(ma)^2)/(abs(ar)^2);
end
%test
% clear;
% close all;
% a = [1, -1.3817, 1.5632, -0.8843, 0.4096];
% b = [1, 0.3544, 0.3508, 0.1736, 0.2401];
% [f, p] = spectrum(a, b, 128);
% plot(f, p);