function [Y,H] = dcells(X,zpf,W,aliasbands,debug) % DCELLS - Break up complex spectrum X into dyadic cells % Origin: http://ccrma.stanford.edu/~jos/sasp/dcells.m % % Reference: % http://ccrma.stanford.edu/~jos/sasp/Audio_Filter_Banks.html % % Example Usage: % http://ccrma.stanford.edu/~jos/sasp/Spectral_Rotation_Real_Signals.html if nargin<5, debug=0; end if nargin<4, aliasbands=1; end if nargin<3, W=1; end % (1 means "impulse in freq domain") if nargin<2, zpf=1; end if zpf>1 && W==1 error('dcells: Must compute asinc kernel for given zero-padding factor'); end X = X(:); W = W(:)/sum(W); N0zp = length(X); N0 = N0zp/zpf; NS = round(log2(N0)); % Number of octaves Y = cell(1,NS+1); % extra cell for lowest band H = zeros(NS+1,N0zp); % collection of filters (spectral windows) N = N0 * 2 .^ -(1:NS+1); % Spectral widths in (critically sampled) bins N(NS+1)=N(NS); % dc channel looks like next band up Nzp = N*zpf; for k=1:NS+1 if k==NS+1, Hideal = []; else Hideal = zeros(Nzp(k),1); end Hideal = [Hideal;ones(Nzp(k),1)]; if k>1, Hideal = [Hideal;zeros(N0zp-length(Hideal),1)]; end if W==1 && zpf==1 && aliasbands!=0 % no overlap-add decomposition kernel specified, so % simply partition the FFT bins dyadically: if k == (NS+1) Y{k} = X(1); % dc band else Y{k} = X(N(k)+1:2*N(k)) .'; % top half spectrum = top octave end H(k,:) = Hideal; else % Need to build an overlap-add spectral window. % Note that these filters H can be precomputed: H(k,:) = cconvr(W,Hideal,N0zp); % real circular convolution % cconvr() available as http://ccrma.stanford.edu/~jos/sasp/cconvr.m toshow = [W/W(1),Hideal,H(k,:).']; if debug>1 plot(toshow,'-*'); legend('Kernel','Hideal','H'); title(sprintf('Filter for Band %d (linear)',k)); grid('on'); pause; plot(db(toshow,-100),'-*'); grid('on'); legend('Kernel','Hideal','H'); ... title(sprintf('Filter for Band %d (dB)',k)); grid('on'); pause; end Yek = (X.') .* H(k,:); % Apply band filter as a spectral window % disp(sprintf("N(k) = %d",N(k))); if aliasbands Y{k} = alias(Yek,N0zp/Nzp(k)); % Aliasing onto a block of length Nzp(k) % alias() available as http://ccrma.stanford.edu/~jos/sasp/alias.m else Y{k} = Yek; % Filter bank retains perfect reconstruction property end if debug>2 plot(abs(Y{k}),'-*'); title(sprintf('|Spectral Band| %d',k)); grid('on'); pause; end end end Hsum = sum(H); fbserr = norm(Hsum-ones(size(Hsum))); disp(sprintf('Filter bank summation L2 error = %0.2f',fbserr)); if debug && fbserr > 1E-7 plot(Hsum,'-*'); title('Filter Bank Sum'); disp('PAUSING'); pause; end if debug>1 plot(H','-*',"linewidth",10); title('Complex Octave Filter Bank',"fontsize",18); axis([1, length(H), 0, max(max(abs(H)))*1.1]); bandslegend(NS+1,"'east'"); xlabel('Spectral Sample Number',"fontsize",18); ylabel('Filter Gain',"fontsize",18); disp('PAUSING'); pause; end