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