new support in pfilt and filterbank and friends for filters defined i…

…n the frequency domain. Added firfilter and blfilter.
Kevin2599 · Jun 7, 2013 · 34f3ea9 · 34f3ea9
1 parent 80bb6b7
commit 34f3ea9
Show file tree

Hide file tree

Showing 32 changed files with 869 additions and 285 deletions.
diff --git a/comp/arg_pfilt.m b/comp/arg_pfilt.m
@@ -0,0 +1,5 @@
+function definput=arg_pfilt(definput)
+
+  definput.keyvals.crossover=120;
+
+
diff --git a/comp/assert_sigreshape_pre.m b/comp/assert_sigreshape_pre.m
@@ -17,7 +17,7 @@
 
 
 if ~isnumeric(f)
-  error('%s: f must numeric.',callfun);
+  error('%s: The input must be numeric.',callfun);
 end;
 
 D=ndims(f);

diff --git a/comp/comp_fourierwindow.m b/comp/comp_fourierwindow.m
@@ -2,7 +2,7 @@
 %COMP_FOURIERWINDOW  Compute the window from numeric, text or cell array.
 %   Usage: [g,info] = comp_fourierwindow(g,a,M,L,wilson,callfun);
 %
-%   [g,info]=COMP_FOURIERWINDOW(g,L,callfun) will compute the window
+%   `[g,info]=comp_fourierwindow(g,L,callfun)` will compute the window
 %   from a text description or a cell array containing additional
 %   parameters.
 %
@@ -14,7 +14,6 @@
 
 % Default values.
 info.gauss=0;
-info.wasrow=0;
 info.isfir=0;
 
 % Manually get the list of window names
@@ -69,19 +68,40 @@
     else
       % g was a row vector.
       g=g(:);
-      info.wasrow=1;
     end;
   end;
-end;
-
-% Information to be determined post creation.
-info.wasreal = isreal(g);
-info.gl      = length(g);
+  g_time=g;
+  g=struct();
+  g.h=fftshift(g_time);
+  info.gl=numel(g_time);
+  g.offset=-floor(info.gl/2);  
+  g.centre=0;
+  g.realonly=0;
+  info.wasreal=isreal(g.h);
+else
 
-if (~isempty(L) && (info.gl<L))
-  info.isfir=1;
+    if isstruct(g)
+        if isfield(g,'h')
+            info.wasreal=isreal(g.h);
+            info.gl=length(g.h);
+            info.isfir=1;
+        else
+            info.wasreal=g.realonly;
+            info.gl=[];        
+        end;    
+    else
+        % Information to be determined post creation.
+        info.wasreal = isreal(g);
+        info.gl      = length(g);
+
+        if (~isempty(L) && (info.gl<L))
+            info.isfir=1;
+        end;
+
+    end;
+
 end;
-
+    
 function complain_L(L,callfun)
 
   if isempty(L)

diff --git a/comp/comp_pfilt.m b/comp/comp_pfilt.m
@@ -0,0 +1,74 @@
+function h=comp_pfilt(f,g,a,do_time);
+%COMP_PFILT  Compute filtering
+%
+%   If *do_time* is 1, the routine will use the time-side algorithm for
+%   FIR filters, otherwise it will always do the multiplication in the
+%   frequency domain.
+
+[L,W]=size(f);
+N=L/a;
+
+h=zeros(N,W,assert_classname(f));
+
+l=(0:L-1).'/L;
+
+if isfield(g,'h') && do_time
+    % Use a direct algorithm
+    g_time=circshift(postpad(g.h,L),g.offset).*...
+           exp(2*pi*1i*round(g.centre*L/2)*l);            
+    if g.realonly
+        g_time=real(g_time);
+    end;
+    g_time=conj(involute(g_time));
+    for n=0:N-1
+        h(n+1,:)=sum(bsxfun(@times,f,circshift(g_time,a*n)));
+    end;
+
+else
+
+    % Zero-extend and use a full length fft algorithm
+    % This case can be further optimized
+    G=comp_transferfunction(g,L);
+
+    for w=1:W
+        F=fft(f(:,w));
+        h(:,w)=ifft(sum(reshape(F.*G,N,a),2))/a;
+    end;
+
+end;
+
+
+% Old code
+if 0
+    if ~isstruct(g)
+        [g,info] = comp_fourierwindow(g,L,'PFILT');
+
+        h=squeeze(comp_ufilterbank_fft(f,g,a));
+
+        % FIXME: This check should be removed when comp_ufilterbank_fft{.c/.cc}
+        % have been fixed.
+        if isreal(f) && isreal(g)
+            h=real(h);
+        end;
+
+    else
+        N=L/a;
+        G=fir2long(g.filter(L),L);
+
+        for w=1:W
+            F=fft(f(:,w));
+            h(:,w)=ifft(sum(reshape(F.*G,N,a),2))/a;
+        end;
+
+        % Insert check for window being
+        if isreal(f) && g.isreal
+            h=real(h);
+        end;
+
+    end;
+end;
+
+
+
+
+
diff --git a/comp/comp_transferfunction.m b/comp/comp_transferfunction.m
@@ -0,0 +1,27 @@
+function H=comp_transferfunction(g,L)
+%COMP_TRANSFERFUNCTION  Compute the transfer function
+
+l=(0:L-1).'/L;
+if isfield(g,'h')
+
+    g_time=circshift(postpad(g.h,L),g.offset).*...
+           exp(2*pi*1i*round(g.centre*L/2)*l);
+
+    if g.realonly
+        g_time=real(g_time);
+    end;
+
+    H=fft(g_time);
+
+else
+    if ~isnumeric(g.H)
+        g.H=g.H(L);
+    end;
+
+    H=circshift(postpad(g.H,L),g.foff(L)).*exp(-2*pi*1i*round(g.delay)*l);
+
+    if g.realonly
+        H=(H+involute(H))/2;
+    end;
+
+end;
diff --git a/filterbank/filterbank.m b/filterbank/filterbank.m
@@ -36,6 +36,7 @@
   error('%s: Too few input parameters.',upper(mfilename));
 end;
 
