1function [AF,BF]=filterbankrealbounds(g,a,L);
2%-*- texinfo -*-
3%@deftypefn {Function} filterbankrealbounds
4%@verbatim
5%FILTERBANKREALBOUNDS  Frame bounds of filter bank for real signals only
6%   Usage: fcond=filterbankrealbounds(g,a,L);
7%          [A,B]=filterbankrealbounds(g,a,L);
8%          [...]=filterbankrealbounds(g,a);
9%
10%   FILTERBANKREALBOUNDS(g,a,L) calculates the ratio B/A of the frame
11%   bounds of the filterbank specified by g and a for a system of length
12%   L. The ratio is a measure of the stability of the system.  Use this
13%   function on the common construction where the filters in g only covers
14%   the positive frequencies.
15%
16%   FILTERBANKREALBOUNDS(g,a) does the same, but the filters must be FIR
17%   filters, as the transform length is unspecified. L will be set to
18%   next suitable length equal or bigger than the longest impulse response
19%   such that L=filterbanklength(gl_longest,a).
20%
21%   [A,B]=FILTERBANKREALBOUNDS(g,a) returns the lower and upper frame
22%   bounds explicitly.
23%
24%@end verbatim
25%@strong{Url}: @url{http://ltfat.github.io/doc/filterbank/filterbankrealbounds.html}
26%@seealso{filterbank, filterbankdual}
27%@end deftypefn
28
29% Copyright (C) 2005-2016 Peter L. Soendergaard <peter@sonderport.dk>.
30% This file is part of LTFAT version 2.3.1
31%
32% This program is free software: you can redistribute it and/or modify
33% it under the terms of the GNU General Public License as published by
34% the Free Software Foundation, either version 3 of the License, or
35% (at your option) any later version.
36%
37% This program is distributed in the hope that it will be useful,
38% but WITHOUT ANY WARRANTY; without even the implied warranty of
39% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
40% GNU General Public License for more details.
41%
42% You should have received a copy of the GNU General Public License
43% along with this program.  If not, see <http://www.gnu.org/licenses/>.
44
45complainif_notenoughargs(nargin,2,'FILTERBANKREALBOUNDS');
46if nargin<3
47    L = [];
48end
49
50[g,asan,info]=filterbankwin(g,a,L,'normal');
51if isempty(L)
52    if info.isfir
53        % Pick shortest possible length for FIR filterbank
54        L = filterbanklength(info.longestfilter,asan);
55    else
56        % Just thow an error, nothing reasonable can be done without L
57        error(['%s: L must be specified when not working with FIR ',...'
58               'filterbanks.'], upper(mfilename));
59    end
60end
61M=info.M;
62
63if L~=filterbanklength(L,asan)
64    error(['%s: Specified length L is incompatible with the length of ' ...
65           'the time shifts.'],upper(mfilename));
66end;
67
68AF=Inf;
69BF=0;
70
71% Prioritize painless over uniform algorithm
72if info.isuniform && info.ispainless
73    info.isuniform = 0;
74end
75
76if info.isuniform
77  % Uniform filterbank, use polyphase representation
78  a=a(1);
79
80  N=L/a;
81
82  % G1 is done this way just so that we can determine the data type.
83  G1=comp_transferfunction(g{1},L);
84  thisclass=assert_classname(G1);
85  G=zeros(L,M,thisclass);
86  G(:,1)=G1;
87  for ii=2:M
88    G(:,ii)=cast(comp_transferfunction(g{ii},L),thisclass);
89  end;
90
91  Ha=zeros(a,M,thisclass);
92  Hb=zeros(a,M,thisclass);
93
94  for w=0:N-1
95    idx_a = mod(w-(0:a-1)*N,L)+1;
96    idx_b = mod((0:a-1)*N-w,L)+1;
97    Ha = G(idx_a,:);
98    Hb = conj(G(idx_b,:));
99
100    % A 'real' is needed here, because the matrices are known to be
101    % Hermitian, but sometimes Matlab/Octave does not recognize this.
102    work=real(eig(real(Ha*Ha'+Hb*Hb')));
103
104    AF=min(AF,min(work));
105    BF=max(BF,max(work));
106
107  end;
108
109  AF=AF/a;
110  BF=BF/a;
111
112else
113    if info.ispainless
114        % Compute the diagonal of the frame operator.
115        f=comp_filterbankresponse(g,asan,L,1);
116
117        AF=min(f);
118        BF=max(f);
119    else
120        error(['%s: There is no fast method to find the frame bounds of ' ...
121               'this filterbank as it is neither uniform nor painless. ' ...
122               'Please see FRAMEBOUNDS for an iterative method that can ' ...
123               'solve the problem.'],upper(mfilename));
124    end;
125end;
126
127if nargout<2
128  % Avoid the potential warning about division by zero.
129  if AF==0
130    AF=Inf;
131  else
132    AF=BF/AF;
133  end;
134end;
135
136
137
138