Mercurial > hg > lots

changeset 0:9c6ab153f1f2

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Initial repository from previous disjointness project
author danieleb@code.soundsoftware.ac.uk
date Tue, 14 Jun 2011 14:44:20 +0100
parents
children 856ac6a53070
files invlappedorthotransform.m lappedorthobasis.m lappedorthoplot.m lappedorthotransform.m lappedwindow.m
diffstat 5 files changed, 423 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/invlappedorthotransform.m	Tue Jun 14 14:44:20 2011 +0100
@@ -0,0 +1,59 @@
+function x = invlappedorthotransform(y,wLength,orthoBasis,tailLength,tailFunc)
+%>INVLAPPEDORTHOTRANSFORM Inverse lapped orthogonal transform (iLOT)
+%>
+%> x = lappedorthotransform(y,wLength,orthoBasis,tailLength,tailFunc)
+%>
+%> returns the inverse lapped orthogonal transform of the signal y 
+%>
+%> INPUT:
+%> - y: vector containing the LOT
+%> - wLength: vector containing the lengths of the overlapping windows (the
+%> sum of wLength must equal the length of x)
+%> - orthobasis: either a string or a function handle corresponding to an
+%> orhogonal transform to be used in each overlapping window. The default
+%> is DCT-IV (which is globally referred as MDCT)
+%> - tailLength: length of the tail of the overlapping windows. Two
+%> consecutive windows overlap in a portion of dimension 2*tailLength. 
+%> The maximum tailLength cannot exceed half the length of the smallest 
+%> window. 
+%> - tailFunc: either a string or a function handle used to define the tails
+%> of the overlapping windows. The default is 'sin2' (twicely differentiable
+%> sinusoidal).
+%>
+%> OUTPUT:
+%> - x: inverse LOT
+%>
+%> REFERENCE:
+%> S. Mallat, A Wavelet Tour of Signal Processing
+%>
+%> SEE ALSO
+%> lappedorthobasis, lappedwindow
+%>
+%> ----------------------------------------------------------- %
+%>  Daniele Barchiesi, daniele.barchiesi@eecs.qmul.ac.uk       %
+%>  Centre for Digital Music, Queen Mary University of London  %
+%>  May 2011                                                  %
+%> ----------------------------------------------------------- %
+%%> Check input & defaults
+error(nargchk(2, 5, nargin, 'struct'));
+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
+if length(x)<sum(wLength), error('you must pad the signal!'); end
+
+if ischar(orthoBasis)
+    switch orthoBasis
+        case 'mdct'
+            orthoFun = @(x) dcti(x,'IV');
+        otherwise
+            error('invalid basis function');
+    end
+elseif isa(orthoBasis,'function_handle')
+    orthoFun = orthoBasis;
+else
+    error('invalid basis function');
+end
+
+%%> Compute transform
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lappedorthobasis.m	Tue Jun 14 14:44:20 2011 +0100
@@ -0,0 +1,100 @@
+function phi = lappedorthobasis(wLength,orthobasis,tailLength,tailFunc)
+%LAPPEDORTHOBASIS Lapped orthogonal basis matrix
+%
+% phi = lappedorthobasis(wLength,orthobasis,tailLength,tailFunc)
+%
+% returns a matrix whose columns contain the atoms of a lapped orhogonal
+% basis function.
+%
+% INPUT:
+% - wLength: vector containing the lengths of the overlapping windows
+% - orthobasis: either a string or a function handle corresponding to an
+% orhogonal transform to be used in each overlapping window. The default
+% is DCT-IV (which is globally referred as MDCT)
+% - tailLength: length of the tail of the overlapping windows. Two
+% consecutive windows overlap in a portion of dimension 2*tailLength. 
+% The maximum tailLength cannot exceed half the length of the smallest 
+% window. 
+% - tailFunc: either a string or a function handle used to define the tails
+% of the overlapping windows. The default is 'sin2' (twicely differentiable
+% sinusoidal).
+%
+% OUTPUT:
+% - phi: matrix containing the lapped orthogonal dictionary.
+%
+% REFERENCE:
+% S. Mallat, A Wavelet Tour of Signal Processing
+%
+% SEE ALSO
+% lappedorthotransform, lappedwindow
+%
+% ----------------------------------------------------------- %
+%  Daniele Barchiesi, daniele.barchiesi@eecs.qmul.ac.uk       %
+%  Centre for Digital Music, Queen Mary University of London  %
+%  Apr. 2011                                                  %
+% ----------------------------------------------------------- %
+%% check inputs & defaults
+error(nargchk(1, 4, nargin, 'struct'));
+if ~exist('tailFunc','var')   || isempty(tailFunc), tailFunc = 'sin2'; end
+if ~exist('tailLength','var') || isempty(tailLength), tailLength = floor(min(wLength)/2); end
+if ~exist('orthobasis','var') || isempty(orthobasis), orthobasis = 'mdct'; end
+
+% select local orthonormal transform
+if ischar(orthobasis)
+    switch orthobasis
+        case 'mdct'
+            basisFun = @(x) dcti(x,'IV');
+        otherwise
+            error('invalid basis function');
+    end
+elseif isa(orthobasis,'function_handle')
+    basisFun = orthobasis;
+else
+    error('invalid basis function');
+end
+
+%% calculate dictionary
+nWindows  = length(wLength);
+sigLength = sum(wLength);                 
+phi       = zeros(sigLength);
+suppLen   = @(p) wLength(p)+2*tailLength;
+
+% first frame
+g = lappedwindow(wLength(1),tailLength,tailFunc,'first');
+orthoMat = basisFun(eye(wLength(1)));
+orthoMatExtended = orthomatextend(orthoMat,tailLength,'first');
+phi(1:wLength(1)+tailLength,1:wLength(1)) = diag(sparse(g))*orthoMatExtended;
+
+% central frames
+for p=2:nWindows-1
+    g = lappedwindow(wLength(p),tailLength,tailFunc);
+    orthoMat = basisFun(eye(wLength(p)));
+    orthoMatExtended = orthomatextend(orthoMat,tailLength);
+    phi(sum(wLength(1:p-1))-tailLength+(1:suppLen(p)),sum(wLength(1:p-1))+(1:wLength(p)))...
+        = diag(sparse(g))*orthoMatExtended;
+end
+
+%last frame
+g = lappedwindow(wLength(end),tailLength,tailFunc,'last');
+if strcmpi(orthobasis,'mdct')
+    orthoMat = dctmatrix(wLength(end),[],'I');
+else
+    orthoMat = basisFun(eye(wLength(end)));
+end
+orthoMatExtended = orthomatextend(orthoMat,tailLength,'last');
+phi(sum(wLength(1:end-1))-tailLength+(1:wLength(end)+tailLength),sum(wLength(1:end-1))+(1:wLength(end)))...
+    = diag(sparse(g))*orthoMatExtended;
+
+function orthomatextended = orthomatextend(orthomat,tailLength,special)
+% Periodic extension of local orthonormal basis
+
+if ~exist('special','var') || isempty(special), special = []; end 
+
+if strcmpi(special,'first')
+    orthomatextended = [orthomat; -orthomat(end:-1:end-tailLength+1,:)];
+elseif strcmpi(special,'last')
+    orthomatextended = [flipud(orthomat(1:tailLength,:)); orthomat];
+else
+    orthomatextended = [flipud(orthomat(1:tailLength,:)); orthomat; ...
+        -orthomat(end:-1:end-tailLength+1,:)];
+end
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lappedorthoplot.m	Tue Jun 14 14:44:20 2011 +0100
@@ -0,0 +1,51 @@
+function varargout = lappedorthoplot(y,wLength,Fs)
+%LAPPEDORTHOPLOT plot of the lapped orthogonal transform
+%
+% varargout = lappedorthoplot(y,wLength,timeVec)
+%
+% plots a spectrogram-like visualisation of the lapped orthogonal
+% transform coefficients contained in y
+%
+% INPUT:
+% - y: vector containing the lapped orthogonal transform coefficients
+% - wLength: vector containing the lengths of the overlapping windows
+% - Fs: sampling frequency (default 44100)
+%
+% OUTPUT:
+% - varargout{1}: optionally, the function outputs the matrix containing 
+% the spectrogram-like representation.
+%
+% SEE ALSO
+% lappedorthotransform
+%
+% ----------------------------------------------------------- %
+%  Daniele Barchiesi, daniele.barchiesi@eecs.qmul.ac.uk       %
+%  Centre for Digital Music, Queen Mary University of London  %
+%  Apr. 2011                                                  %
+% ----------------------------------------------------------- %
+nWindows = length(wLength);
+maxLength = max(wLength);
+sigLength = sum(wLength);
+
+if ~exist('Fs','var') || isempty(Fs), Fs = 44100; end
+
+timeVec = linspace(0,sigLength/Fs,10*ceil(sigLength/min(wLength))); 
+
+S = zeros(maxLength,length(timeVec));
+timeSupp  = @(p) floor(wLength(p)*length(timeVec)/sigLength);
+interpFun = @(p,v) interp1(1:wLength(p),abs(v),linspace(1,wLength(p),maxLength))';
+S(:,1:timeSupp(1)) = repmat(interpFun(1,y(1:wLength(1))),1,timeSupp(1));
+for i=2:nWindows
+    S(:,sum(timeSupp(1:i-1))+(1:timeSupp(i))) = ...
+        repmat(interpFun(i,y(sum(wLength(1:i-1))+(1:wLength(i)))),1,timeSupp(i));
+end
+
+SmagdB = mag2db(S);
+if nargout, varargout{1} = SmagdB; end
+
+surf(timeVec,1:maxLength,SmagdB,'EdgeColor','none');
+axis xy; axis tight;
+colormap(jet);
+view(0,90);
+set(gca,'YScale','log');
+title(['l1/l2 = ' num2str(l1overl2(y))])
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lappedorthotransform.m	Tue Jun 14 14:44:20 2011 +0100
@@ -0,0 +1,117 @@
+function y = lappedorthotransform(x,wLength,orthoBasis,tailLength,tailFunc)
+%LAPPEDORTHOTRANSFORM Lapped orthogonal transform
+%
+% y = lappedorthotransform(x,wLength,orthoBasis,tailLength,tailFunc)
+%
+% returns the lapped orthogonal transform of the signal x 
+%
+% INPUT:
+% - x: vector containing the input signal 
+% - wLength: vector containing the lengths of the overlapping windows (the
+% sum of wLength must equal the length of x)
+% - orthobasis: either a string or a function handle corresponding to an
+% orhogonal transform to be used in each overlapping window. The default
+% is DCT-IV (which is globally referred as MDCT)
+% - tailLength: length of the tail of the overlapping windows. Two
+% consecutive windows overlap in a portion of dimension 2*tailLength. 
+% The maximum tailLength cannot exceed half the length of the smallest 
+% window. 
+% - tailFunc: either a string or a function handle used to define the tails
+% of the overlapping windows. The default is 'sin2' (twicely differentiable
+% sinusoidal).
+%
+% OUTPUT:
+% - y: vector of coefficients of the lapped orthogonal transform
+%
+% REFERENCE:
+% S. Mallat, A Wavelet Tour of Signal Processing
+%
+% SEE ALSO
+% lappedorthobasis, lappedwindow
+%
+% ----------------------------------------------------------- %
+%  Daniele Barchiesi, daniele.barchiesi@eecs.qmul.ac.uk       %
+%  Centre for Digital Music, Queen Mary University of London  %
+%  Apr. 2011                                                  %
+% ----------------------------------------------------------- %
+%% Check input & defaults
+error(nargchk(2, 5, nargin, 'struct'));
+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
+if isscalar(wLength)
+    nWindows = floor(length(x)/wLength);
+    wLength = wLength*ones(nWindows,1);
+end
+if length(x)<sum(wLength), error('you must pad the signal!'); end
+if ischar(orthoBasis)
+    switch orthoBasis
+        case 'mdct'
+            orthoFun = @(x) dcti(x,'IV');
+        otherwise
+            error('invalid basis function');
+    end
+elseif isa(orthoBasis,'function_handle')
+    orthoFun = orthoBasis;
+else
+    error('invalid basis function');
+end
+
+%% Compute transform
+nWindows  = length(wLength);
+sigLength = sum(wLength);
+y         = zeros(sigLength,1);
+suppLen   = @(p) wLength(p)+2*tailLength;
+
+% first frame
+h = lappedprojector(x(1:wLength(1)+tailLength),wLength(1),tailLength,tailFunc,'first');
+y(1:wLength(1)) = orthoFun(h(1:wLength(1)));
+
+% central frames
+for p=2:nWindows-1
+    h = lappedprojector(x(sum(wLength(1:p-1))-tailLength+(1:suppLen(p))),...
+        wLength(p),tailLength,tailFunc);
+    y(sum(wLength(1:p-1))+(1:wLength(p))) = orthoFun(h(tailLength+(1:wLength(p))));
+end
+    
+%last frame
+h = lappedprojector(x(sum(wLength(1:end-1))-tailLength+(1:wLength(end)+tailLength)),...
+    wLength(end),tailLength,tailFunc,'last');
+if strcmpi(orthoBasis,'mdct')
+    y(sum(wLength(1:end-1))+(1:wLength(end))) = dcti(h(tailLength+(1:wLength(end))),'I');
+else
+    y(sum(wLength(1:end-1))+(1:wLength(end))) = orthoFun(h(tailLength+(1:wLength(end))));
+end
+
+function h = lappedprojector(f,wLength,tailLength,tailFunc,special)
+% Projects the input f into the space of functions with support contained
+% in the interval defined by the lapped window defined according to the
+% input parameteres
+
+% check inputs
+error(nargchk(2, 5, nargin, 'struct'));
+if ~exist('special','var') || isempty(special), special = []; end
+if ~exist('tailFunx','var') || isempty(tailFunc), tailFunc = 'sin2'; end
+if ~exist('tailLength','var') || isempty(tailLength)
+    tailLength = floor(min(wLength)/2);
+end
+
+% generate the window and project
+g = lappedwindow(wLength,tailLength,tailFunc,special);
+h = f;
+if ~isempty(special)
+    if strcmpi(special,'first')
+        O = wLength-tailLength+(1:2*tailLength);
+        h(O) = g(O).*f(O) - g(fliplr(O)).*f(fliplr(O));
+    elseif strcmpi(special,'last')
+        O = 1:2*tailLength;
+        h(O) = g(O).*f(O) + g(fliplr(O)).*f(fliplr(O));
+    end
+else
+    Os = 1:2*tailLength;
+    Oe = wLength+(1:2*tailLength);
+    h(Os) = g(Os).*f(Os) + g(fliplr(Os)).*f(fliplr(Os));
+    h(Oe) = g(Oe).*f(Oe) - g(fliplr(Oe)).*f(fliplr(Oe));
+end
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lappedwindow.m	Tue Jun 14 14:44:20 2011 +0100
@@ -0,0 +1,96 @@
+function g = lappedwindow(wLength,tailLength,tailFunc,special)
+%LAPPEDWINDOW Lapped window used for lapped orthogonal transform
+%
+% g = lappedwindow(wLength,tailLength,tailFunc,special)
+%
+% returns the window g used in the lapped orthogonal transform
+%
+% INPUT:
+% - wLength: length of the window
+% - tailLength: length of the tail of the overlapping window. 
+%
+%     ......................
+%   |.                      .|
+%   .                        .
+%  .|                        |.
+% . |---------wLength--------| .
+%                            |-| 
+%                         tailLength
+%
+% The maximum tailLength cannot exceed half the window length. 
+% - tailFunc: either a string or a function handle used to define the tails
+% of the overlapping windows. The default is 'sin2' (twicely differentiable
+% sinusoidal).
+%
+% OUTPUT:
+% - g: vector containing the window
+%
+% REFERENCE:
+% S. Mallat, A Wavelet Tour of Signal Processing
+%
+% SEE ALSO
+% lappedorthobasis, lappedorthotransform
+%
+% ----------------------------------------------------------- %
+%  Daniele Barchiesi, daniele.barchiesi@eecs.qmul.ac.uk       %
+%  Centre for Digital Music, Queen Mary University of London  %
+%  Apr. 2011                                                  %
+% ----------------------------------------------------------- %
+
+%% Check inputs
+error(nargchk(2, 4, nargin, 'struct'));
+if ~exist('special','var') || isempty(special), special = []; end
+if ~exist('tailFunx','var') || isempty(tailFunc), tailFunc = 'sin2'; end
+
+%check that the window is not too short
+if wLength<2*tailLength
+    error('the window length must be at least twice as long as the tail length');
+end
+
+%manage tail function
+if ischar(tailFunc)
+    [fname iter] = parsefunctionname(tailFunc);
+    switch lower(fname)
+        case 'sin'
+            tailFunc = @(x)sin((pi/4)*(x+1));
+        otherwise
+            error('invalid tail function');
+    end
+    if iter
+        for i=1:iter
+            tailFun = @(x)tailFun(sin(pi*x/2));
+        end
+    end
+else %if a function handler, check that it corresponds to a valid tail.
+    t = linspace(-1,1,1000);
+    test = norm((tailFunc(t).^2)+(tailFunc(fliplr(t)).^2)-ones(1,length(t)));
+    if test>1e-8
+        error('invalid tail function');
+    end
+end
+
+%% construct the window
+if ~isempty(special)
+    g = ones(wLength+tailLength,1);
+    g(1:2*tailLength) = tailFunc((-tailLength+0.5:tailLength-0.5)'/tailLength);
+    if strcmpi(special,'first')
+        g = flipud(g);
+    end
+else
+    g = ones(wLength+2*tailLength,1);
+    g(1:2*tailLength) = tailFunc((-tailLength+0.5:tailLength-0.5)'/tailLength);
+    g(wLength+1:wLength+2*tailLength) = ...
+        flipud(tailFunc((-tailLength+0.5:tailLength-0.5)'/tailLength));
+end
+
+function [fname iter] = parsefunctionname(tailFunc)
+nameLen = length(tailFunc);
+for i=1:nameLen
+    if real(str2double(tailFunc(i)))>0
+        fname = tailFunc(1:i-1);
+        iter = uint16(str2double(tailFunc(i:end)));
+        return
+    end
+end
+fname = tailFunc;
+iter = [];
\ No newline at end of file