function [ S, F_center, F_start, F_end, Freq_Resp, Trans_Band, S_element ] = CCDST(x, N0, b, r)

% Function to compute the CC-DST of a given signal.

% x: signal to be analyzed by the CC-DST.
% N0: number of frequency samples in the first spectral subband (without those in the transition bands).
% b: parameter to control the frequency spacing between the spectral subbands (b>1).
% r: roll-off factor for the asymmetric raised-cosine filters (r=0 gives brick-wall filters).
% S: cell array contains the CC-DST voices.
% F_Center: matrix contains the voices' center frequencies.
% F_start: matrix contains lower-edge frequencies of the voices.
% F_end: matrix contains the upper-edge frequency of the voices.
% Freq_Resp: cell array contains the frequency responses of the filters.
% Trans_Band: matrix contains the number of frequency samples of the left-hand transition bands.
% S_element: number of the time-frequency coefficients returned by the CC-DST.

% example: CCDST(randn(1,1024), 64, 1.2, 0.8)

if N0<2
    error('N0 must be larger than or equal to 2.')
end
if b<=1
    error('b must be larger than 1.')
end

if iscolumn(x)
    x = x.';
end

N  = length(x);
if mod(N, 2) 
    x(end) = [];
    N = N-1;
end

Xf = fftshift(fft(x));
 
f0  = @(b,j)round(b^(j-1));           
j0 = 1 + log((N0-1)/(b-1))/log(b);     % identity of first partition
jend = log((f0(b, j0) + N/2))/log(b);  % identity of last partition
j = j0:round(jend-j0)+j0+1;
 
for jj=1:length(j)
     
    [ Freq_Resp, f_center, f_start, f_end, Trans_Band_L, Trans_Band_R] = FreqPartition(b, r, j(jj), j0);   
    f_start_pos = f_start + N/2 +1; % freq index
    f_end_pos = f_end + N/2 +1;  
    
    if f_start_pos<=N % this condition is because depending on N and b, it might be enough to have max(j)=round(jend-j0)+j0
        if f_end_pos<=N
            Xf_part = Xf(f_start_pos:f_end_pos);
        else
            Xf_part = [Xf(f_start_pos:N), Xf(1:mod(f_end_pos,N))];
        end
        Xf_wind = Xf_part.*Freq_Resp;
        NormFactor = length(Xf_wind)/N;
        S_pos{jj} = NormFactor * ifft(ifftshift(Xf_wind)); 
        F_center_pos(jj) = f_center;
        F_start_pos(jj) = f_start;
        F_end_pos(jj) = f_end;
        Freq_Resp_pos{jj} = Freq_Resp;
        Trans_Band_pos(jj) = Trans_Band_L;
               
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Negative Freq %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Note: this part can be skipped if the signal is real, thanks to the CC-DST complex conjugate symmetry 
         
        Freq_Resp = flip(Freq_Resp);  % to mirror the freq response
        f_start_neg = -f_end + N/2 + 1;
        f_end_neg = -f_start + N/2 + 1;
        if f_start_neg>0
            Xf_part = Xf(f_start_neg:f_end_neg);
        elseif f_start_neg<0
            Xf_part = [Xf(max(1,mod(f_start_neg, N)):N), Xf(1:f_end_neg)];
        else
            Xf_part = [Xf(end), Xf(1:f_end_neg)];
        end     
        Xf_wind = Xf_part.*Freq_Resp;    
        if mod(length(Xf_wind), 2)== 0  % to preserve Fourier symmetry
            Xf_wind = circshift(Xf_wind, [0,1]);
        end
        NormFactor = length(Xf_wind)/N;
        S_neg{jj} = NormFactor * ifft(ifftshift(Xf_wind));
        F_center_neg(jj)= -f_center;
        F_start_neg(jj) = -f_end;
        F_end_neg(jj) = -f_start;
        Freq_Resp_neg{jj} = Freq_Resp;
        Trans_Band_neg(jj) = Trans_Band_R;  
    end
end

S = [flip(S_neg), S_pos]; 
S_element = sum(cellfun('length',S));
F_center = [flip(F_center_neg), F_center_pos];
F_start = [flip(F_start_neg), F_start_pos];
F_end = [flip(F_end_neg), F_end_pos];
Freq_Resp = [flip(Freq_Resp_neg), Freq_Resp_pos];
Trans_Band = [flip(Trans_Band_neg), Trans_Band_pos];

end