diff ilot.m @ 4:e23a23349e31

Added inverse transform

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