[39b069]: gabor / tconv.m  Maximize  Restore  History

Download this file

88 lines (69 with data), 2.3 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
function h=tconv(f,g)
%TCONV Twisted convolution
% Usage: h=tconv(f,g);
%
% `tconv(f,g)` computes the twisted convolution of the square matrices
% *f* and *g*.
%
% Let `h=tconv(f,g)` for *f,g* being $L \times L$ matrices. Then *h* is given by
%
% .. L-1 L-1
% h(m+1,n+1) = sum sum f(k+1,l+1)*g(m-k+1,n-l+1)*exp(-2*pi*i*(m-k)*l/L);
% l=0 k=0
%
% .. math:: h\left(m+1,n+1\right)=\sum_{l=0}^{L-1}\sum_{k=0}^{L-1}f\left(k+1,l+1\right)g\left(m-k+1,n-l+1\right)e^{-2\pi i(m-k)l/L}
%
% where $m-k$ and $n-l$ are computed modulo *L*.
%
% If both *f* and *g* are of class `sparse` then *h* will also be a sparse
% matrix. The number of non-zero elements of *h* is usually much larger than
% the numbers for *f* and *g*. Unless *f* and *g* are very sparse, it can be
% faster to convert them to full matrices before calling `tconv`.
%
% The routine |spreadinv|_ can be used to calculate an inverse convolution.
% Define *h* and *r* by::
%
% h=tconv(f,g);
% r=tconv(spreadinv(f),h);
%
% then *r* is equal to *g*.
%
% See also: spreadop, spreadfun, spreadinv
% AUTHOR: Peter Soendergaard
% TESTING: TEST_SPREAD
% REFERENCE: REF_TCONV
error(nargchk(2,2,nargin));
if any(size(f)~=size(g))
error('Input matrices must be same size.');
end;
if size(f,1)~=size(f,2)
error('Input matrices must be square.');
end;
L=size(f,1);
if issparse(f) && issparse(g)
% Version for sparse matrices.
% precompute the Lth roots of unity
% Optimization note : the special properties and symmetries of the
% roots of unity could be exploited to reduce this computation.
% Furthermore here we precompute every possible root if some are
% unneeded.
temp=exp((-i*2*pi/L)*(0:L-1)');
[rowf,colf,valf]=find(f);
[rowg,colg,valg]=find(g);
h=sparse(L,L);
for indf=1:length(valf)
for indg=1:length(valg)
m=mod(rowf(indf)+rowg(indg)-2, L);
n=mod(colf(indf)+colg(indg)-2, L);
h(m+1,n+1)=h(m+1,n+1)+valf(indf)*valg(indg)*temp(mod((m-(rowf(indf)-1))*(colf(indf)-1),L)+1);
end
end
else
% The conversion to 'full' is in order for Matlab to work.
f=ifft(full(f))*L;
g=ifft(full(g))*L;
Tf=comp_col2diag(f);
Tg=comp_col2diag(g);
Th=Tf*Tg;
h=spreadfun(Th);
end;

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

Sign up for the SourceForge newsletter:

JavaScript is required for this form.





No, thanks