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

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; ```