leastsquarep2.m

基于普通最小二乘法的ARMA模型谐波频率估计

% -----------------------------------------------------------------
%------------------王程  2008-12-10--------------------------------
%------------------仿真作业一：p=2---------------------------------
%-----------------------------------------------------------------
q = 1 ;
num = 1000 ;        %重复计算次数
T = 128 ;
M = 40 ;
p = 2 ;
n = 1 : T ;
A = zeros(40, 2*p) ;
b = zeros(40, 1) ;
y = ones(2*p+1, 1);
lumda = 1 ;

while q < num                       %重复1k次

x = sqrt(20) * sin(2*pi*0.2*n) + sqrt(2) * sin(2*pi*0.213*n) + randn(1, T);  % generate signals
r = xcorr(x, 'biased');             %有偏的自相关函数
for m = 1 : M
    A(m, : ) = r(T+m+2*p-1 : -1 : T+m) ;
    b(m) = -r(T+2*p+m) ;
end
a = inv(A' * A) * A' * b ;          %用最小二乘法求解方程组
y(1) = a(4) ; 
y(2) = a(3) ; 
y(3) = a(2) ; 
y(4) = a(1) ; 

z = roots(y) ; 
for i = 1 : 2 * p
    z(i) = 1 / z(i) ;               %求出极点
end

fi1 = atan( abs( imag(z(1)) / real(z(1)) ) ) / (2*pi) ;        %求解谐波频率
fi2 = atan( abs( imag(z(3)) / real(z(3)) ) ) / (2*pi) ;


if fi1 ~= 0 && fi2 ~= 0
    if fi1 < fi2
        f1(lumda) = fi1 ;         %f1是小频率
        f2(lumda) = fi2 ;
    else
        f1(lumda) = fi2 ;
        f2(lumda) = fi1 ;
    end
    lumda = lumda + 1 ;
end
q = q + 1 ;
end
err = 0;                          %误差平方和
for i = 1 : lumda-1
    err = err + (f1(i)-0.2)^2 + (f2(i)-0.213)^2 ;
end
err = err / lumda ;

figure, hist([f1',f2'],  100) ;