Mercurial > hg > smallbox
changeset 160:e3035d45d014 danieleb
Added support classes
| author | Daniele Barchiesi <daniele.barchiesi@eecs.qmul.ac.uk> |
|---|---|
| date | Wed, 31 Aug 2011 10:53:10 +0100 |
| parents | 23763c5fbda5 |
| children | 88578ec2f94a |
| files | util/classes/@audio/audio.m util/classes/@audio/plot.m util/classes/@dictionary/cumcoherence.m util/classes/@dictionary/dictionary.m util/classes/@dictionary/mtimes.m util/classes/@dictionary/plotcumcoherence.m util/classes/@dictionary/spark.m util/classes/dictionaryMatrices/dctmatrix.m util/classes/dictionaryMatrices/grassmannian.m util/classes/util/regularisedictionary.m |
| diffstat | 10 files changed, 410 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /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
--- /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
--- /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
--- /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
--- /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
--- /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
--- /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
--- /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
--- /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);
--- /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