1function c=comp_nonsepdgt_shear(f,g,a,M,s0,s1,br); 2%-*- texinfo -*- 3%@deftypefn {Function} comp_nonsepdgt_shear 4%@verbatim 5%COMP_NONSEPDGT_SHEAR Compute Non-separable Discrete Gabor transform 6% Usage: c=nonsepdgt(f,g,a,M,lt); 7% 8% This is a computational subroutine, do not call it directly. 9%@end verbatim 10%@strong{Url}: @url{http://ltfat.github.io/doc/comp/comp_nonsepdgt_shear.html} 11%@end deftypefn 12 13% Copyright (C) 2005-2016 Peter L. Soendergaard <peter@sonderport.dk>. 14% This file is part of LTFAT version 2.3.1 15% 16% This program is free software: you can redistribute it and/or modify 17% it under the terms of the GNU General Public License as published by 18% the Free Software Foundation, either version 3 of the License, or 19% (at your option) any later version. 20% 21% This program is distributed in the hope that it will be useful, 22% but WITHOUT ANY WARRANTY; without even the implied warranty of 23% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 24% GNU General Public License for more details. 25% 26% You should have received a copy of the GNU General Public License 27% along with this program. If not, see <http://www.gnu.org/licenses/>. 28 29% AUTHOR : Nicki Holighaus and Peter L. Soendergaard 30% TESTING: TEST_NONSEPDGT 31% REFERENCE: REF_NONSEPDGT 32 33% Assert correct input. 34 35L=size(f,1); 36W=size(f,2); 37b=L/M; 38N=L/a; 39 40c=zeros(M,N,W,assert_classname(f,g)); 41 42ar = a*b/br; 43Mr = L/br; 44Nr = L/ar; 45 46 47if s1 ~= 0 48 p = comp_pchirp(L,s1); 49 g = p.*g; 50 f = bsxfun(@times,p,f); 51end 52 53if s0 == 0 54 55 c_rect = comp_dgt_long(f,g,a,M); 56 tmp1=mod(s1*a*(L+1),2*N); 57 58 for k=0:Nr-1 59 phsidx= mod(mod(tmp1*k,2*N)*k,2*N); 60 phs = exp(pi*1i*phsidx/N); 61 idx2 = floor(mod(-s1*k*a+(0:Mr-1)*b,L)/b); 62 63 c(idx2+1,k+1,:) = c_rect(:,k+1,:).*phs; 64 end; 65 66else 67 twoN=2*N; 68 cc1=ar/a; 69 cc2=mod(-s0*br/a,twoN); 70 cc3=mod(a*s1*(L+1),twoN); 71 cc4=mod(cc2*br*(L+1),twoN); 72 cc5=mod(2*cc1*br,twoN); 73 cc6=mod((s0*s1+1)*br,L); 74 75 p = comp_pchirp(L,-s0); 76 g = p.*fft(g)/L; 77 f = bsxfun(@times,p,fft(f)); 78 79 c_rect = comp_dgt_long(f,g,br,Nr); 80 81 for k=0:Nr-1 82 for m=0:Mr-1 83 sq1=mod(k*cc1+cc2*m,twoN); 84 phsidx = mod(mod(cc3*sq1.^2,twoN)-mod(m*(cc4*m+k*cc5),twoN),twoN); 85 phs = exp(pi*1i*phsidx/N); 86 87 idx1 = mod( cc1*k+cc2*m,N); 88 idx2 = floor(mod(-s1*ar*k+cc6*m,L)/b); 89 90 c(idx2+1,idx1+1,:) = c_rect(mod(-k,Nr)+1,m+1,:).*phs; 91 end; 92 end; 93end; 94 95 96% This is some old code just kept around for historical/documentational 97% purposes 98if 0 99 100 ind = [ar 0; 0 br]*[kron((0:L/ar-1),ones(1,L/br));kron(ones(1,L/ar), ... 101 (0:L/br-1))]; 102 103 104 phs = reshape(mod((s1*(ind(1,:)-s0*ind(2,:)).^2+s0*ind(2,:).^2)*(L+1) ... 105 -2*(s0 ~= 0)*ind(1,:).*ind(2,:),2*L),Mr,Nr); 106 107 ind_final = [1 0;-s1 1]*[1 -s0;0 1]*ind; 108 109 phs = exp(pi*1i*phs/L); 110 111 ind_final = mod(ind_final,L); 112 113 if s1 ~= 0 114 g = comp_pchirp(L,s1).*g; 115 f = bsxfun(@times,comp_pchirp(L,s1),f); 116 end 117 118 if s0 ~= 0 119 g = comp_pchirp(L,-s0).*fft(g)/L; 120 f = bsxfun(@times,comp_pchirp(L,-s0),fft(f)); 121 122 c_rect = comp_dgt_long(f,g,br,Nr); 123 124 for w=0:W-1 125 126 c(floor(ind_final(2,:)/b)+1+(ind_final(1,:)/a)*M+w*M*N) = ... 127 c_rect(ind(1,[1:Mr,end:-1:Mr+1])/ar+1+(ind(2,:)/br)*Nr+w*M*N).* ... 128 phs(ind(2,:)/br+1+(ind(1,:)/ar)*Mr); 129 end; 130 else 131 132 c_rect = comp_dgt_long(f,g,ar,Mr); 133 134 for w=1:W 135 c_rect(:,:,w) = phs.*c_rect(:,:,w); 136 end; 137 138 % The code line below this comment executes the commented for-loop 139 % using Fortran indexing. 140 % 141 % for jj = 1:size(ind,2) 142 % c(floor(ind_final(2,jj)/b)+1, ind_final(1,jj)/a+1) = ... 143 % c_rect(ind(2,jj)/br+1, ind(1,jj)/ar+1); 144 % end 145 for w=0:W-1 146 c(floor(ind_final(2,:)/b)+1+(ind_final(1,:)/a)*M+w*M*N) = ... 147 c_rect(ind(2,:)/br+1+(ind(1,:)/ar)*Mr+w*M*N); 148 end; 149 end; 150 151end; 152 153 154 155