cumprod([var 1-rc(2:L+1).^2] ); end

function [rc,n_obs,resf] = burg_su(x,L)

%BURG_SU Burg type AR estimator for multiple segments of unequal length.
%   [RC,N_OBS,RES] = BURG_SU(SIG,MAX_ORDER) estimates a single vector of AR-
%   reflectioncoefficients RC up to order MAX_ORDER, using data from 
%   multiple segments simultaneously. The segments may be of unequal
%   length.
%
%   x can be a cell array, where each element contains a single signal
%   in a column, or a matrix of several segments of equal length.
%
%   Example:
%   - 2 signals x and y of 5 observations and a signal z of
%     100 observations can be denoted as:
%     x = {[x(1)  y(1)    [z(1)
%           x(2)  y(2)     z(2)
%           ...   ...      ...
%           x(5)  y(5)]    ...
%                          ...
%                          z(100)]
%
%   For segments of equal length, the algorithm in BURG_S is more efficient.
%
%   See also: BURG_S, BURG, DATA_SEGMENTS.

if ~iscell(x),
    if size(x,1) == 1, x = x'; end
    x = {x};
end   
ns = length(x);

f = x; b = x;
rc(1) = 1;
for j = 1:L,
    num = 0; den = 0;
    v1 = cell(ns,1); v2 = cell(ns,1);
    for i = 1:ns,
        v1{i} = b{i}(1:size(b{i},1)-1,:); v2{i} = f{i}(2:end,:);
        num = num+sum(dot(v1{i},v2{i}));
        den = den+sum(dot(v1{i},v1{i}))+sum(dot(v2{i},v2{i}));
    end  %for i = 1:ns,
    k = -2*num/den;
    rc = [rc k];
    for i = 1:ns,
        f{i} = v2{i}+k*v1{i}; b{i} = v1{i}+k*v2{i};
    end %for i = 1:ns,
end %for j = 1:L,

switch nargout
case 2
    n_obs = [];
    for i = 1:ns
        sz = size(x{i});
        n_obs = [n_obs sz(1)*ones(1,sz(2))];
    end
case 3
    var = 0;
    n_obs = [];
    for i = 1:ns,
        sz = size(x{i});
        n_obs = [n_obs sz(1)*ones(1,sz(2))];
        var = var+sum(sum(x{i}.^2));
    end
    var = var/sum(n_obs);
    resf = cumprod([var 1-rc(2:L+1).^2] );
end