danieleb@4: function x = ilot(y,wLength,orthoBasis,tailLength,tailFunc)
danieleb@4: %INVLAPPEDORTHOTRANSFORM Inverse lapped orthogonal transform
danieleb@4: %
danieleb@4: % USAGE
danieleb@4: % x = lappedorthotransform(y,wLength,orthoBasis,tailLength,tailFunc)
danieleb@4: %
danieleb@4: % INPUT:
danieleb@4: % y: vector or matrix containing the LOT coefficients. If matrix,
danieleb@4: % the function acts columnwise.
danieleb@4: %
danieleb@4: % wLength: scalar or vector containing the lengths of the overlapping
danieleb@4: % windows (if vector, the sum of wLength must equal the length of y)
danieleb@4: %
danieleb@4: % orthobasis: (optional, default='mdct') either a string or a function
danieleb@4: % handle corresponding to an orhogonal transform to be used in each
danieleb@4: % overlapping window.
danieleb@4: %
danieleb@4: % tailLength: (optional, default=floor(min(wLength)/2)) length of the tail
danieleb@4: % of the overlapping windows. Two consecutive windows overlap in a interval
danieleb@4: % of dimension 2*tailLength. The maximum tailLength cannot exceed half the
danieleb@4: % length of the smallest window.
danieleb@4: %
danieleb@4: % tailFunc: (optional, default='sin2') either a string or a function handle
danieleb@4: % used to define the tails of the overlapping windows (see LAPPEDWINDOW for
danieleb@4: % more details).
danieleb@4: %
danieleb@4: % OUTPUT:
danieleb@4: % x: vector or matrix of coefficients of the inverse LOT.
danieleb@4: %
danieleb@4: % REFERENCES:
danieleb@4: % S. Mallat, A Wavelet Tour of Signal Processing
danieleb@4: %
danieleb@4: % SEE ALSO
danieleb@4: % LAPPEDORTHOBASIS, LAPPEDWINDOW, LAPPEDORTHOTRANSFORM
danieleb@4: 
danieleb@4: % Author(s): Daniele Barchiesi
danieleb@4: % Copyright QMUL
danieleb@4: % $Revision: 1
danieleb@4: 
danieleb@4: %% Check input & defaults
danieleb@0: error(nargchk(2, 5, nargin, 'struct'));
danieleb@4: 
danieleb@4: if isscalar(wLength) %create vector of fixed window lengths
danieleb@4:     nWindows = floor(length(y)/wLength);
danieleb@4:     wLength = wLength*ones(nWindows,1);
danieleb@4: end
danieleb@4: 
danieleb@0: if ~exist('orthoBasis','var') || isempty(orthoBasis), orthoBasis = 'mdct'; end
danieleb@0: if ~exist('tailFunc','var')   || isempty(tailFunc), tailFunc = 'sin2'; end
danieleb@0: if ~exist('tailLength','var') || isempty(tailLength)
danieleb@4:     tailLength = floor(min(wLength)/2);
danieleb@0: end
danieleb@0: 
danieleb@4: % check wLength
danieleb@4: if length(y)<sum(wLength), error('you must pad the signal!'); end
danieleb@4: % define function handle to the forward local transform
danieleb@0: if ischar(orthoBasis)
danieleb@0:     switch orthoBasis
danieleb@0:         case 'mdct'
danieleb@4:             invOrthoFun = @(y) idcti(y,'IV');
danieleb@0:         otherwise
danieleb@0:             error('invalid basis function');
danieleb@0:     end
danieleb@0: elseif isa(orthoBasis,'function_handle')
danieleb@4:     invOrthoFun = orthoBasis;
danieleb@0: else
danieleb@0:     error('invalid basis function');
danieleb@0: end
danieleb@0: 
danieleb@4: %% Matrix case: act columnwise
danieleb@4: if size(y,2)>1
danieleb@4:     x = zeros(size(y));
danieleb@4:     for iCol=1:size(y,2)
danieleb@4:         x(:,iCol) = invlappedorthotransform(y(:,iCol),wLength,orthoBasis,...
danieleb@4:             tailLength,tailFunc);
danieleb@4:     end
danieleb@4:     return
danieleb@4: end
danieleb@4: 
danieleb@4: %% Compute transform
danieleb@4: nWindows  = length(wLength);
danieleb@4: sigLength = sum(wLength);
danieleb@4: x         = zeros(sigLength,1);
danieleb@5: p = 1;
danieleb@4: 
danieleb@4: % first frame
danieleb@4: g = lappedwindow(wLength(1),tailLength,tailFunc,'first');
danieleb@4: h = invOrthoFun(y(1:wLength(1))); % folded ilot
danieleb@4: C = 1:wLength(1)-tailLength;      % central interval
danieleb@4: Om = wLength(1)-tailLength+(1:tailLength);  %overlapping interval (end)
danieleb@5: if strcmpi(orthoBasis,'mdct') && p==nWindows-1 % if mdct use idct I
danieleb@5:     hp = idcti(y(sum(wLength(1:end-1))+(1:wLength(end))),'I');
danieleb@5: else
danieleb@5:     hp = invOrthoFun(y(sum(wLength(1:end-1))+(1:wLength(end))));
danieleb@5: end
danieleb@4: x(C) = h(C);
danieleb@4: % unfold
danieleb@4: x(Om) = g(Om).*h(Om) + g(end:-1:end-tailLength+1).*hp(tailLength:-1:1);
danieleb@4: 
danieleb@4: % central frames
danieleb@4: for p=2:nWindows-1
danieleb@4:     hm = h; % previous folded ilot
danieleb@4:     h = hp; % current folded ilot
danieleb@4:     % next folded ilot
danieleb@4:     if strcmpi(orthoBasis,'mdct') && p==nWindows-1 % if mdct use idct I
danieleb@4:         hp = idcti(y(sum(wLength(1:end-1))+(1:wLength(end))),'I');
danieleb@4:     else
danieleb@4:         hp = invOrthoFun(y(sum(wLength(1:end-1))+(1:wLength(end))));
danieleb@4:     end
danieleb@4:     g = lappedwindow(wLength(p),tailLength,tailFunc);
danieleb@4:     % overlapping interval (start)
danieleb@4:     Op = sum(wLength(1:p-1))+(1:tailLength);
danieleb@4:     % central interval
danieleb@4:     C  = sum(wLength(1:p-1))+(tailLength+1:wLength(p)-tailLength);
danieleb@4:     % overlapping interval (end)
danieleb@4:     Om = sum(wLength(1:p-1))+ wLength(p)-tailLength+(1:tailLength);
danieleb@4:     % unfold
danieleb@4:     x(Op) = g(tailLength+(1:tailLength)).*h(1:tailLength) - ...
danieleb@4:         g(tailLength:-1:1).*hm(end:-1:end-tailLength+1);
danieleb@4:     x(C) = h(tailLength+1:wLength(p)-tailLength);
danieleb@4:     x(Om) = g(wLength(p)+(1:tailLength)).*h(wLength(p)-tailLength+(1:tailLength)) + ...
danieleb@4:         g(end:-1:end-tailLength+1).*hp(tailLength:-1:1);
danieleb@4: end
danieleb@4: 
danieleb@4: %last frame
danieleb@4: hm = h; % previous folded ilot
danieleb@4: h = hp; % current folded ilot
danieleb@4: g = lappedwindow(wLength(end),tailLength,tailFunc);
danieleb@4: % overlapping interval (start)
danieleb@4: Op = sum(wLength(1:end-1))+(1:tailLength);
danieleb@4: % central interval
danieleb@4: C  = sum(wLength(1:end-1))+(tailLength+1:wLength(end));
danieleb@4: % unfold
danieleb@4: x(Op) = g(tailLength+(1:tailLength)).*h(1:tailLength) - ...
danieleb@5:     g(tailLength:-1:1).*hm(end:-1:end-tailLength+1);
danieleb@4: x(C) = h(tailLength+1:end);