generate_ar.m

surrogates,替代数据生成算法源码

% Generate an AR-based surrogate
% AR-order is determined using the minimum description length (MDL)
% plot (MDL(:,1),MDL(:,2))     yields the MDL plot
%
% H0: signal is generated by a linear Gaussian stochastic process
% set random seed using randn('seed',s)
%
% Usage: [Xs,MDL] = generate_AR (X, max_ord)

function [Xs,MDL] = generate_AR (X, max_ord)

if (nargin<1)
	Xs = [];
	return;
end
if (nargin<2)
	max_ord = 20;
end

m1 = mean(X);
s1 = std(X);
X = (X-m1)./s1;

M = 2:max([2 max_ord]);
n_it = length(M);
pp = length(X);

E = zeros (n_it,1);
for n=1:n_it
	[AA{n} E(n)] = aryule (X, M(n));
end

MDL = [M(:) pp.*log10(E)+M(:).*log10(pp)];

[a b] = min(MDL(:,2));
m_opt = M(b);
A = -AA{b};
A(1) = [];

Xs1 = randn(1,pp+m_opt);
Xs1(1:m_opt) = randn(m_opt,1);
for p=m_opt+1:pp
	Xs1(p) = Xs1(p) + sum(A.*Xs1(p-1:-1:p-m_opt));
end

Xs1 = Xs1((m_opt+1):end);
Xs = X;
Xs(:) = Xs1(:);
Xs = Xs./std(Xs)*s1+m1;