Skip to content

Commit

Permalink
Simplified backend code in the filterbank routines
Browse files Browse the repository at this point in the history 
A duplicate code was moved from filterbankwin (comp_fourierwindow) to
comp_filterbank_pre. filterbankwin now no longer instantizes g.H(L), and
requires L just in special cases.
Related routines can now handle more filter formats like 'gauss',
{'hann',16} etc.
  • Loading branch information
@susnak
susnak committed Nov 28, 2014
1 parent e89d896 commit f22afa0
Show file tree
Hide file tree
Showing 19 changed files with 284 additions and 337 deletions.
4 changes: 2 additions & 2 deletions blockproc/blockframeaccel.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,9 +101,9 @@
case {'dwilt','wmdct'}
[~, info] = wilwin(F.g,F.M,[],upper(mfilename));
case {'filterbank','ufilterbank'}
[~, info] = filterbankwin(F.g,F.a);
[~, ~,info] = filterbankwin(F.g,F.a);
case {'filterbankreal','ufilterbankreal'}
[~, info] = filterbankwin(F.g,F.a,'real');
[~, ~,info] = filterbankwin(F.g,F.a,'real');
case 'fwt'
winLen = (F.g.a(1)^F.J-1)/(F.g.a(1)-1)*(numel(F.g.g{1}.h)-1)+1;
end;
Expand Down
44 changes: 33 additions & 11 deletions comp/comp_filterbank_pre.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,8 +23,7 @@
mFreq(mTime) = [];

% Prepare time-domain filters
for mId=1:numel(mTime)
m = mTime(mId);
for m = mTime

% Handle .fc parameter
if isfield(g{m},'fc') && g{m}.fc~=0
Expand All @@ -41,12 +40,12 @@
% Do zero padding when the offset is big enough so the initial imp. resp.
% support no longer intersects with zero
if g{m}.offset > 0
g{m}.h = [zeros(g{m}.offset,1);g{m}.h(:)];
g{m}.h = [zeros(g{m}.offset,1,class(g{m}.h));g{m}.h(:)];
g{m}.offset = 0;
end

if g{m}.offset < -(numel(g{m}.h)-1)
g{m}.h = [g{m}.h(:);zeros(-g{m}.offset-numel(g{m}.h)+1,1)];
g{m}.h = [g{m}.h(:);zeros(-g{m}.offset-numel(g{m}.h)+1,1,class(g{m}.h))];
g{m}.offset = -(numel(g{m}.h)-1);
end

Expand All @@ -57,20 +56,43 @@
for m=mFreq

if isfield(g{m},'h')
g{m}.H=comp_transferfunction(g{m},L);
tmpg = circshift(postpad(g{m}.h,L),g{m}.offset);
g{m}=rmfield(g{m},'h');
g{m}=rmfield(g{m},'offset');
if isfield(g{m},'fc'), g{m}=rmfield(g{m},'fc'); end;
if isfield(g{m},'fc')
l=(0:L-1).'/L;
tmpg = tmpg.*exp(2*pi*1i*round(g{m}.fc*L/2)*l);
g{m}=rmfield(g{m},'fc');
end;
g{m}.H = fft(tmpg);
% The following parameters have to be set to zeros, because they
% have already been incorporated in the freq. resp. calculation.
g{m}.foff = 0;
% Store the length used
g{m}.L = L;
elseif isfield(g{m},'H') && ~isnumeric(g{m}.H)
g{m}.H=g{m}.H(L);
g{m}.foff=g{m}.foff(L);
% Store the length used
g{m}.L = L;
elseif isfield(g{m},'H')
if isnumeric(g{m}.H)
if isfield(g{m},'L')
if g{m}.L~=L
% middlepad in the time domain. This will break
% g.H = fft(middlepad(ifft(circshift(postpad(g.H(:),g.L),g.foff)),L));
% g.foff = 0;
% g.L = L;
error(['%s: g.H was already instantialized with L=%i, but ',...
'it is now used with L=%i.'],upper(mfilename),g.L,L);
end
else
% We do not know which L was g.H created with, there is no way
% how to handle this properly.
error('%s: g.H is already a numeric vector, but g.L was not specified.',...
upper(mfilename));
end
else
g{m}.H=g{m}.H(L);
g{m}.foff=g{m}.foff(L);
% Store the length used
g{m}.L = L;
end
end

