Code covered by the BSD License

# 2-D Fast DOST Decompotion and Reconstruction

### Yanwei Wang (view profile)

30 Oct 2009 (Updated )

This is the 2-D Fast DOST Decomposition and reconstruction. The computational complexity is O(NlnN)

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```