comp_dgt_fb.m

linear time－frequency toolbox

function [coef]=comp_dgt_fb(f,g,a,M)
%COMP_DGT_FB  Filter bank DGT
%   Usage:  c=comp_dgt_fb(f,g,a,M);
%  
%   This is a computational routine. Do not call it directly.
%
%   See help on DGT.

%   Author : Peter Soendergaard.

% Calculate the parameters that was not specified.
L=size(f,1);
b=L/M;
N=L/a;
gl=length(g);
W=size(f,2);      % Number of columns to apply the transform to.

g=fftshift(g);

glhalf=ceil(gl/2);

coef=zeros(M*N,W);

% Calculate indexing needed for the splitting.
no_split_s=ceil((gl/2/a));
no_split_e=floor((L-gl/2)/a);

nrange=(0:N-1)';
spr=rem(nrange*a-floor(gl/2)+L,L)+1;       % Starting point
epr=rem(nrange*a-floor(gl/2)+gl-1+L,L)+1;  % ending point.
nsplitr=L+1-spr;

% Initialize fp_g, because first reference might be slicing.
fp_g=zeros(gl,W);

% Replicate g when multiple columns should be transformed.
gw=repmat(g,1,W);

% This covers the the cases when splitting is necessary:
Nsplit=[0:no_split_s-1,no_split_e+1:N-1];
% This n is one-indexed, to avoid to many +1
for n=Nsplit+1

  % Multiply the correct part of the signal by g
  fp_g(1:nsplitr(n),:)    =f(spr(n):L,:).*gw(1:nsplitr(n),:);
  fp_g(nsplitr(n)+1:gl,:) =f(1:epr(n),:).*gw(nsplitr(n)+1:gl,:);

  % Periodize the result prior to the fft.
  fpp=zeros(M,W);
  for ii=0:gl/M-1
    fpp=fpp+fp_g(ii*M+1:(ii+1)*M,:);
  end;
  
  % Shift back again.
  coef((n-1)*M+1:n*M,:)=circshift(fpp,(n-1)*a-glhalf);
    
end;

% This cover the cases when splitting is not necessary:
% This n is one-indexed, to avoid to many +1
for n=no_split_s+1:no_split_e+1
  fp_g=f(spr(n):epr(n),:).*gw;

  % Periodize the result prior to the fft.
  fpp=zeros(M,W);
  for ii=0:gl/M-1
    fpp=fpp+fp_g(ii*M+1:(ii+1)*M,:);
  end;
  
  % Shift back again.
  coef((n-1)*M+1:n*M,:)=circshift(fpp,(n-1)*a-glhalf);

end;

coef=reshape(coef,M,N,W);
coef=fft(coef);


% For later reference.
%
% The "fast" summation works using
%
%
% ind_per=reshape(1:gl,M,gl/M);
% fpp=sum(fp_g(ind_per),2);