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