Mercurial > hg > smallbox
changeset 8:33850553b702
(none)
| author | idamnjanovic |
|---|---|
| date | Mon, 22 Mar 2010 10:56:54 +0000 |
| parents | 0151f1ea080d |
| children | 28f2b5fe3483 |
| files | util/AMT_analysis.m util/Pierre_reconstruct.m util/SL_A.m util/SMALL_AMT_plot.m util/SMALL_AudioDeNoiseResult.m util/SMALL_ImgDeNoiseResult.m util/SMALL_denoise.m util/SMALL_init_DL.m util/SMALL_init_solver.m util/SMALL_learn.m util/SMALL_midiGenerate.m util/SMALL_playAudio.m util/SMALL_plot.m util/SMALL_solve.m util/SMALL_swipe.m |
| diffstat | 15 files changed, 744 insertions(+), 43 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/AMT_analysis.m Mon Mar 22 10:56:54 2010 +0000 @@ -0,0 +1,49 @@ +function AMT_res = AMT_analysis(Problem, solver) +%%% Automatic Music Transcription results analysis +% Ivan Damnjanovic 2009 +% +% If wav file that is transcribed is generated from midi file (i.e. if +% groundtruth exists) transcription is comapred to the original notes and +% AMT_res structure is generated. It contains following fields: +% - tp_notes - true positive notes (notes corectly transcribed) +% - oe_notes - octave errors (erroes due to imperfect pitch estimation) +% - fn_notes_wo_oe - false negative notes without octave errors +% (notes that were not detected) +% - fp_notes_wo_oe - false positive notes without octave erors +% - TP - number of true positives +% - FN - number of false negatives +% - FP - number of false positives + +timeOr=Problem.notesOriginal(:,5); +noteOr=Problem.notesOriginal(:,3); +timeTr=solver.reconstructed.notes(:,5); +noteTr=solver.reconstructed.notes(:,3); +n=size(timeOr,1); +m=size(timeTr,1); + +% tolerance (ts) is set to one window before and after the reference offset +% time + +ts=(Problem.windowSize)/Problem.fs; + +Hits=[]; +OE=[]; + +for i=1:n + Hit= find((noteTr(:)==noteOr(i))&(abs(timeOr(i)-timeTr(:))<ts)); + if size(Hit,1)>1 Hit=Hit(1); end + if ~isempty(Hit) Hits=[Hits; i , Hit]; + else + OctErr=find(((noteTr(:)==noteOr(i)-12)|(noteTr(:)==noteOr(i)-24))&(abs(timeOr(i)-timeTr(:))<ts), 1); + if ~isempty(OctErr) OE=[OE; i , OctErr]; end + end +end + +AMT_res.tp_notes = [Problem.notesOriginal(Hits(:,1),[3 5]) solver.reconstructed.notes(Hits(:,2),[3 5])]; +AMT_res.oe_notes = [Problem.notesOriginal(OE(:,1),[3 5]) solver.reconstructed.notes(OE(:,2),[3 5])]; +AMT_res.fn_notes_wo_oe = Problem.notesOriginal(setdiff([1:n],union(Hits(:,1),OE(:,1))),[3 5]); +AMT_res.fp_notes_wo_oe = solver.reconstructed.notes(setdiff([1:m],union(Hits(:,2),OE(:,2))),[3 5]); +AMT_res.TP=size(Hits,1); +AMT_res.FN=n-AMT_res.TP; +AMT_res.FP=m-AMT_res.TP; +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/Pierre_reconstruct.m Mon Mar 22 10:56:54 2010 +0000 @@ -0,0 +1,22 @@ +function reconstructed=Pierre_reconstruct(y, Problem) +%% Pierre Villars Example - reconstruction function +% This example is based on the experiment suggested by Professor Pierre +% Vandergheynst on the SMALL meeting in Villars. + +% using sparse representation y in dictionary Problem.A reconstruct the +% patches from the target image + +imout=Problem.A*y; + +% combine the patches into reconstructed image + +im=col2im(imout,Problem.blocksize,size(Problem.imageTrg),'disctint'); + +% bound the pixel values to [0,255] range +im(im<0)=0; +im(im>255)=255; + +%% output structure image+psnr %% +reconstructed.image=im; +reconstructed.psnr = 20*log10(Problem.maxval * sqrt(numel(Problem.imageTrg(:))) / norm(Problem.imageTrg(:)-im(:))); +end \ No newline at end of file
--- a/util/SL_A.m Mon Mar 22 10:47:17 2010 +0000 +++ b/util/SL_A.m Mon Mar 22 10:56:54 2010 +0000 @@ -1,4 +1,4 @@ - function y = SL_A(mode, m, n, x, I, dim) + function y = SL_A(A, mode, m, n, x, I, dim) % Ivan Damnjanovic 2009 % This is auxilary function to allow implicit matrices from SPARCO % to be used with SparsLab solvers @@ -8,11 +8,11 @@ u = zeros(dim, 1); u(I) = x; - y = SMALL.Problem.A(u,1); + y = A(u,1); elseif (mode == 2) - x2 = SMALL.Problem.A(x,2); + x2 = A(x,2); y = x2(I); end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/SMALL_AMT_plot.m Mon Mar 22 10:56:54 2010 +0000 @@ -0,0 +1,21 @@ +function figAMT=SMALL_AMT_plot(SMALL, AMT_res) +% Ivan Damnjanovic 2010 +% Function gets as input SMALL structure and plots AMT + + + figAMT=figure('Name','Automatic Music Transcription'); + + m=size(AMT_res,2); + +for i =1:m + subplot(m,1, i);plot(AMT_res(i).tp_notes(:,2), AMT_res(i).tp_notes(:,1),'ko', ... + AMT_res(i).tp_notes(:,4), AMT_res(i).tp_notes(:,3),'gx', ... + AMT_res(i).oe_notes(:,2), AMT_res(i).oe_notes(:,1),'bo', ... + AMT_res(i).oe_notes(:,4), AMT_res(i).oe_notes(:,3),'bx', ... + AMT_res(i).fn_notes_wo_oe(:,2), AMT_res(i).fn_notes_wo_oe(:,1),'ro', ... + AMT_res(i).fp_notes_wo_oe(:,2), AMT_res(i).fp_notes_wo_oe(:,1),'rx') + title(sprintf('%s dictionary in %.2f s - TP=%d FN=%d (Octave Errors = %d) FP=%d', SMALL.DL(i).name, SMALL.DL(i).time, AMT_res(i).TP, AMT_res(i).FN, size(AMT_res(i).oe_notes,1), AMT_res(i).FP)); + xlabel('Time') + ylabel('Note Number') + +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/SMALL_AudioDeNoiseResult.m Mon Mar 22 10:56:54 2010 +0000 @@ -0,0 +1,36 @@ +function SMALL_AudioDeNoiseResult(SMALL) + +fMain=figure('Name', sprintf('File %s (training set size- %d, sigma - %d)',SMALL.Problem.name, SMALL.Problem.n, SMALL.Problem.sigma)); +m=size(SMALL.solver,2); +maxval=SMALL.Problem.maxval; +au=SMALL.Problem.Original; +aunoise=SMALL.Problem.Noisy; + +subplot(2, m, 1); plot(au/maxval); +title('Original audio'); + +subplot(2,m,2); plot(aunoise/maxval); +title(sprintf('Noisy audio, PSNR = %.2fdB', 20*log10(maxval * sqrt(numel(au)) / norm(au(:)-aunoise(:))) )); + +for i=1:m + params=SMALL.solver(i).param; + sWav=subplot(2, m, m+i, 'Parent', fMain); plot(SMALL.solver(i).reconstructed.Image/maxval, 'Parent', sWav); + title(sprintf('%s Denoised audio, PSNR: %.2fdB', SMALL.DL(i).name, SMALL.solver(i).reconstructed.psnr),'Parent', sWav ); + if strcmpi(SMALL.DL(i).name,'ksvds') + D = kron(SMALL.Problem.basedict{2},SMALL.Problem.basedict{1})*SMALL.DL(i).D; + else + D = SMALL.DL(i).D; + end + figure('Name', sprintf('%s dictionary in %.2f s', SMALL.DL(i).name, SMALL.DL(i).time)); + imshow(D*255); +% n= size(D,2); +% sqrtn=round(sqrt(size(D,2))); +% for j=1:n +% subplot(sqrtn,sqrtn,j); plot(D(:,j)); +% end +% dictimg = showdict(D,[params.blocksize 1],round(sqrt(size(D,2))),round(sqrt(size(D,2))),'lines','highcontrast'); +% +% subplot(2,m,m+i);imshow(imresize(dictimg,2,'nearest')); +% title(sprintf('%s dictionary in %.2f s', SMALL.DL(i-1).name, SMALL.DL(i-1).time)); + +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/SMALL_ImgDeNoiseResult.m Mon Mar 22 10:56:54 2010 +0000 @@ -0,0 +1,38 @@ +function SMALL_ImgDeNoiseResult(SMALL) +% Ivan Damnjanovic 2010 +% Function gets as input SMALL structure and plots Image Denoise +% results: Original Image, Noisy Image and for learned dictionaries and +% denoised images + + +figure('Name', sprintf('Image %s (training set size- %d, sigma - %d)',SMALL.Problem.name, SMALL.Problem.n, SMALL.Problem.sigma)); + +m=size(SMALL.solver,2)+1; +maxval=SMALL.Problem.maxval; +im=SMALL.Problem.Original; +imnoise=SMALL.Problem.Noisy; + +subplot(2, m, 1); imshow(im/maxval); +title('Original image'); + +subplot(2,m,m+1); imshow(imnoise/maxval); +title(sprintf('Noisy image, PSNR = %.2fdB', 20*log10(maxval * sqrt(numel(im)) / norm(im(:)-imnoise(:))) )); + +for i=2:m + params=SMALL.solver(i-1).param; + subplot(2, m, i); imshow(SMALL.solver(i-1).reconstructed.Image/maxval); + title(sprintf('%s Denoised image, PSNR: %.2f dB in %.2f s',... + SMALL.DL(i-1).name, SMALL.solver(i-1).reconstructed.psnr, SMALL.solver(i-1).time )); + if strcmpi(SMALL.DL(i-1).name,'ksvds') + D = kron(SMALL.Problem.basedict{2},SMALL.Problem.basedict{1})*SMALL.DL(i-1).D; + else + D = SMALL.DL(i-1).D; + end + dictimg = showdict(D,params.blocksize,... + round(sqrt(size(D,2))),round(sqrt(size(D,2))),'lines','highcontrast'); + + subplot(2,m,m+i);imshow(imresize(dictimg,2,'nearest')); + title(sprintf('%s dictionary in %.2f s',... + SMALL.DL(i-1).name, SMALL.DL(i-1).time)); + +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/SMALL_denoise.m Mon Mar 22 10:56:54 2010 +0000 @@ -0,0 +1,77 @@ +function solver=SMALL_denoise(Problem, solver) + +%%% SMALL denoising +% Ivan Damnjanovic 2009 +% Function gets as input SMALL structure that contains SPARCO problem to +% be solved, name of the toolbox and solver, and parameters file for +% particular solver. +% It is based on omp/omps denoising by Ron Rubenstein (KSVD toolbox) +% +% Outputs are denoised image, psnr, number of non-zero coeficients and +% time spent +%% + + %%%%% denoise the signal %%%%% + fprintf('\nStarting %s... \n', solver.name); + start=cputime; + %% + if strcmpi(solver.toolbox,'ompbox') + if (~isfield(solver.param,'lambda')) + solver.param.lambda = Problem.maxval/(10*Problem.sigma); + end + + solver.param = Problem; + %solver.param.memusage = 'high'; + solver.param = rmfield(solver.param, {'Noisy' 'Original' 'b' 'm' 'n' 'p' 'initdict'}); + solver.param.x = Problem.Noisy; + solver.param.dict = Problem.A; + p = Problem.signalDim; + msgdelta=5; + % call the appropriate ompdenoise function + if (p==1) + [y,nz] = ompdenoise1(solver.param,msgdelta); + elseif (p==2) + [y,nz] = ompdenoise2(solver.param,msgdelta); + elseif (p==3) + [y,nz] = ompdenoise3(solver.param,msgdelta); + else + [y,nz] = ompdenoise(solver.param,msgdelta); + end + elseif strcmpi(solver.toolbox,'ompsbox') + if (~isfield(solver.param,'lambda')) + solver.param.lambda = Problem.maxval/(10*Problem.sigma); + end + + solver.param = Problem; + %solver.param.memusage = 'high'; + solver.param = rmfield(solver.param, {'Noisy' 'Original' 'b' 'm' 'n' 'p' 'initdict'}); + solver.param.x = Problem.Noisy; + if issparse(Problem.A) + solver.param.A = Problem.A; + else + solver.param.A = sparse(Problem.A); + end + p = Problem.signalDim; + msgdelta=5; + % call the appropriate ompdenoise function + if (p==1) + [y,nz] = ompsdenoise1(solver.param,msgdelta); + elseif (p==2) + [y,nz] = ompsdenoise2(solver.param,msgdelta); + elseif (p==3) + [y,nz] = ompsdenoise3(solver.param,msgdelta); + else + [y,nz] = ompsdenoise(solver.param,msgdelta); + end + else + printf('\nToolbox has not been registered. Please change SMALL_learn file.\n'); + return + end + %% + solver.time = cputime - start; + fprintf('\n%s finished task in %2f seconds. \n', solver.name, solver.time); + solver.reconstructed.Image=y; + solver.reconstructed.psnr=20*log10(Problem.maxval * sqrt(numel(Problem.Original)) / norm(Problem.Original(:)-solver.reconstructed.Image(:))); + solver.reconstructed.nz=nz; + +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/SMALL_init_DL.m Mon Mar 22 10:56:54 2010 +0000 @@ -0,0 +1,15 @@ +function DL = SMALL_init_DL(varargin) +% Ivan Damnjanovic 2010 + % Function initialise SMALL structure for Dictionary Learning. + % Optional input variables: + % toolbox - name of Dictionary Learning toolbox you want to use + % name - name of the algorithm from above toolbox + % param - parameters you want to set + %% + + DL.toolbox=[]; + DL.name=[]; + DL.param=[]; + DL.D=[]; + DL.time=[]; + end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/SMALL_init_solver.m Mon Mar 22 10:56:54 2010 +0000 @@ -0,0 +1,16 @@ +function solver = SMALL_init_solver(varargin) +% Ivan Damnjanovic 2010 + % Function initialise SMALL structure for Dictionary Learning. + % Optional input variables: + % toolbox - name of Dictionary Learning toolbox you want to use + % name - name of the algorithm from above toolbox + % param - parameters you want to set + %% + + solver.toolbox=[]; + solver.name=[]; + solver.param=[]; + solver.solution=[]; + solver.reconstructed=[]; + solver.time=[]; +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/SMALL_learn.m Mon Mar 22 10:56:54 2010 +0000 @@ -0,0 +1,72 @@ +function DL = SMALL_learn(Problem,DL) +%%% SMALL Dictionary Learning +% Ivan Damnjanovic 2009 +% Function gets as input Problem and Dictionary Learning (DL) structures +% In Problem structure field b with the training set needs to be defined +% In DL fields with name of the toolbox and solver, and parameters file +% for particular dictionary learning technique needs to be present. +% +% Outputs are Learned dictionary and time spent as a part of DL structure +%% + + fprintf('\nStarting Dictionary Learning %s... \n', DL.name); + start=cputime; + + if strcmpi(DL.toolbox,'KSVD') + param=DL.param; + param.data=Problem.b; + + D = eval([DL.name,'(param, ''t'', 5);']); + elseif strcmpi(DL.toolbox,'KSVDS') + param=DL.param; + param.data=Problem.b; + + D = eval([DL.name,'(param, ''t'', 5);']); + elseif strcmpi(DL.toolbox,'SPAMS') + + X = Problem.b; + param=DL.param; + + D = eval([DL.name,'(X, param);']); + % As some versions of SPAMS does not produce unit norm column + % dictionaries, we need to make sure that columns are normalised to + % unit lenght. + + for i = 1: size(D,2) + D(:,i)=D(:,i)/norm(D(:,i)); + end + +% To introduce new dictionary learning technique put the files in +% your Matlab path. Next, unique name <TolboxID> for your toolbox needs +% to be defined and also prefferd API for toolbox functions <Preffered_API> +% +% elseif strcmpi(DL.toolbox,'<ToolboxID>') +% % This is an example of API that can be used: +% % - get training set from Problem part of structure +% % - assign parameters defined in the main program +% +% X = Problem.b; +% param=DL.param; +% +% % - Evaluate the function (DL.name - defined in the main) with +% % parameters given above +% +% D = eval([DL.name,'(<Preffered_API>);']); + + else + printf('\nToolbox has not been registered. Please change SMALL_learn file.\n'); + return + end + +%% +% Dictionary Learning time + + DL.time = cputime - start; + fprintf('\n%s finished task in %2f seconds. \n', DL.name, DL.time); + +% If dictionary is given as a sparse matrix change it to full + + DL.D = full(D); + +end + \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/SMALL_midiGenerate.m Mon Mar 22 10:56:54 2010 +0000 @@ -0,0 +1,109 @@ +function reconstructed=SMALL_midiGenerate(V, Problem) +%%% Reconstraction of midi file from representation in the given dictionary +% Ivan Damnjanovic 2009 +% +% SMALL_midiGenerate is a part of SMALLbox and can be use to reconstruct +% a midi file given representation of the training set (V) in the +% dictionary Problem.A. +% Output is reconstructed structure with two fields: +% - reconstructed.notes - matrix with transcribed notes +% - reconstructed.midi - midi representation of transcription + + +U=Problem.A; % Dictionary used for representation +fs=Problem.fs; % Sampling rate +f=Problem.f; % vector of frequencies at wihch spectrogram is computed + +ts=(Problem.windowSize*Problem.overlap)/fs; %size of an analysis frame in seconds + +%% +% Components pitch estimation using modified SWIPE algorithm by Arthuro +% Camacho +% +% Columns of matrix U are spectrograms of the notes learned from the +% training set. We are estimating pitches of these notes by also +% restricting pitch values to the one of the 88 piano notes. + +pitch=zeros(size(U,2),1); + +for i=1:size(U,2) + + pitch(i) = SMALL_swipe(U(:,i),fs, f, [27.50 8192], 1/12); + +end + +%% +% If some of columns of U have the same pitch, their contribution to the +% score (matrix V) is summed. + +[Ps,idx]=sort(pitch); +ndp=1; +Pd(ndp)=Ps(1); +Vnew(ndp,:)=V(idx(1),:); +for i=2:88 + if Ps(i)> Ps(i-1) + + ndp=ndp+1; + Vnew(ndp,:)=V(idx(i),:); + Pd(ndp)=Ps(i); + + else + Vnew(ndp,:)=Vnew(ndp,:)+V(idx(i),:); + end +end +%% +% Generate midi matrix + +midx=0; +for i=1:ndp + + % Threshold for finding onsets and offsets of notes + + thr=mean(Vnew(i,:));%+std(Vnew(i,:)); + + if(Pd(i)~=0) + for j=1:size(Vnew,2) + if Vnew(i,j)<thr + Vnew(i,j)=0; + if(j>1) + if (Vnew(i,j-1)==1) + try + M(midx,6)=(j-1)*ts; + if (M(midx,6)-M(midx,5))<2*ts + midx=midx-1; + end + catch + pause; + end + end + end + else + Vnew(i,j)=1; + if(j>1) + if (Vnew(i,j-1)==0) + midx=midx+1; + M(midx,1)=1; + M(midx,2)=1; + M(midx,3)=69 +round( 12 *log2(Pd(i)/440)); + M(midx,4)=80; + M(midx,5)=(j-1)*ts; + end + else + midx=midx+1; + M(midx,1)=1; + M(midx,2)=1; + M(midx,3)=69 + round(12 *log2(Pd(i)/440)); + M(midx,4)=80; + M(midx,5)=0; + end + end + end + if M(midx,6)==0 + M(midx,6)=(j-1)*ts; + end + end +end + +M=sortrows(M,5); +reconstructed.notes=M; +reconstructed.midi = matrix2midi(M); \ No newline at end of file
--- a/util/SMALL_playAudio.m Mon Mar 22 10:47:17 2010 +0000 +++ b/util/SMALL_playAudio.m Mon Mar 22 10:56:54 2010 +0000 @@ -4,6 +4,7 @@ % the reconstructed signal +SMALL.solver.reconstructed = SMALL.Problem.reconstruct(SMALL.solver.solution); ch=''; while 1 request = input('\nWhat do you want to hear? \n 1. Original signal \n 2. Mixed \n 3. Reconstructed signal \n 4. Quit player\n','s');
--- a/util/SMALL_plot.m Mon Mar 22 10:47:17 2010 +0000 +++ b/util/SMALL_plot.m Mon Mar 22 10:56:54 2010 +0000 @@ -4,20 +4,22 @@ % reconstructed signal figure; - - plot(1:length(SMALL.solver.solution), SMALL.solver.solution, 'b'); - title(['Coefficients of the reconstructed signal using ', SMALL.solver.name, ' with parameters ', SMALL.solver.param]) + + m=size(SMALL.solver,2); + n=size(SMALL.solver(1).reconstructed,2)+1; + for i =1:m + + subplot(m,n, (i-1)*n+1); plot(1:length(SMALL.solver(i).solution), SMALL.solver(i).solution, 'b') + title([ SMALL.solver(i).name,'(', SMALL.solver(i).param,')']) xlabel('Coefficient') - % Use SMALL.solution to reconstruct the signal. - % Plot reconstructed signal against original - for i=1:size(SMALL.solver.reconstructed,2) - figure; - plot(1:length(SMALL.solver.reconstructed(:,i)), SMALL.solver.reconstructed(:,i) ,'b.-', 1:length(SMALL.Problem.signal(:,i)), SMALL.Problem.signal(:,i),'r--'); - legend(SMALL.solver.name,'Original signal'); + for j=2:n + + subplot(m,n,(i-1)*n+j); plot(1:length(SMALL.solver(i).reconstructed(:,j-1)), SMALL.solver(i).reconstructed(:,j-1) ,'b.-', 1:length(SMALL.Problem.signal(:,j-1)), SMALL.Problem.signal(:,j-1),'r--') + %legend(SMALL.solver(i).name,'Original signal',0); title('Reconstructed and original signals'); end end
--- a/util/SMALL_solve.m Mon Mar 22 10:47:17 2010 +0000 +++ b/util/SMALL_solve.m Mon Mar 22 10:56:54 2010 +0000 @@ -1,32 +1,77 @@ -function SMALL = SMALL_solve(SMALL) - % Ivan Damnjanovic 2009 - % Function gets as input SMALL structure that contains SPARCO problem to - % be solved, name of the toolbox and solver, and parameters file for particular solver - % Outputs are solution, reconstructed signal and time spent - %% - - A = @(x) SMALL.Problem.A(x,1); % The operator - AT = @(y) SMALL.Problem.A(y,2); % and its transpose. - b = SMALL.Problem.b; % The right-hand-side vector. - m = SMALL.Problem.sizeA(1); % m is the no. of rows. - n = SMALL.Problem.sizeA(2); % n is the no. of columns. - - fprintf('\nStarting solver %s... \n', SMALL.solver.name); - start=cputime; - %% - if strcmpi(SMALL.solver.toolbox,'sparselab') - y = eval([SMALL.solver.name,'(''SL_A'', b, n,',SMALL.solver.param,');']); - elseif strcmpi(SMALL.solver.toolbox,'sparsify') - y = eval([SMALL.solver.name,'(b, A, n, ''P_trans'', AT,',SMALL.solver.param,');']); - elseif (strcmpi(SMALL.solver.toolbox,'spgl1')||strcmpi(SMALL.solver.toolbox,'gpsr')) - y = eval([SMALL.solver.name,'(b, A,',SMALL.solver.param,');']); - else - y = eval([SMALL.solver.name,'(A, b, n,',SMALL.solver.param,',AT);']); - end - %% - SMALL.solver.time = cputime - start; - fprintf('Solver %s finished task in %2f seconds. \n', SMALL.solver.name, SMALL.solver.time); - SMALL.solver.solution = full(y); - SMALL.solver.reconstructed = SMALL.Problem.reconstruct(SMALL.solver.solution); +function solver = SMALL_solve(Problem, solver) +%%% SMALL sparse solver +% Ivan Damnjanovic 2009 +% Function gets as input SMALL structure that contains SPARCO problem to +% be solved, name of the toolbox and solver, and parameters file for +% particular solver. +% +% Outputs are solution, reconstructed signal and time spent +%% + +if isa(Problem.A,'float') + A = Problem.A; + SparseLab_A=Problem.A; + m = size(Problem.A,1); % m is the no. of rows. + n = size(Problem.A,2); % n is the no. of columns. +else + A = @(x) Problem.A(x,1); % The operator + AT = @(y) Problem.A(y,2); % and its transpose. + SparseLab_A =@(mode, m, n, x, I, dim) SL_A(Problem.A, mode, m, n, x, I, dim); + m = Problem.sizeA(1); % m is the no. of rows. + n = Problem.sizeA(2); % n is the no. of columns. + end - \ No newline at end of file +b = Problem.b; % The right-hand-side vector. +%% +fprintf('\nStarting solver %s... \n', solver.name); +start=cputime; + +if strcmpi(solver.toolbox,'sparselab') + y = eval([solver.name,'(SparseLab_A, b, n,',solver.param,');']); +elseif strcmpi(solver.toolbox,'sparsify') + y = eval([solver.name,'(b, A, n, ''P_trans'', AT,',solver.param,');']); +elseif (strcmpi(solver.toolbox,'spgl1')||strcmpi(solver.toolbox,'gpsr')) + y = eval([solver.name,'(b, A,',solver.param,');']); +elseif (strcmpi(solver.toolbox,'SPAMS')) + y = eval([solver.name,'(b, A, solver.param);']); +elseif (strcmpi(solver.toolbox,'SMALL')) + if isa(Problem.A,'float') + y = eval([solver.name,'(A, b, n,',solver.param,');']); + else + y = eval([solver.name,'(A, b, n,',solver.param,',AT);']); + end + +% To introduce new sparse representation algorithm put the files in +% your Matlab path. Next, unique name <TolboxID> for your toolbox and +% prefferd API <Preffered_API> needs to be defined. +% +% elseif strcmpi(solver.toolbox,'<ToolboxID>') +% +% % - Evaluate the function (solver.name - defined in the main) with +% % parameters given above +% +% y = eval([solver.name,'(<Preffered_API>);']); + +else + printf('\nToolbox has not been registered. Please change SMALL_solver file.\n'); + return +end + +%% +% Sparse representation time + +solver.time = cputime - start; +fprintf('Solver %s finished task in %2f seconds. \n', solver.name, solver.time); + +% geting around out of memory problem when converting big matrix from +% sparse to full... + +if (Problem.sparse==1) + solver.solution = y; +else + solver.solution = full(y); +end + +% Reconstruct the signal from the solution +solver.reconstructed = Problem.reconstruct(solver.solution); +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/SMALL_swipe.m Mon Mar 22 10:56:54 2010 +0000 @@ -0,0 +1,198 @@ +function [p,s] = SMALL_swipe(X,fs, f, plim,dlog2p,dERBs,woverlap,sTHR) +% +% Ivan Damnjanovic 2010 +% +% This is modified swipep MATLAB code that is working directly in spectral +% domain and uses only one window size. The results are suboptimal +% comparing to original code. It is also converted to SWIPE which uses all +% the harmonics of the signal. +% +%SWIPEP Pitch estimation using SWIPE'. +% P = SWIPEP(X,Fs,[PMIN PMAX],DT,DLOG2P,DERBS,STHR) estimates the pitch +% of the vector signal X every DT seconds. The sampling frequency of +% the signal is Fs (in Hertz). The spectrum is computed using a Hann +% window with an overlap WOVERLAP between 0 and 1. The spectrum is +% sampled uniformly in the ERB scale with a step size of DERBS ERBs. The +% pitch is searched within the range [PMIN PMAX] (in Hertz) with samples +% distributed every DLOG2P units on a base-2 logarithmic scale of Hertz. +% The pitch is fine-tuned using parabolic interpolation with a resolution +% of 1 cent. Pitch estimates with a strength lower than STHR are treated +% as undefined. +% +% [P,T,S] = SWIPEP(X,Fs,[PMIN PMAX],DT,DLOG2P,DERBS,WOVERLAP,STHR) +% returns the times T at which the pitch was estimated and the pitch +% strength S of every pitch estimate. +% +% P = SWIPEP(X,Fs) estimates the pitch using the default settings PMIN = +% 30 Hz, PMAX = 5000 Hz, DT = 0.001 s, DLOG2P = 1/48 (48 steps per +% octave), DERBS = 0.1 ERBs, WOVERLAP = 0.5, and STHR = -Inf. +% +% P = SWIPEP(X,Fs,...,[],...) uses the default setting for the parameter +% replaced with the placeholder []. +% +% REMARKS: (1) For better results, make DLOG2P and DERBS as small as +% possible and WOVERLAP as large as possible. However, take into account +% that the computational complexity of the algorithm is inversely +% proportional to DLOG2P, DERBS and 1-WOVERLAP, and that the default +% values have been found empirically to produce good results. Consider +% also that the computational complexity is directly proportional to the +% number of octaves in the pitch search range, and therefore , it is +% recommendable to restrict the search range to the expected range of +% pitch, if any. (2) This code implements SWIPE', which uses only the +% first and prime harmonics of the signal. To convert it into SWIPE, +% which uses all the harmonics of the signal, replace the word +% PRIMES with a colon (it is located almost at the end of the code). +% However, this may not be recommendable since SWIPE' is reported to +% produce on average better results than SWIPE (Camacho and Harris, +% 2008). +% +% EXAMPLE: Estimate the pitch of the signal X every 10 ms within the +% range 75-500 Hz using the default resolution (i.e., 48 steps per +% octave), sampling the spectrum every 1/20th of ERB, using a window +% overlap factor of 50%, and discarding samples with pitch strength +% lower than 0.2. Plot the pitch trace. +% [x,Fs] = wavread(filename); +% [p,t,s] = swipep(x,Fs,[75 500],0.01,[],1/20,0.5,0.2); +% plot(1000*t,p) +% xlabel('Time (ms)') +% ylabel('Pitch (Hz)') +% +% REFERENCES: Camacho, A., Harris, J.G, (2008) "A sawtooth waveform +% inspired pitch estimator for speech and music," J. Acoust. Soc. Am. +% 124, 1638-1652. +if ~ exist( 'plim', 'var' ) || isempty(plim), plim = [30 5000]; end +%if ~ exist( 'dt', 'var' ) || isempty(dt), dt = 0.001; end +if ~ exist( 'dlog2p', 'var' ) || isempty(dlog2p), dlog2p = 1/48; end +if ~ exist( 'dERBs', 'var' ) || isempty(dERBs), dERBs = 0.05; end +% if ~ exist( 'woverlap', 'var' ) || isempty(woverlap) +% woverlap = 0.5; +% elseif woverlap>1 || woverlap<0 +% error('Window overlap must be between 0 and 1.') +% end +if ~ exist( 'sTHR', 'var' ) || isempty(sTHR), sTHR = -Inf; end +%t = [ 0: dt: length(x)/fs ]'; % Times +% Define pitch candidates +log2pc = [ log2(plim(1)): dlog2p: log2(plim(2)) ]'; +pc = 2 .^ log2pc; +S = zeros( length(pc), 1 ); % Pitch strength matrix +% Determine P2-WSs +%logWs = round( log2( 8*fs ./ plim ) ); +ws = [2822];%2.^[ logWs(1): -1: logWs(2) ]; % P2-WSs +pO = 8 * fs ./ ws; % Optimal pitches for P2-WSs +% Determine window sizes used by each pitch candidate +d = 1 + log2pc - log2( 8*fs./ws(1) ); +% Create ERB-scale uniformly-spaced frequencies (in Hertz) +fERBs = erbs2hz([ hz2erbs(min(pc)/4): dERBs: hz2erbs(fs/2) ]'); +for i = 1 : length(ws) + %dn = max( 1, round( 8*(1-woverlap) * fs / pO(i) ) ); % Hop size + % Zero pad signal + %xzp = [ zeros( ws(i)/2, 1 ); x(:); zeros( dn + ws(i)/2, 1 ) ]; + % Compute spectrum + %w = hanning( ws(i) ); % Hann window + %o = max( 0, round( ws(i) - dn ) ); % Window overlap + %[ X, f, ti ] = specgram( xzp, ws(i), fs, w, o ); + % Select candidates that use this window size + if length(ws) == 1 + j=[1:size(pc)]'; k = []; + elseif i == length(ws) + j=find(d-i>-1); k=find(d(j)-i<0); + elseif i==1 + j=find(d-i<1); k=find(d(j)-i>0); + else + j=find(abs(d-i)<1); k=1:length(j); + end + % Compute loudness at ERBs uniformly-spaced frequencies + fERBs = fERBs( find( fERBs > pc(1)/4, 1, 'first' ) : end ); + L = sqrt( max( 0, interp1( f, X, fERBs, 'spline', 0) ) ); + % Compute pitch strength + Si = pitchStrengthAllCandidates( fERBs, L, pc ); + % Interpolate pitch strength at desired times +% if size(Si,2) > 1 +% warning off MATLAB:interp1:NaNinY +% Si = interp1( ti, Si', t, 'linear', NaN )'; +% warning on MATLAB:interp1:NaNinY +% else +% Si = repmat( NaN, length(Si),1 ); +% end + % Add pitch strength to combination +% lambda = d( j(k) ) - i; + mu = ones( size(j) ); +% mu(k) = 1 - abs( lambda ); + S(j,:) = S(j,:) + repmat(mu,1,size(Si,2)) .* Si; +end +% Fine tune pitch using parabolic interpolation +p = repmat( NaN, size(S,2), 1 ); +s = repmat( NaN, size(S,2), 1 ); +for j = 1 : size(S,2) + [ s(j), i ] = max( S(:,j), [], 1 ); + if s(j) < sTHR, continue, end + if i == 1 || i == length(pc) + p(j) = pc(i); + else + I = i-1 : i+1; + tc = 1 ./ pc(I); + ntc = ( tc/tc(2) - 1 ) * 2*pi; + c = polyfit( ntc, S(I,j), 2 ); + ftc = 1 ./ 2.^[ log2(pc(I(1))): 1/12/100: log2(pc(I(3))) ]; + nftc = ( ftc/tc(2) - 1 ) * 2*pi; + [s(j) k] = max( polyval( c, nftc ) ); + p(j) = 2 ^ ( log2(pc(I(1))) + (k-1)/12/100 ); +% if (p(j)-pc(I(1)))<0.75*abs(p(j)-pc(I(2))) +% p(j)=pc(I(1)); +% elseif (pc(I(3))-p(j))<0.75*abs(p(j)-pc(I(2))) +% p(j)=pc(I(3)); +% else + p(j)=pc(I(2)); +% end + end +end + +function S = pitchStrengthAllCandidates( f, L, pc ) +% Create pitch strength matrix +S = zeros( length(pc), size(L,2) ); +% Define integration regions +k = ones( 1, length(pc)+1 ); +for j = 1 : length(k)-1 + k(j+1) = k(j) - 1 + find( f(k(j):end) > pc(j)/4, 1, 'first' ); +end +k = k(2:end); +% Create loudness normalization matrix +N = sqrt( flipud( cumsum( flipud(L.*L) ) ) ); +for j = 1 : length(pc) + % Normalize loudness + warning off MATLAB:divideByZero + NL = L(k(j):end,:) ./ repmat( N(k(j),:), size(L,1)-k(j)+1, 1); + warning on MATLAB:divideByZero + % Compute pitch strength + + S(j,:) = pitchStrengthOneCandidate( f(k(j):end), NL, pc(j) ); +end + +function S = pitchStrengthOneCandidate( f, NL, pc ) +n = fix( f(end)/pc - 0.75 ); % Number of harmonics +if n==0, S=NaN; return, end +k = zeros( size(f) ); % Kernel +% Normalize frequency w.r.t. candidate +q = f / pc; +% Create kernel +for i = [ 1:n] % primes(n) ] + a = abs( q - i ); + % Peak's weigth + p = a < .25; + k(p) = cos( 2*pi * q(p) ); + % Valleys' weights + v = .25 < a & a < .75; + k(v) = k(v) + cos( 2*pi * q(v) ) / 2; +end +% Apply envelope +k = k .* sqrt( 1./f ); +% K+-normalize kernel +k = k / norm( k(k>0) ); +% Compute pitch strength +S = k' * NL; + +function erbs = hz2erbs(hz) +erbs = 21.4 * log10( 1 + hz/229 ); + +function hz = erbs2hz(erbs) +hz = ( 10 .^ (erbs./21.4) - 1 ) * 229;