+definput.import={'pfilt'};
 definput.keyvals.L=[];
 [flags,kv,L]=ltfatarghelper({'L'},definput,varargin);
 
@@ -67,17 +68,22 @@
 
 c=cell(M,1);
 for m=1:M
-  c{m}=zeros(N(m),W,assert_classname(f,g{1}));
+  c{m}=comp_pfilt(f,g{m},a(m),info.gl(m)<kv.crossover);
+  %zeros(N(m),W,assert_classname(f,g{1}));
 end;
 
-G=zeros(L,M,assert_classname(f,g{1}));
-for ii=1:M
-  G(:,ii)=fft(fir2long(g{ii},L));
-end;
+% Old, but more efficient code from before the struct filters
+if 0
+    G=zeros(L,M,assert_classname(f,g{1}));
+    for ii=1:M
+        G(:,ii)=fft(fir2long(g{ii},L));
+    end;
+
+    for w=1:W
+        F=fft(f(:,w),L);
+        for m=1:M
+            c{m}(:,w)=ifft(sum(reshape(F.*G(:,m),N(m),a(m)),2))/a(m);
+        end;
+    end;
 
-for w=1:W
-  F=fft(f(:,w),L);
-  for m=1:M
-    c{m}(:,w)=ifft(sum(reshape(F.*G(:,m),N(m),a(m)),2))/a(m);
-  end;
-end;
+end;
diff --git a/filterbank/filterbankbounds.m b/filterbank/filterbankbounds.m
@@ -1,30 +1,29 @@
-function [AF,BF]=filterbankbounds(g,a,varargin);
+function [AF,BF]=filterbankbounds(g,a,L);
 %FILTERBANKBOUNDS  Frame bounds of filter bank
 %   Usage: fcond=filterbankbounds(g,a);
 %          [A,B]=filterbankbounds(g,a);
 %
-%   `filterbankbounds(g,a)` calculates the ratio $B/A$ of the frame bounds of
-%   the filterbank specified by *g* and *a*. The ratio is a measure of the
-%   stability of the system.
+%   `filterbankbounds(g,a,L)` calculates the ratio $B/A$ of the frame bounds
+%   of the filterbank specified by *g* and *a* for a system of length
+%   *L*. The ratio is a measure of the stability of the system.
 %
-%   `[A,B]=filterbankbounds(g,a)` returns the lower and upper frame bounds
+%   `[A,B]=filterbankbounds(g,a,L)` returns the lower and upper frame bounds
 %   explicitly.
 %
 %   See also: filterbank, filterbankdual
 
-if nargin<2
+if nargin<3
   error('%s: Too few input parameters.',upper(mfilename));
 end;
 
-[a,M,longestfilter,lcm_a]=assert_filterbankinput(g,a);
-
-definput.keyvals.L=[];
-[flags,kv,L]=ltfatarghelper({'L'},definput,varargin);
-
-if isempty(L)
-  L=ceil(longestfilter/lcm_a)*lcm_a;
+if L~=filterbanklength(L,a)
+    error(['%s: Specified length L is incompatible with the length of ' ...
+           'the time shifts.'],upper(mfilename));
 end;
 
