[34f3ea]: filterbank / filterbankrealdual.m  Maximize  Restore  History

Download this file

87 lines (67 with data), 2.1 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
function gdout=filterbankrealdual(g,a,varargin);
%FILTERBANKREALDUAL Dual filters of filterbank for real signals only
% Usage: gd=filterbankdual(g,a);
%
% `filterabankdual(g,a)` computes the canonical dual filters of *g* for a
% channel subsampling rate of *a* (hop-size). The dual filters work only
% for real-valued signals. Use this function on the common construction
% where the filters in *g* only covers the positive frequencies.
%
% The format of the filters *g* are described in the
% help of |filterbank|.
%
% To actually invert the output of a filterbank, use the dual filters
% together with `2*real(ifilterbank(...))`.
%
% See also: filterbank, ufilterbank, ifilterbank
if nargin<2
error('%s: Too few input parameters.',upper(mfilename));
end;
definput.keyvals.L=[];
[flags,kv,L]=ltfatarghelper({'L'},definput,varargin);
[g,info]=filterbankwin(g,a,L,'normal');
M=info.M;
if (~isempty(L)) && (L~=filterbanklength(L,a))
error(['%s: Specified length L is incompatible with the length of ' ...
'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);
% 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;
% This is the original code
%for k=0:a-1
% Ha(k+1,:) = G(mod(w-k*N,L)+1,:);
% Hb(k+1,:) = conj(G(mod(k*N-w,L)+1,:));
%end;
gd=zeros(N,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,:));
Ha=(Ha*Ha'+Hb*Hb')\Ha;
gd(idx_a,:)=Ha;
end;
gd=ifft(gd)*a;
if isreal(g)
gd=real(gd);
end;
gdout=cell(1,M);
for m=1:M
gdout{m}=cast(gd(:,m),thisclass);
end;
else
error('Not implemented yet.');
end;

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

Sign up for the SourceForge newsletter:





No, thanks