# HG changeset patch
# User Daniele Barchiesi <daniele.barchiesi@eecs.qmul.ac.uk>
# Date 1314784390 -3600
# Node ID e3035d45d0145a03239a9c8efb15ad3f633312c3
# Parent  23763c5fbda5672cd172c1f30f535d3b78c21c68
Added support classes

diff -r 23763c5fbda5 -r e3035d45d014 util/classes/@audio/audio.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/classes/@audio/audio.m	Wed Aug 31 10:53:10 2011 +0100
@@ -0,0 +1,54 @@
+classdef audio
+    %% Audio object
+    properties
+        s  %vector containing the audio signal
+        fs %sampling frequency 
+        nBits % number of bits per sample
+        name % string containing the name of the audio file
+        format % string containing the format of the audio file
+    end
+    
+    methods
+        %% Constructor
+        function obj = audio(varargin)
+            error(nargchk(1,5,nargin));
+            if ischar(varargin{1})
+                [~, obj.name obj.format] = fileparts(varargin{1});
+                switch obj.format
+                    case '.wav'
+                        [obj.s obj.fs obj.nBits] = wavread(varargin{1});
+                    otherwise
+                        error('Unsupported audio format')
+                end
+            else
+                obj.s = varargin{1};
+                if nargin>1, obj.fs = varargin{2}; else obj.fs = []; end
+                if nargin>2, obj.nBits = varargin{3}; else obj.nBits = []; end
+                if nargin>3, obj.name = varargin{4}; else obj.name = []; end
+                if nargin>4, obj.format = varargin{5}; else obj.format = []; end
+            end
+        end
+        
+        %% Playback functions
+        function player = play(obj, player)
+            if ~exist('player','var') || isempty(player)
+                player = audioplayer(obj.s,obj.fs);
+            end
+            play(player);
+        end
+        
+        function player = stop(obj, player)
+            if ~exist('player','var') || isempty(player)
+                player = audioplayer(obj.s,obj.fs);
+            end
+            stop(player)
+        end
+        
+        function player = pause(obj, player)
+            if ~exist('player','var') || isempty(player)
+                player = audioplayer(obj.s,obj.fs);
+            end
+            pause(player)
+        end
+    end
+end
diff -r 23763c5fbda5 -r e3035d45d014 util/classes/@audio/plot.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/classes/@audio/plot.m	Wed Aug 31 10:53:10 2011 +0100
@@ -0,0 +1,58 @@
+function plot(obj)
+
+%% Plot the time domain signal
+s = obj.s;
+fs = obj.fs;
+figure, plot((1:length(s))/fs,s);
+title('Audio signal')
+xlabel('time (s)');
+axis tight
+
+player = audioplayer(s,fs);
+set(player,'StartFcn',@plotTransportBar);
+set(player,'TimerFcn',@updateTransportBar);
+set(player,'StopFcn',@deleteTransportBar);
+
+%% Add playbaack controls
+playButtonH = uicontrol(gcf,'Style','pushbutton','String','play','Units',...
+    'Normalized','Position',[0.02 + 0.39 0 0.1 0.05]);
+stopButtonH = uicontrol(gcf,'Style','pushbutton','String','stop','Units',...
+    'Normalized','Position',[0.12 + 0.39 0 0.1 0.05]);
+
+set(playButtonH,'Callback',@play_callback);
+set(stopButtonH,'Callback',@stop_callback);
+
+    function play_callback(~,~)
+        if strcmpi(get(playButtonH,'String'),'play')
+            play(player,player.CurrentSample);
+            set(playButtonH,'String','pause');
+        else
+            pause(player)
+            set(playButtonH,'String','play');
+        end
+    end
+
+    function stop_callback(~,~)
+        stop(player);
+    end
+
+%% Transport Bar functions
+    function plotTransportBar(~,~)
+        global tbH
+        xLim = get(gca,'Xlim');
+        yLim = get(gca,'YLim');
+        tbH = line([xLim(1) xLim(1)],yLim,'Color','k');
+    end
+        
+    function updateTransportBar(hObject,~)
+        global tbH
+        currentSample = hObject.CurrentSample;
+        pos = currentSample/fs;
+        set(tbH,'XData',[pos pos]);
+    end
+
+    function deleteTransportBar(~,~)
+        global tbH
+        delete(tbH);
+    end
+end
diff -r 23763c5fbda5 -r e3035d45d014 util/classes/@dictionary/cumcoherence.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/classes/@dictionary/cumcoherence.m	Wed Aug 31 10:53:10 2011 +0100
@@ -0,0 +1,19 @@
+function mu = cumcoherence(obj)
+obj = normalize(obj);
+[M N] = size(obj.phi);
+mu = zeros(M,1);
+for m=1:M
+    c = zeros(N);
+    for i=1:N
+        c(:,i) = abs(obj.phi'*obj.phi(:,i));
+        c(i,i) = 0;
+    end
+    c = sort(c,'descend');
+    c = c(1:m,:);
+    if m==1
+        mu(m) = max(c);
+    else
+        mu(m) = max(sum(c));
+    end
+end
+end
diff -r 23763c5fbda5 -r e3035d45d014 util/classes/@dictionary/dictionary.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/classes/@dictionary/dictionary.m	Wed Aug 31 10:53:10 2011 +0100
@@ -0,0 +1,114 @@
+classdef dictionary
+    %% Dictionary for sparse representation
+    properties
+        phi         %Matrix containing the dictionary
+        len         %Length of basis functions
+        nAtoms      %Number of basis function
+        name        %String containing the matrix ensemble from which the dictionary is drawn
+    end
+    properties (Dependent = true)
+        redundancy      %Redundancy of the dictionary: nAtoms/len
+        coherence       %Maximum inner product of different basis
+        isNormalised    %True if the atoms have unit norm
+        rank            %rank of the dictionary
+    end
+    
+    methods
+        %% Constructor
+        function obj = dictionary(phi,len,nAtoms)
+            % obj = dictionary(phi,len,nAtoms)
+            % INPUTS:
+            % - phi: either a string specifying a matrix ensamble or a
+            % matrix defining an explicit dictionary
+            % - len: length of the atoms (only for implicit dictionaries)
+            % - nAtoms: number of atoms (only for implicit dictionaries)
+            if nargin
+                if ~ischar(phi)
+                    [obj.len obj.nAtoms] = size(phi);
+                    obj.phi              = phi;
+                    obj.name             = 'explicit';
+                else
+                    switch lower(phi)
+                        case 'dct'
+                            obj.phi = dctmatrix(len,nAtoms);
+                        case 'grassmanian'
+                            obj.phi = grassmanian(len,nAtoms);
+                        otherwise
+                            obj.phi = MatrixEnsemble(len,nAtoms,phi);
+                    end
+                    obj.len    = len;
+                    obj.nAtoms = nAtoms;
+                    obj.name   = lower(phi);
+                end
+            end
+        end
+        
+        %% Dependent properties
+        function redundancy = get.redundancy(obj)
+            redundancy = obj.nAtoms/obj.len;
+        end
+        function coherence = get.coherence(obj)
+            obj.phi = normcol(obj.phi);
+            G = obj.phi'*obj.phi;
+            G = G - eye(size(G));
+            coherence = max(abs(G(:)));
+        end
+        function isNormalised = get.isNormalised(obj)
+            isNormalised = norm(sum(conj(obj.phi).*obj.phi) - ...
+                ones(1,obj.nAtoms))<1e-9;
+        end
+        function r = get.rank(obj)
+            r = rank(obj.phi);
+        end
+        %% Operations
+        function obj = normalize(obj)
+            obj.phi = normcol(obj.phi);
+        end
+        
+        %% Visualization
+        function image(obj)
+            %Image of the dictionary
+            if isreal(obj.phi)
+                imagesc(obj.phi);
+                title('Dictionary');
+                xlabel('Atom number');
+            else
+                subplot(2,1,1)
+                imagesc(real(obj.phi));
+                title('Real');
+                xlabel('Atom number');
+                subplot(2,1,2)
+                imagesc(imag(obj.phi));
+                title('Imaginary');
+                xlabel('Atom number');
+            end
+        end
+        function imagegram(obj)
+            G = obj.phi'*obj.phi;
+            imagesc(G);
+            title('Gram Matrix')
+        end
+        function plot(obj,n)
+            %Plot of the n-th basis
+            if isreal(obj.phi)
+                plot(obj.phi(:,n));
+                title(['Atom number ' num2str(n) '/' num2str(size(obj.phi,2))]);
+            else
+                subplot(2,1,1)
+                plot(real(obj.phi(:,n)));
+                title(['Atom number ' num2str(n) '/' num2str(size(obj.phi,2)) ' - Real']);
+                subplot(2,1,2)
+                plot(imag(obj.phi(:,n)));
+                title(['Atom number ' num2str(n) '/' num2str(size(obj.phi,2)) ' - Imaginary']);
+            end
+        end
+        
+        function movie(obj)
+            %Movie of the basis
+            for i=1:size(obj.phi,2)
+                obj.plot(i);
+                pause(1/25);
+            end
+        end
+    end
+end
diff -r 23763c5fbda5 -r e3035d45d014 util/classes/@dictionary/mtimes.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/classes/@dictionary/mtimes.m	Wed Aug 31 10:53:10 2011 +0100
@@ -0,0 +1,11 @@
+function C = mtimes(A,B)
+isAdic = strcmp(class(A),'dictionary');
+isBdic = strcmp(class(B),'dictionary');
+if isAdic && ~isBdic
+    C = A.phi*B;
+elseif isBdic && ~isAdic
+    C = A*B.phi;
+elseif isAdic && isBdic
+    C = A.phi*B.phi;
+end
+end
\ No newline at end of file
diff -r 23763c5fbda5 -r e3035d45d014 util/classes/@dictionary/plotcumcoherence.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/classes/@dictionary/plotcumcoherence.m	Wed Aug 31 10:53:10 2011 +0100
@@ -0,0 +1,10 @@
+function plotcumcoherence(obj)
+mu = cumcoherence(obj);
+v = conv(mu,[1 1]);
+ind = find(v<1, 1, 'last');
+plot(1:obj.len,mu);
+hold on
+line([ind ind], [min(mu) max(mu)],'Color','red');
+title(['Minimum allowed sparsity (Tanner):' num2str(ind/obj.len)]);
+grid on
+end
\ No newline at end of file
diff -r 23763c5fbda5 -r e3035d45d014 util/classes/@dictionary/spark.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/classes/@dictionary/spark.m	Wed Aug 31 10:53:10 2011 +0100
@@ -0,0 +1,35 @@
+function varargout = Spark(obj)
+% Calculates the minimum number of linearly dependent atoms in the matrix A
+% WARNING: this function computes a combinatorial function, use only if the
+% size of the problem is small (i.e. <20)
+if nargout <= 1
+    A = obj.phi;
+    k = size(A,2);
+    if k>20
+        warning('This function computes a combinatorial function, use only if thesize of the problem is small (i.e. <20).');
+        choice = input('The calculation of spark will take a long time... do you wish to continue anyway (y/n)','s');
+        if strcmpi( choice,'n')
+            return
+        end
+    end
+    sigma = 2;
+    while true
+        P = nchoosek(1:k,sigma);
+        for i=1:size(P,1)
+            r = rank(A(:,P(i,:)));
+            if r<sigma
+                varargout{1} = sigma;
+                return
+            end
+        end
+        sigma = sigma + 1;
+        if sigma==k
+            varargout{1} = inf;
+            return
+        end
+    end
+else
+    %% TODO: calculate lower and upper bounds on the spark
+    varargout{1} = 2;
+    varargout{2} = inf;
+end
diff -r 23763c5fbda5 -r e3035d45d014 util/classes/dictionaryMatrices/dctmatrix.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/classes/dictionaryMatrices/dctmatrix.m	Wed Aug 31 10:53:10 2011 +0100
@@ -0,0 +1,24 @@
+function D = dctmatrix(N,K,type)
+
+error(nargchk(1,3,nargin,'struct'));
+if ~exist('type','var') || isempty(type), type='II'; end
+if ~exist('K','var') || isempty(K), K=N; end
+
+[c r] = meshgrid(0:K-1,0:N-1);
+switch type
+    case 'I'
+        D = cos(pi*c.*r/(K-1));
+        D(1,:) = D(1,:)/sqrt(2);
+        D(end,:) = D(end,:)/sqrt(2);
+    case 'II'
+        D = cos(pi*(2*c+1).*r/(2*K));
+        D(1,:) = D(1,:)/sqrt(2);
+    case 'III'
+        D = cos(pi*(2*r+1).*c/(2*K));
+        D(:,1) = D(:,1)/sqrt(2);
+    case 'IV'
+        D = cos(pi*(2*r+1+2*c+4*c.*r)/(4*K));
+    otherwise
+        error('unsupported dct type');
+end
+D = normcol(D);
\ No newline at end of file
diff -r 23763c5fbda5 -r e3035d45d014 util/classes/dictionaryMatrices/grassmannian.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/classes/dictionaryMatrices/grassmannian.m	Wed Aug 31 10:53:10 2011 +0100
@@ -0,0 +1,44 @@
+function [A G res muMin] = grassmannian(n,m,nIter,dd1,dd2,initA,verb)
+% grassmanian attempts to create an n by m matrix with minimal mutual
+% coherence using an iterative projection method.
+%
+% [A G res] = grassmanian(n,m,nIter,dd1,dd2,initA)
+%
+%
+%% Parameters and Defaults
+error(nargchk(2,7,nargin));
+
+if ~exist('verb','var')  || isempty(verb),  verb = false; end %verbose output
+if ~exist('initA','var') || isempty(initA), initA = randn(n,m); end %initial matrix
+if ~exist('dd2','var')   || isempty(dd2),   dd2 = 0.95; end %shrinking factor
+if ~exist('dd1','var')   || isempty(dd1),   dd1 = 0.9; end %percentage of coherences to be shrinked
+if ~exist('nIter','var') || isempty(nIter), nIter = 5; end %number of iterations
+
+%% Compute svd and gramian
+A = normc(initA); %normalise columns
+[Uinit Sigma] = svd(A); %calculate svd of the matrix
+G = A'*A; %gramian matrix
+
+muMin = sqrt((m-n)/(n*(m-1)));              %Lower bound on mutual coherence
+res = zeros(nIter,1);
+for iIter = 1:nIter
+   gg  = sort(abs(G(:))); %sort inner products from less to ost correlated
+   pos = find(abs(G(:))>gg(round(dd1*(m^2-m))) & abs(G(:)-1)>1e-6); 
+   G(pos) = G(pos)*dd2;
+   [U S V] = svd(G);
+   S(n+1:end,1+n:end) = 0;
+   G = U*S*V';
+   G = diag(1./abs(sqrt(diag(G))))*G*diag(1./abs(sqrt(diag(G))));
+   gg = sort(abs(G(:)));
+   pos = find(abs(G(:))>gg(round(dd1*(m^2-m))) & abs(G(:)-1)>1e-6);
+   res(iIter) = max(abs(G(pos)));
+   if verb
+    fprintf(1,'%6i  %12.8f   %12.8f  %12.8f \n',...
+            [iIter,muMin,mean(abs(G(pos))),max(abs(G(pos)))]);
+   end
+end
+
+[~, Sigma_gram V_gram] = svd(G); %calculate svd decomposition of gramian
+Sigma_new = sqrt(Sigma_gram(1:n,:)).*sign(Sigma);
+A = Uinit*Sigma_new*V_gram';
+A = normc(A);
diff -r 23763c5fbda5 -r e3035d45d014 util/classes/util/regularisedictionary.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/classes/util/regularisedictionary.m	Wed Aug 31 10:53:10 2011 +0100
@@ -0,0 +1,41 @@
+function Dreg = regularisedictionary(D)
+% REGULARISEDICTIONARY regularise dictionary for sparse representation
+%
+% Example
+% Dreg = regularisedictionary(D)
+%
+% Input
+% - D: initial dictionary matrix.
+%
+% Output
+% - Dreg: regularised dictionary
+%
+% References:
+%
+% See also
+
+% Author(s): Daniele Barchiesi
+% Copyright 2011-2011
+
+[n m] = size(D);	%dimension and number of atoms
+[U S V] = svd(D);	%SVD decomposition of the dictionary
+G = V*(S'*S)*V';	%gram matrix (equivalent to D'*D)
+I = eye(m);			%identity matrix
+lambda = 1;		    %regularisation coefficient
+
+% Optimise the gram matrix
+cvx_begin							%start cvx optimisation problem
+	variable Greg(m,m) symmetric;	%declare optimisation variable to be an mxm symmetric matrix
+	expression residual(m,m)		%declare residual as an intermediate calculation step
+	residual = G-Greg;				%define residual
+	minimise (max(max(abs(Greg-I))));	%define objective function
+	subject to						%declare constraints
+		Greg == semidefinite(m);	%positive semidefinite cone membership
+		diag(Greg) == ones(m,1);	%unit diagonal
+		norm(residual,'fro') <= lambda;	%
+cvx_end
+
+[~, Sgram Vgram] = svd(Greg);		%SVD decomposition of gramian
+Snew = sqrt(Sgram(1:n,:)).*sign(S); %calculate singular values of the regularised dictionary
+Dreg = U*Snew*Vgram';				%calculate regularised dictionary
+Dreg = normc(Dreg);					%normalise columns of regularised dictionary