if isfield(g{m},'H') && isfield(g{m},'delay') && g{m}.delay~=0
Expand Down
40 changes: 11 additions & 29 deletions comp/comp_fourierwindow.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@
%
% `[g,info]=comp_fourierwindow(g,L,callfun)` will compute the window
% from a text description or a cell array containing additional
% parameters.
% parameters.
%
% See also: gabwin, wilwin

Expand Down Expand Up @@ -101,7 +101,7 @@
end
elseif isfield(g,'H') && ...
( isnumeric(g.H) && isvector(g.H) || isa(g.H,'function_handle') )

info.wasreal=isfield(g,'realonly') && g.realonly;
info.gl=[];

Expand All @@ -116,34 +116,14 @@
end

if isa(g.H,'function_handle')
if isempty(L)
error(['%s: You must specify a length L if a window is ',...
'represented in a fourier domain as a function ',...
'handle'],upper(mfilename));
end
g.H=g.H(L);
if ~isa(g.foff,'function_handle')
error('%s: g.foff should be a function handle.',...
upper(mfilename));
callfun);
end
g.foff=g.foff(L);
% And store the length the freq. resp. was instantiated
% with.
g.L = L;
else
% If g.H is already numeric, we must handle transforming it to a
% different length.
if isfield(g,'L')
if g.L~=L
g.H = fft(middlepad(ifft(circshift(postpad(g.H(:),g.L),g.foff)),L));
g.foff = 0;
g.L = L;
end
else
% We do not know which L was g.H created with, there is no way
% how to handle this properly.
error('%s: No way of knowing which L was used for g.H.',...
upper(mfilename));
elseif isnumeric(g.H)
if ~isfield(g,'L')
error('%s: .H is numeric, but g.L was not defined.',...
callfun);
end
end

Expand All @@ -153,6 +133,7 @@
'anonymous fcn) fields.'],callfun);
end;
else
% We can probably never get here
% Information to be determined post creation.
info.wasreal = isreal(g);
info.gl = length(g);
Expand All @@ -168,8 +149,9 @@
function complain_L(L,callfun)

if isempty(L)
error(['%s: You must specify a length L if a window is represented as a ' ...
'text string or cell array.'],callfun);
error(['L:undefined',...
'%s: You must specify a length L if a window is represented as a ' ...
'text string, cell array or anonymous function.'],callfun);
end;


4 changes: 2 additions & 2 deletions comp/comp_painlessfilterbank.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,12 +2,12 @@
%COMP_PAINLESSFILTERBANK
%
% Function computes filterbank dual or tight frame for the painless case.
% All g{ii}.H should already be numeric verctors.
%

M = numel(g);

F=comp_filterbankresponse(g,a,L,do_real);
F = comp_filterbankresponse(g,a,L,do_real);
g = comp_filterbank_pre(g,a,L,0);

if strcmpi(type,'tight')
F = sqrt(F);
Expand Down
54 changes: 10 additions & 44 deletions comp/comp_transferfunction.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,48 +6,14 @@
% handle filters in a proper internal format i.e. already processed by
% |filterbankwin|.

l=(0:L-1).'/L;
if isfield(g,'h')
g_time=circshift(postpad(g.h,L),g.offset);

if isfield(g,'fc')
g_time = g_time.*exp(2*pi*1i*round(g.fc*L/2)*l);
end

if isfield(g,'realonly') && g.realonly
g_time=real(g_time);
end;

H=fft(g_time);
elseif isfield(g,'H')
if ~isnumeric(g.H)
g.H=g.H(L);
g.foff=g.foff(L);
else
% If g.H is already numeric, we must handle transforming it to a
% different length.
if isfield(g,'L')
if g.L~=L
g.H = fft(middlepad(ifft(circshift(postpad(g.H(:),g.L),g.foff)),L));
g.foff = 0;
end
else
% We do not know which L was g.H created with, there is no way
% how to handle this properly.
error('%s: No way of knowing which L was used for g.H.',...
upper(mfilename));
end
end;

