annotate ilot.m @ 4:e23a23349e31

Added inverse transform
author danieleb@code.soundsoftware.ac.uk
date Thu, 16 Jun 2011 13:16:29 +0100
parents invlappedorthotransform.m@9c6ab153f1f2
children 084000643ad1
rev   line source
danieleb@4 1 function x = ilot(y,wLength,orthoBasis,tailLength,tailFunc)
danieleb@4 2 %INVLAPPEDORTHOTRANSFORM Inverse lapped orthogonal transform
danieleb@4 3 %
danieleb@4 4 % USAGE
danieleb@4 5 % x = lappedorthotransform(y,wLength,orthoBasis,tailLength,tailFunc)
danieleb@4 6 %
danieleb@4 7 % INPUT:
danieleb@4 8 % y: vector or matrix containing the LOT coefficients. If matrix,
danieleb@4 9 % the function acts columnwise.
danieleb@4 10 %
danieleb@4 11 % wLength: scalar or vector containing the lengths of the overlapping
danieleb@4 12 % windows (if vector, the sum of wLength must equal the length of y)
danieleb@4 13 %
danieleb@4 14 % orthobasis: (optional, default='mdct') either a string or a function
danieleb@4 15 % handle corresponding to an orhogonal transform to be used in each
danieleb@4 16 % overlapping window.
danieleb@4 17 %
danieleb@4 18 % tailLength: (optional, default=floor(min(wLength)/2)) length of the tail
danieleb@4 19 % of the overlapping windows. Two consecutive windows overlap in a interval
danieleb@4 20 % of dimension 2*tailLength. The maximum tailLength cannot exceed half the
danieleb@4 21 % length of the smallest window.
danieleb@4 22 %
danieleb@4 23 % tailFunc: (optional, default='sin2') either a string or a function handle
danieleb@4 24 % used to define the tails of the overlapping windows (see LAPPEDWINDOW for
danieleb@4 25 % more details).
danieleb@4 26 %
danieleb@4 27 % OUTPUT:
danieleb@4 28 % x: vector or matrix of coefficients of the inverse LOT.
danieleb@4 29 %
danieleb@4 30 % REFERENCES:
danieleb@4 31 % S. Mallat, A Wavelet Tour of Signal Processing
danieleb@4 32 %
danieleb@4 33 % SEE ALSO
danieleb@4 34 % LAPPEDORTHOBASIS, LAPPEDWINDOW, LAPPEDORTHOTRANSFORM
danieleb@4 35
danieleb@4 36 % Author(s): Daniele Barchiesi
danieleb@4 37 % Copyright QMUL
danieleb@4 38 % $Revision: 1
danieleb@4 39
danieleb@4 40 %% Check input & defaults
danieleb@0 41 error(nargchk(2, 5, nargin, 'struct'));
danieleb@4 42
danieleb@4 43 if isscalar(wLength) %create vector of fixed window lengths
danieleb@4 44 nWindows = floor(length(y)/wLength);
danieleb@4 45 wLength = wLength*ones(nWindows,1);
danieleb@4 46 end
danieleb@4 47
danieleb@0 48 if ~exist('orthoBasis','var') || isempty(orthoBasis), orthoBasis = 'mdct'; end
danieleb@0 49 if ~exist('tailFunc','var') || isempty(tailFunc), tailFunc = 'sin2'; end
danieleb@0 50 if ~exist('tailLength','var') || isempty(tailLength)
danieleb@4 51 tailLength = floor(min(wLength)/2);
danieleb@0 52 end
danieleb@0 53
danieleb@4 54 % check wLength
danieleb@4 55 if length(y)<sum(wLength), error('you must pad the signal!'); end
danieleb@4 56 % define function handle to the forward local transform
danieleb@0 57 if ischar(orthoBasis)
danieleb@0 58 switch orthoBasis
danieleb@0 59 case 'mdct'
danieleb@4 60 invOrthoFun = @(y) idcti(y,'IV');
danieleb@0 61 otherwise
danieleb@0 62 error('invalid basis function');
danieleb@0 63 end
danieleb@0 64 elseif isa(orthoBasis,'function_handle')
danieleb@4 65 invOrthoFun = orthoBasis;
danieleb@0 66 else
danieleb@0 67 error('invalid basis function');
danieleb@0 68 end
danieleb@0 69
danieleb@4 70 %% Matrix case: act columnwise
danieleb@4 71 if size(y,2)>1
danieleb@4 72 x = zeros(size(y));
danieleb@4 73 for iCol=1:size(y,2)
danieleb@4 74 x(:,iCol) = invlappedorthotransform(y(:,iCol),wLength,orthoBasis,...
danieleb@4 75 tailLength,tailFunc);
danieleb@4 76 end
danieleb@4 77 return
danieleb@4 78 end
danieleb@4 79
danieleb@4 80 %% Compute transform
danieleb@4 81 nWindows = length(wLength);
danieleb@4 82 sigLength = sum(wLength);
danieleb@4 83 x = zeros(sigLength,1);
danieleb@4 84
danieleb@4 85 % first frame
danieleb@4 86 g = lappedwindow(wLength(1),tailLength,tailFunc,'first');
danieleb@4 87 h = invOrthoFun(y(1:wLength(1))); % folded ilot
danieleb@4 88 C = 1:wLength(1)-tailLength; % central interval
danieleb@4 89 Om = wLength(1)-tailLength+(1:tailLength); %overlapping interval (end)
danieleb@4 90 hp = invOrthoFun(y(wLength(1)+(1:wLength(2)))); % next folded ilot
danieleb@4 91 x(C) = h(C);
danieleb@4 92 % unfold
danieleb@4 93 x(Om) = g(Om).*h(Om) + g(end:-1:end-tailLength+1).*hp(tailLength:-1:1);
danieleb@4 94
danieleb@4 95 % central frames
danieleb@4 96 for p=2:nWindows-1
danieleb@4 97 hm = h; % previous folded ilot
danieleb@4 98 h = hp; % current folded ilot
danieleb@4 99 % next folded ilot
danieleb@4 100 if strcmpi(orthoBasis,'mdct') && p==nWindows-1 % if mdct use idct I
danieleb@4 101 hp = idcti(y(sum(wLength(1:end-1))+(1:wLength(end))),'I');
danieleb@4 102 else
danieleb@4 103 hp = invOrthoFun(y(sum(wLength(1:end-1))+(1:wLength(end))));
danieleb@4 104 end
danieleb@4 105 g = lappedwindow(wLength(p),tailLength,tailFunc);
danieleb@4 106 % overlapping interval (start)
danieleb@4 107 Op = sum(wLength(1:p-1))+(1:tailLength);
danieleb@4 108 % central interval
danieleb@4 109 C = sum(wLength(1:p-1))+(tailLength+1:wLength(p)-tailLength);
danieleb@4 110 % overlapping interval (end)
danieleb@4 111 Om = sum(wLength(1:p-1))+ wLength(p)-tailLength+(1:tailLength);
danieleb@4 112 % unfold
danieleb@4 113 x(Op) = g(tailLength+(1:tailLength)).*h(1:tailLength) - ...
danieleb@4 114 g(tailLength:-1:1).*hm(end:-1:end-tailLength+1);
danieleb@4 115 x(C) = h(tailLength+1:wLength(p)-tailLength);
danieleb@4 116 x(Om) = g(wLength(p)+(1:tailLength)).*h(wLength(p)-tailLength+(1:tailLength)) + ...
danieleb@4 117 g(end:-1:end-tailLength+1).*hp(tailLength:-1:1);
danieleb@4 118 end
danieleb@4 119
danieleb@4 120 %last frame
danieleb@4 121 hm = h; % previous folded ilot
danieleb@4 122 h = hp; % current folded ilot
danieleb@4 123 g = lappedwindow(wLength(end),tailLength,tailFunc);
danieleb@4 124 % overlapping interval (start)
danieleb@4 125 Op = sum(wLength(1:end-1))+(1:tailLength);
danieleb@4 126 % central interval
danieleb@4 127 C = sum(wLength(1:end-1))+(tailLength+1:wLength(end));
danieleb@4 128 % unfold
danieleb@4 129 x(Op) = g(tailLength+(1:tailLength)).*h(1:tailLength) - ...
danieleb@4 130 g(tailLength:-1:1).*hm(end:-1:end-tailLength+1);
danieleb@4 131 x(C) = h(tailLength+1:end);