Mentions légales du service

Skip to content
Snippets Groups Projects
Commit 699dee5b authored by hhakim's avatar hhakim
Browse files

Enable butterfly matrix optimization in matfaust.circ.

parent bddf85f9
No related branches found
No related tags found
No related merge requests found
...@@ -3,6 +3,8 @@ ...@@ -3,6 +3,8 @@
%> %>
%> @param c: the vector to define the circulant Faust. %> @param c: the vector to define the circulant Faust.
%> @param 'dev', str: 'gpu' or 'cpu' to create the Faust on CPU or GPU ('cpu' is the default). %> @param 'dev', str: 'gpu' or 'cpu' to create the Faust on CPU or GPU ('cpu' is the default).
%> @param 'diag_opt', logical: if true then the returned Faust is optimized using
%> matfaust.opt_butterfly_faust (because the DFT is used to implement circ).
%> %>
%> @b Example: %> @b Example:
%> %>
...@@ -58,6 +60,7 @@ function C = circ(c, varargin) ...@@ -58,6 +60,7 @@ function C = circ(c, varargin)
import matfaust.Faust import matfaust.Faust
dev = 'cpu'; dev = 'cpu';
argc = length(varargin); argc = length(varargin);
diag_opt = false;
if(argc > 0) if(argc > 0)
for i=1:2:argc for i=1:2:argc
if(argc > i) if(argc > i)
...@@ -71,6 +74,13 @@ function C = circ(c, varargin) ...@@ -71,6 +74,13 @@ function C = circ(c, varargin)
else else
dev = tmparg; dev = tmparg;
end end
case 'diag_opt'
if(argc == i || ~ islogical(tmparg))
error('diag_opt keyword argument is not followed by a logical')
else
diag_opt = tmparg;
end
otherwise otherwise
if((isstr(varargin{i}) || ischar(varargin{i})) && ~ strcmp(tmparg, 'cpu') && ~ startsWith(tmparg, 'gpu')) if((isstr(varargin{i}) || ischar(varargin{i})) && ~ strcmp(tmparg, 'cpu') && ~ startsWith(tmparg, 'gpu'))
error([ tmparg ' unrecognized argument']) error([ tmparg ' unrecognized argument'])
...@@ -92,22 +102,38 @@ function C = circ(c, varargin) ...@@ -92,22 +102,38 @@ function C = circ(c, varargin)
error('c must be a vector') error('c must be a vector')
end end
n = numel(c); n = numel(c);
F = matfaust.dft(n, 'normed', false); F = matfaust.dft(n, 'normed', false, 'diag_opt', false);
FH = F'; nf = numfactors(F);
P = factors(F, nf); % bitrev perm
if (size(c, 1) < size(c, 2)) if (size(c, 1) < size(c, 2))
c = c.'; c = c.';
end end
S = sparse(diag(FH*(c/n))); D = diag(F'*(c/n));
S = sparse(P * D * P');
FwP = left(F, nf - 1); % F without permutation
% C = F * matfaust.Faust(S*factors(FH, 1)) * right(FH, 2); % C = F * matfaust.Faust(S*factors(FH, 1)) * right(FH, 2);
nf = numfactors(F); if ~ matfaust.isFaust(FwP)
if(nf > 3) FwP = matfaust.Faust(FwP);
C = left(F, nf-1) * Faust(factors(F, nf) * S * factors(FH, 1) * factors(FH, 2)) * right(FH, 3); end
elseif(nf > 2) if diag_opt
C = left(F, nf-1) * Faust(factors(F, nf) * S * factors(FH, 1) * factors(FH, 2)) * Faust(right(FH, 3)); r = matfaust.opt_butterfly_faust(FwP)';
else else
C = Faust(left(F, nf-1)) * Faust(factors(F, nf) * S * factors(FH, 1) * factors(FH, 2)); r = FwP';
end
nfwp = numfactors(FwP);
l = left(FwP, nfwp-1);
if ~ matfaust.isFaust(l)
l = matfaust.Faust(l);
end end
l = l * matfaust.Faust(factors(FwP, nfwp) * S);
if diag_opt
l = matfaust.opt_butterfly_faust(l);
end
C = l * r;
if startsWith(dev, 'gpu') if startsWith(dev, 'gpu')
if diag_opt
error('diag_opt on GPU Faust is not yet implemented')
end
C = clone(C, 'dev', 'gpu'); C = clone(C, 'dev', 'gpu');
end end
end end
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment