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 |