1function c = block_fwt( f, w, J)
2%-*- texinfo -*-
3%@deftypefn {Function} block_fwt
4%@verbatim
5%BLOCK_FWT FWT func. wrapper for a block processing
6%   Usage: c = block_fwt( f, w, J);
7%
8%   Input parameters:
9%         f     : Input data.
10%         w     : Analysis Wavelet Filterbank.
11%         J     : Number of filterbank iterations.
12%
13%   Output parameters:
14%         c      : Coefficient vector.
15%
16%   c = BLOCK_FWT(f,w,J) accepts suitably extended block of data f*
17%   and produces correct coefficients using the SegDWT algorithm (based on
18%   overlap-save block convolution) with wavelet filters defined by w
19%   and J levels. f is expected to be a column vector or a matrix and
20%   the processing is done column-wise.
21%
22%   Do not call this function directly. The function is called from
23%   BLOCKANA when used with frame type 'fwt' and 'segola' block transform
24%   handling see BLOCKFRAMEACCEL.
25%
26%   Function should be independent of block_interface.
27%
28%
29%   References:
30%     Z. Průša. Segmentwise Discrete Wavelet Transform. PhD thesis, Brno
31%     University of Technology, Brno, 2012.
32%
33%
34%@end verbatim
35%@strong{Url}: @url{http://ltfat.github.io/doc/blockproc/private/block_fwt.html}
36%@seealso{block, block_ifwt}
37%@end deftypefn
38
39% Copyright (C) 2005-2016 Peter L. Soendergaard <peter@sonderport.dk>.
40% This file is part of LTFAT version 2.3.1
41%
42% This program is free software: you can redistribute it and/or modify
43% it under the terms of the GNU General Public License as published by
44% the Free Software Foundation, either version 3 of the License, or
45% (at your option) any later version.
46%
47% This program is distributed in the hope that it will be useful,
48% but WITHOUT ANY WARRANTY; without even the implied warranty of
49% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
50% GNU General Public License for more details.
51%
52% You should have received a copy of the GNU General Public License
53% along with this program.  If not, see <http://www.gnu.org/licenses/>.
54
55if nargin<3
56  error('%s: Too few input parameters.',upper(mfilename));
57end;
58
59% Initialize the wavelet filters structure
60%h = fwtinit(h,'ana');
61
62if any(w.a~=w.a(1))
63   error('%s: Non-equal subsampling factors are not supported.',upper(mfilename));
64end
65
66w = fwtinit(w);
67% Extended block length
68Ls = size(f,1);
69% Low-pass filter length
70m = numel(w.h{1}.h);
71% Low-pass subsampling factor
72a = w.a(1);
73% Extension length
74rred = (a^J-1)/(a-1)*(m-a);
75% Block boundaries
76blocksize=w.a(1)^J;
77% Input signal samples to be processed
78
79% This is effectivelly the "negative" right extension described in chapter
80% 4.1.4 in the reference.
81L=rred+floor((Ls-rred)/blocksize)*blocksize;
82
83levelLen = L;
84filtNo = length(w.h);
85subbNo = (filtNo-1)*J+1;
86Lc = zeros(subbNo,1);
87runPtr = 0;
88for jj=1:J
89   for ff=filtNo:-1:2
90      Lc(end-runPtr) = floor((levelLen-m-1)/w.a(ff));
91      runPtr = runPtr + 1;
92   end
93   levelLen = floor((levelLen-m-1)/w.a(1));
94end
95Lc(1)=levelLen;
96
97%
98%[Lc, L] = fwtclength(Ls,h,J,'valid');
99
100% Crop to the right length
101if(Ls>L)
102   f=postpad(f,L);
103end
104
105if Ls<rred+a^J
106   error('%s: Insufficient input signal length for the %s flag. Minimum is %i.',upper(mfilename),'''valid''',rred+a^J);
107end
108
109c = comp_fwt(f,w.h,w.a,J,'valid');
110
111% Do the cropping
112runPtr = 0;
113for jj=1:J-1
114   for ff=filtNo:-1:2
115      cstart = (a^(J-jj)-1)/(a-1)*(m-a);
116      c{end-runPtr} = c{end-runPtr}(cstart+1:end,:);
117      runPtr = runPtr + 1;
118   end
119end
120
121% To the pack format
122c = cell2mat(c);
123
124