H=circshift(postpad(g.H(:),L),g.foff);

if isfield(g,'delay')
H = H.*exp(-2*pi*1i*round(g.delay)*l);
end

if isfield(g,'realonly') && g.realonly
H=(H+involute(H))/2;
end;
else
error('%s: Unrecognized filter format. The struct should have either .h or .H field.',upper(mfilename));
% Setting crossover to 0 ensures FIR filters to be transformed to
% full-length Frequency-domain defined filters with g.H and g.foff fields.
g = comp_filterbank_pre({g},1,L,0);
% Band-limited filters have to be made full-length
H = circshift(postpad(g{1}.H(:),L),g{1}.foff);

% Realonly has to be treated separatelly for band-limited filters
if isfield(g,'realonly') && g.realonly
H=(H+involute(H))/2;
end;

22 changes: 10 additions & 12 deletions filterbank/filterbank.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,23 +47,21 @@
end;

if isempty(L)
L=filterbanklength(Ls,a);
L=filterbanklength(Ls,a);
end;

[g,info]=filterbankwin(g,a,L,'normal');
[g,asan]=filterbankwin(g,a,L,'normal');

if size(a,1)>1
if size(a,1)~=info.M
error(['%s: The number of entries in "a" must match the number of ' ...
'filters.'],upper(mfilename));
end;
end;

% if size(a,1)>1
% if size(a,1)~= numel(g);
% error(['%s: The number of entries in "a" must match the number of ' ...
% 'filters.'],upper(mfilename));
% end;
% end;

f=postpad(f,L);

g=comp_filterbank_pre(g,info.a,L,kv.crossover);

c=comp_filterbank(f,g,info.a);
g=comp_filterbank_pre(g,asan,L,kv.crossover);
c=comp_filterbank(f,g,asan);


43 changes: 30 additions & 13 deletions filterbank/filterbankbounds.m
View file
Original file line number Diff line number Diff line change
@@ -1,29 +1,46 @@
function [AF,BF]=filterbankbounds(g,a,L);
%FILTERBANKBOUNDS Frame bounds of filter bank
% Usage: fcond=filterbankbounds(g,a);
% [A,B]=filterbankbounds(g,a);
function [AF,BF]=filterbankbounds(g,a,L)
%FILTERBANKBOUNDS Frame bounds of a filterbank
% Usage: fcond=filterbankbounds(g,a,L);
% [A,B]=filterbankbounds(g,a,L);
% [...]=filterbankbounds(g,a);
%
% `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,L)` returns the lower and upper frame bounds
% `filterbankbounds(g,a)` does the same, but the filters must be FIR
% filters, as the transform length is unspecified. *L* will be set to
% next suitable length equal or bigger than the longest impulse response
% such that `L=filterbanklength(gl_longest,a)`.
%
% `[A,B]=filterbankbounds(...)` returns the lower and upper frame bounds
% explicitly.
%
% See also: filterbank, filterbankdual

complainif_notenoughargs(nargin,2,'FILTERBANKBOUNDS');
if nargin<3
error('%s: Too few input parameters.',upper(mfilename));
end;
L = [];
end

if L~=filterbanklength(L,a)
[g,asan,info]=filterbankwin(g,a,L,'normal');
if isempty(L)
if info.isfir
% Pick shortest possible length for FIR filterbank
L = filterbanklength(info.longestfilter,asan);
else
% Just thow an error, nothing reasonable can be done without L
error(['%s: L must be specified when not working with FIR ',...'
'filterbanks.'], upper(mfilename));
end
end
M=info.M;

if L~=filterbanklength(L,asan)
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;

Expand Down Expand Up @@ -69,7 +86,7 @@

if info.ispainless
% Compute the diagonal of the frame operator.
f=comp_filterbankresponse(g,info.a,L,0);
f=comp_filterbankresponse(g,asan,L,0);

AF=min(f);
BF=max(f);
Expand All @@ -84,7 +101,7 @@

if nargout<2
% Avoid the potential warning about division by zero.
if AF==0
if AF==0
AF=Inf;
else
AF=BF/AF;
Expand Down
Loading

0 comments on commit f22afa0

Please sign in to comment.