annotate lot.m @ 17:0da5eb32e49d tip

minor changes
author Daniele Barchiesi <daniele.barchiesi@eecs.qmul.ac.uk>
date Thu, 10 Nov 2011 15:53:06 +0000
parents 88af68016e8a
children
rev   line source
danieleb@4 1 function y = lot(x,wLength,orthoBasis,tailLength,tailFunc)
danieleb@15 2 % LOT lapped orthogonal transform
danieleb@0 3 %
danieleb@14 4 % Example
danieleb@0 5 % y = lappedorthotransform(x,wLength,orthoBasis,tailLength,tailFunc)
danieleb@0 6 %
danieleb@14 7 % Input
danieleb@14 8 % - x: vector or matrix containing the input signal. If matrix, the function
danieleb@4 9 % acts columnwise.
danieleb@14 10 % - wLength: scalar or vector containing the lengths of the overlapping
danieleb@4 11 % windows (if vector, the sum of wLength must equal the length of x)
danieleb@14 12 % - orthobasis: (optional, default='mdct') either a string or a function
danieleb@4 13 % handle corresponding to an orhogonal transform to be used in each
danieleb@4 14 % overlapping window.
danieleb@14 15 % - tailLength: (optional, default=floor(min(wLength)/2)) length of the tail
danieleb@4 16 % of the overlapping windows. Two consecutive windows overlap in a portion
danieleb@4 17 % of dimension 2*tailLength. The maximum tailLength cannot exceed half the
danieleb@4 18 % length of the smallest window.
danieleb@14 19 % - tailFunc: (optional, default='sin2') either a string or a function handle
danieleb@4 20 % used to define the tails of the overlapping windows (see LAPPEDWINDOW for
danieleb@4 21 % more details).
danieleb@0 22 %
danieleb@14 23 % Output
danieleb@14 24 % - y: vector or matrix of coefficients of the lapped orthogonal transform.
danieleb@0 25 %
danieleb@14 26 % References:
danieleb@0 27 % S. Mallat, A Wavelet Tour of Signal Processing
danieleb@0 28 %
danieleb@14 29 % See also LAPPEDWINDOW, ilot
danieleb@4 30
danieleb@4 31 % Author(s): Daniele Barchiesi
danieleb@14 32 % Copyright 2011-2011
danieleb@4 33
danieleb@0 34 %% Check input & defaults
danieleb@0 35 error(nargchk(2, 5, nargin, 'struct'));
danieleb@4 36
danieleb@4 37 if isscalar(wLength) %create vector of fixed window lengths
danieleb@4 38 nWindows = floor(length(x)/wLength);
danieleb@4 39 wLength = wLength*ones(nWindows,1);
danieleb@4 40 end
danieleb@4 41
danieleb@0 42 if ~exist('orthoBasis','var') || isempty(orthoBasis), orthoBasis = 'mdct'; end
danieleb@0 43 if ~exist('tailFunc','var') || isempty(tailFunc), tailFunc = 'sin2'; end
danieleb@0 44 if ~exist('tailLength','var') || isempty(tailLength)
danieleb@0 45 tailLength = floor(min(wLength)/2);
danieleb@0 46 end
danieleb@4 47
danieleb@4 48 % check wLength
danieleb@0 49 if length(x)<sum(wLength), error('you must pad the signal!'); end
danieleb@4 50 % define function handle to the forward local transform
danieleb@8 51 orthoFun = orthohandle(orthoBasis,'for');
danieleb@0 52
danieleb@4 53 %% Matrix case: act columnwise
danieleb@4 54 if size(x,2)>1
danieleb@4 55 y = zeros(size(x));
danieleb@4 56 for iCol=1:size(x,2)
danieleb@6 57 y(:,iCol) = lot(x(:,iCol),wLength,orthoBasis,...
danieleb@4 58 tailLength,tailFunc);
danieleb@4 59 end
danieleb@4 60 return
danieleb@4 61 end
danieleb@4 62
danieleb@0 63 %% Compute transform
danieleb@0 64 nWindows = length(wLength);
danieleb@0 65 sigLength = sum(wLength);
danieleb@0 66 y = zeros(sigLength,1);
danieleb@4 67 % length of the support of the p-th window
danieleb@0 68 suppLen = @(p) wLength(p)+2*tailLength;
danieleb@0 69
danieleb@0 70 % first frame
danieleb@4 71 % project input signal onto orthogonal subspace associated with the first
danieleb@4 72 % window
danieleb@0 73 h = lappedprojector(x(1:wLength(1)+tailLength),wLength(1),tailLength,tailFunc,'first');
danieleb@4 74 % transform projected signal
danieleb@0 75 y(1:wLength(1)) = orthoFun(h(1:wLength(1)));
danieleb@0 76
danieleb@0 77 % central frames
danieleb@0 78 for p=2:nWindows-1
danieleb@4 79 % project input signal onto orthogonal subspace associated with the
danieleb@4 80 % p-th window
danieleb@0 81 h = lappedprojector(x(sum(wLength(1:p-1))-tailLength+(1:suppLen(p))),...
danieleb@0 82 wLength(p),tailLength,tailFunc);
danieleb@4 83 % transform projected signal
danieleb@0 84 y(sum(wLength(1:p-1))+(1:wLength(p))) = orthoFun(h(tailLength+(1:wLength(p))));
danieleb@0 85 end
danieleb@0 86
danieleb@0 87 %last frame
danieleb@4 88 % project input signal onto orthogonal subspace associated with the
danieleb@4 89 % last window
danieleb@0 90 h = lappedprojector(x(sum(wLength(1:end-1))-tailLength+(1:wLength(end)+tailLength)),...
danieleb@0 91 wLength(end),tailLength,tailFunc,'last');
danieleb@4 92 % transform projected signal
danieleb@4 93 if strcmpi(orthoBasis,'mdct') %id mdct use DCT I at last frame to avoid artefacts
danieleb@0 94 y(sum(wLength(1:end-1))+(1:wLength(end))) = dcti(h(tailLength+(1:wLength(end))),'I');
danieleb@0 95 else
danieleb@0 96 y(sum(wLength(1:end-1))+(1:wLength(end))) = orthoFun(h(tailLength+(1:wLength(end))));
danieleb@0 97 end
danieleb@0 98
danieleb@0 99 function h = lappedprojector(f,wLength,tailLength,tailFunc,special)
danieleb@0 100 % Projects the input f into the space of functions with support contained
danieleb@4 101 % in the interval defined by the lapped window.
danieleb@0 102
danieleb@0 103 % check inputs
danieleb@0 104 error(nargchk(2, 5, nargin, 'struct'));
danieleb@0 105 if ~exist('special','var') || isempty(special), special = []; end
danieleb@4 106 if ~exist('tailFunc','var') || isempty(tailFunc), tailFunc = 'sin2'; end
danieleb@0 107 if ~exist('tailLength','var') || isempty(tailLength)
danieleb@0 108 tailLength = floor(min(wLength)/2);
danieleb@0 109 end
danieleb@0 110
danieleb@0 111 % generate the window and project
danieleb@0 112 g = lappedwindow(wLength,tailLength,tailFunc,special);
danieleb@0 113 h = f;
danieleb@0 114 if ~isempty(special)
danieleb@0 115 if strcmpi(special,'first')
danieleb@4 116 % overlapping region
danieleb@0 117 O = wLength-tailLength+(1:2*tailLength);
danieleb@4 118 % project
danieleb@0 119 h(O) = g(O).*f(O) - g(fliplr(O)).*f(fliplr(O));
danieleb@0 120 elseif strcmpi(special,'last')
danieleb@4 121 % overlapping portion
danieleb@0 122 O = 1:2*tailLength;
danieleb@4 123 % project
danieleb@0 124 h(O) = g(O).*f(O) + g(fliplr(O)).*f(fliplr(O));
danieleb@0 125 end
danieleb@0 126 else
danieleb@4 127 % start overlapping region
danieleb@0 128 Os = 1:2*tailLength;
danieleb@4 129 % end overlapping region
danieleb@0 130 Oe = wLength+(1:2*tailLength);
danieleb@4 131 % project
danieleb@0 132 h(Os) = g(Os).*f(Os) + g(fliplr(Os)).*f(fliplr(Os));
danieleb@0 133 h(Oe) = g(Oe).*f(Oe) - g(fliplr(Oe)).*f(fliplr(Oe));
danieleb@0 134 end