function [Un,len_filter] = PhaSpa2VoltCoef(xn,p)
% 構造 Volterra 自適應 FIR 濾波器的輸入信號矢量 Un
% [Un, len_filter] = PhaSpa2VoltCoef(xn, p)
% 輸入參數：
% xn           相空間中的點序列(每一列為一個點)
% p            Volterra 級數階數
% 輸出參數：
% Un           Volterra 自適應 FIR 濾波器的輸入信號矢量 Un

[m,xn_cols] = size(xn);         % m 為重構維數，xn_cols 為訓練樣本個數
Un(1,:) = ones(1,xn_cols);      % FIR 濾波器的抽頭輸入信號矢量 Un (第一個系數為 1)
h_cols_1 = 0;

for k = 1:p

    clear h;
    h(1,:) = zeros(1,k);        % k 階 Volterra 核的下標 (第一個為 0,0,0... )
    i = 1;
    
    % 當上一行最后一列數值達到 m-1 結束循環
    while h(i,end)<m-1
        i = i + 1;
        % 從后往前考察上一行每一列
        for j = k:-1:1
            % 當上一行第 j 列數值達到 m-1 時，這一行的第 1 至第 j+1 列的數值均為上一行第 j+1 列數值加 1，其余不變
            if (h(i-1,j)==m-1)
                h(i,1:j+1) = ones(1,j+1) * (h(i-1,j+1)+1);
                h(i,j+2:end) = h(i-1,j+2:end);
                break;
            end
            
            % 當上一行數值都沒有達到 m-1 時，這一行第 1 列數值加 1，其余不變
            h(i,1) = h(i-1,1)+1;
            h(i,2:end) = h(i-1,2:end);            
        end
    end
    %disp('------------------')
    h;                          % 構造 k 階 Volterra 核的下標
    h_cols_1 = [h_cols_1;h(:,1)];
    
    index = m - h;
    [index_rows,index_cols] = size(index);
    
    un = zeros(index_rows,xn_cols);
    
    %for q = 1:xn_cols
    %    vector = xn(:,q);
    %    array = vector(index);    % 從列向量中提取出x(n),x(n-tau),x(n-2*tau)...
    %    un(:,q) = prod(array,2);    % 計算x(n),x(n-tau),x(n-2*tau)...之間的乘積      
    %end
   
    %------------------------------------------------
    % 上面是原始算法，下面是優化算法
    
    for j = 1:index_rows
        xn_rows_index = index(j,:);
        xn_rows = xn(xn_rows_index,:);
        un(j,:) = prod(xn_rows,1);
    end
    
    Un = [Un;un];       % FIR 濾波器的抽頭輸入信號矢量 Un
    clear un;
end
      
len_filter = size(Un,1);