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