+[g,info]=filterbankwin(g,a,L,'normal');
+M=info.M;
+
 AF=Inf;
 BF=0;
 
@@ -33,13 +32,17 @@
   a=a(1);  
 
   N=L/a;
-
-  G=zeros(L,M,assert_classname(g{1}));
-  for ii=1:M
-    G(:,ii)=fft(fir2long(g{ii},L));
+
+  % G1 is done this way just so that we can determine the data type.
+  G1=comp_transferfunction(g{1},L);
+  thisclass=assert_classname(G1);
+  G=zeros(L,M,thisclass);
+  G(:,1)=G1;
+  for ii=2:M
+    G(:,ii)=comp_transferfunction(g{ii},L);
   end;
 
-  H=zeros(a,M,assert_classname(g{1}));
+  H=zeros(a,M,thisclass);
 
   for w=0:N-1
     idx = mod(w-(0:a-1)*N,L)+1;

diff --git a/filterbank/filterbankdual.m b/filterbank/filterbankdual.m
@@ -1,15 +1,15 @@
-function gdout=filterbankdual(g,a,varargin);
+function gdout=filterbankdual(g,a,L,varargin);
 %FILTERBANKDUAL  Dual filters
 %   Usage:  gd=filterbankdual(g,a);
 %           gd=filterbankdual(g,a,L);
 %
-%   `filterabankdual(g,a)` computes the canonical dual filters of *g* for a
+%   `filterbankdual(g,a)` computes the canonical dual filters of *g* for a
 %   channel subsampling rate of *a* (hop-size).
 %
 %   The input and output format of the filters *g* are described in the
 %   help of |filterbank|.
 %
-%   `filterabankdual(g,a,L)` computes canonical dual filters for a system
+%   `filterbankdual(g,a,L)` computes canonical dual filters for a system
 %   of length *L*. If *L* is not specified, the shortest possible
 %   transform length is choosen.
 %
@@ -22,33 +22,38 @@
   error('%s: Too few input parameters.',upper(mfilename));
 end;
 
-[a,M,longestfilter,lcm_a]=assert_filterbankinput(g,a);
+%definput.keyvals.L=[];
+%[flags,kv,L]=ltfatarghelper({'L'},definput,varargin);
 
-definput.keyvals.L=[];
-[flags,kv,L]=ltfatarghelper({'L'},definput,varargin);
+[g,info]=filterbankwin(g,a,L,'normal');
+M=info.M;
 
-if isempty(L)
-  L=ceil(longestfilter/lcm_a)*lcm_a;
-else
-  if rem(L,lcm_a)>0
+if (~isempty(L)) && (L~=filterbanklength(L,a))
     error(['%s: Specified length L is incompatible with the length of ' ...
-           'the time shifts. L = %i, lcm_a = %i'],upper(mfilename),L,lcm_a);
-  end;
+           'the time shifts.'],upper(mfilename));
 end;
 
-
 if all(a==a(1))
+
   % Uniform filterbank, use polyphase representation
+  if isempty(L)
+      error('%s: You need to specify L.',upper(mfilename));
+  end;
+
   a=a(1);
 
-  G=zeros(L,M,assert_classname(g{1}));
-  for ii=1:M
-    G(:,ii)=fft(fir2long(g{ii},L));
+  % G1 is done this way just so that we can determine the data type.
+  G1=comp_transferfunction(g{1},L);
+  thisclass=assert_classname(G1);
+  G=zeros(L,M,thisclass);
+  G(:,1)=G1;
+  for ii=2:M
+    G(:,ii)=comp_transferfunction(g{ii},L);
   end;
 
   N=L/a;
 
-  gd=zeros(N,M,assert_classname(g{1}));
+  gd=zeros(N,M,thisclass);
 
   for w=0:N-1
     idx = mod(w-(0:a-1)*N,L)+1;
@@ -67,11 +72,22 @@
 
   gdout=cell(1,M);
   for m=1:M
-    gdout{m}=cast(gd(:,m),assert_classname(g{1}));
+    gdout{m}=cast(gd(:,m),thisclass);
   end;
 
 else
-
-  error('Not implemented yet.');
-
+
+    F=filterbankframediag(g,a,L);
+
+    gdout=cell(1,M);
+    for m=1:M
+        thisgd=struct();
+        thisgd.H=comp_transferfunction(g{m},L)./F;
+        thisgd.foff=@(L) 0;
+        thisgd.realonly=0;
+        thisgd.delay=0;
+
+        gdout{m}=thisgd;
+    end;
+
 end;
Original file line number	Diff line number	Diff line change
		@@ -0,0 +1,5 @@
		function definput=arg_pfilt(definput)

		definput.keyvals.crossover=120;