[931198]: filterbank / filterbankrealbounds.m  Maximize  Restore  History

Download this file

93 lines (74 with data), 2.4 kB

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
function [AF,BF]=filterbankrealbounds(g,a,L);
%FILTERBANKREALBOUNDS Frame bounds of filter bank for real signals only
% Usage: fcond=filterbankrealbounds(g,a);
% [A,B]=filterbankrealbounds(g,a);
%
% `filterbankrealbounds(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. Use this
% function on the common construction where the filters in *g* only covers
% the positive frequencies.
%
% `[A,B]=filterbankrealbounds(g,a)` returns the lower and upper frame
% bounds explicitly.
%
% See also: filterbank, filterbankdual
if nargin<3
error('%s: Too few input parameters.',upper(mfilename));
end;
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;
if all(a==a(1))
% Uniform filterbank, use polyphase representation
a=a(1);
N=L/a;
% 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;
Ha=zeros(a,M,thisclass);
Hb=zeros(a,M,thisclass);
for w=0:N-1
idx_a = mod(w-(0:a-1)*N,L)+1;
idx_b = mod((0:a-1)*N-w,L)+1;
Ha = G(idx_a,:);
Hb = conj(G(idx_b,:));
% A 'real' is needed here, because the matrices are known to be
% Hermitian, but sometimes Matlab/Octave does not recognize this.
work=real(eig(real(Ha*Ha'+Hb*Hb')));
AF=min(AF,min(work));
BF=max(BF,max(work));
end;
AF=AF/a;
BF=BF/a;
else
if info.ispainless
% Compute the diagonal of the frame operator.
f=filterbankresponse(g,a,L,'real');
AF=min(f);
BF=max(f);
else
error(['%s: There is no fast method to find the frame bounds of ' ...
'this filterbank as it is neither uniform nor painless. ' ...
'Please see FRAMEBOUNDS for an iterative method that can ' ...
'solve the problem.'],upper(mfilename));
end;
end;
if nargout<2
% Avoid the potential warning about division by zero.
if AF==0
AF=Inf;
else
AF=BF/AF;
end;
end;

Get latest updates about Open Source Projects, Conferences and News.

Sign up for the SourceForge newsletter:





No, thanks