FOST_2D_DEC(S_input)

function D2_coeff = FOST_2D_DEC(S_input)
    % 2-Dimensional DOST:
%     tic
    si = size(S_input);
    NN = max(si);
    partitions = zeros(1,NN);
    D2_coeffM = zeros(si);
    D2_coeff = zeros(si); 
    
    %rows first
    for ii = 1:NN
        timeseries = S_input(ii, :);
        D2_coeffM(ii,:) = FOST_1D_DEC(timeseries);
    end

    for jj = 1:NN
        timeseries = D2_coeffM(:, jj);
        D2_coeff(:, jj) = FOST_1D_DEC(timeseries);
    end
%     Decomposition_Time = toc
end

function bas2 = FOST_1D_DEC(S_input)

    si = size(S_input);
    NN = max(si);
    [vv, bb, tt, num_it, compact_b] = get_vbt(NN);
    
    %Calculate the FFT
    fftdata = transpose(fftshift(fft(S_input)));

%     tic
    %Use the inverse FFT to calculate the DOST coefficients
    bas2 = zeros(si); 
    init_p = 0;
    coof = -pi*i;
    for ll = 1:num_it
        step = compact_b(ll);
        %Construct the corresponding Ramp matrix
        Ramp = zeros(1, step);
        for ii = 1:step
%             vv(init_p + ii);
            Ramp(ii) = (-1)^ii;
%             Ramp(ii) = exp((ii-1)*coof*(1-(2*vv(init_p + ii)-step)/step));            
        end
        temp_matrix = ifft(fftdata(init_p+1:init_p+step))*sqrt(step);
        if size(Ramp)==size(temp_matrix)
            bas2(init_p+1:init_p+step) = Ramp.*temp_matrix;
        else
            bas2(init_p+1:init_p+step) = Ramp.*transpose(temp_matrix);
        end
        clear temp_matrix;
        clear Ramp;
        init_p = init_p + step;
    end
%     FOST_decom_time = toc
end


function [v, b, t, num, com_b] = get_vbt(NN);
    count=0;
    end_P = log2(NN)-1;
    b = zeros(1,NN);
    v = zeros(1,NN);
    t = zeros(1,NN);
    com_b = zeros(1, 2*end_P + 2);
    num = 0;
        
    
    %%%%%%NNegative frequency
    
    count = count+1;
    v(count) = -NN/2 + .5; b(count) = 1; t(count) = 0;    
    num = num + 1;
    com_b(num) = b(count);
        
    for p = end_P:-1:2
        freq = -2^(abs(p)-2)*3 + 1;
        beta = 2^(abs(p)-1);
        for tau = beta-1:-1:0
%         for tau = 0: beta-1
            count = count+1;
            v(count) = freq;    b(count) = beta;   t(count) = tau;
        end
        num = num + 1;
        com_b(num) = beta;
    end
    
    count=count+1;
    v(count) = -0.5; b(count) = 1; t(count) = 0;
    num = num + 1;    
    com_b(num) = b(count);
    
    %%%%%%%Positive frequency
    count=count+1;
    v(count) = 0.5; b(count) = 1; t(count) = 0;
    num = num + 1;
    com_b(num) = b(count);
    
    count=count+1;
    v(count) = 1.5; b(count) = 1; t(count) = 0;
    num = num + 1;    
    com_b(num) = b(count);
    
    for p=2:end_P
        freq = 2^(abs(p)-2)*3;
        beta = 2^(abs(p)-1);        
        for tau=0:beta-1
            count=count+1;
            v(count) = freq;    b(count) = beta;   t(count) = tau;
        end
        num = num + 1;
        com_b(num) = beta;
    end   
end