function [ W, H, cnmfY, costFunc ] = NMFconv( V, parameter ) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Name: NMFconv % Date: March 2018 % Programmer: Christian Dittmar % % Convolutive Non-Negative Matrix Factorization with Beta-Divergence and % optional regularization parameters as described in chapter 3.7 of [1]. % The averaged activation updates are computed via the compact algorithm % given in paragraph 3.7.3. For the sake of consistency, we use the notation % from [2] instead of the one from the book. % % References: % [1] Andrzej Cichocki, Rafal Zdunek, Anh Huy Phan, and Shun-ichi Amari % "Nonnegative Matrix and Tensor Factorizations: Applications to % Exploratory Multi-Way Data Analysis and Blind Source Separation" % John Wiley and Sons, 2009. % % [2] Christian Dittmar and Meinard Müller "Reverse Engineering the Amen % Break -- Score-informed Separation and Restoration applied to Drum % Recordings" IEEE/ACM Transactions on Audio, Speech, and Language Processing, % 24(9): 1531-1543, 2016. % % Input: V matrix that shall be decomposed (typically a % magnitude spectrogram of dimension % numBins x numFrames) % parameter. % numComp number of NMFD components (denoted as R in [2]) % numIter number of NMFD iterations (denoted as L in [2]) % numTemplateFrames number of time frames for the 2D-template (denoted as T in [2]) % initW an initial estimate for the templates (denoted as W^(0) in [2]) % initH an initial estimate for the gains (denoted as H^(0) in [2]) % beta the beta parameter of the divergence: % -1 -> equals Itakura Saito divergence % 0 -> equals Kullback Leiber divergence % 1 -> equals Euclidean distance % sparsityWeight strength of the activation sparsity % uncorrWeight strength of the template uncorrelatedness % % Output: W cell-array with the learned templates % H matrix with the learned activations % cnmfV cell-array with approximated component % spectrograms % costFunc the approximation quality per iteration % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % If you use the 'NMF toolbox' please refer to: % [3] Patricio López-Serrano, Christian Dittmar, Yiğitcan Özer, and Meinard % Müller % NMF Toolbox: Music Processing Applications of Nonnegative Matrix % Factorization % In Proceedings of the International Conference on Digital Audio Effects % (DAFx), 2019. % % License: % This file is part of 'NMF toolbox'. % https://www.audiolabs-erlangen.de/resources/MIR/NMFtoolbox/ % 'NMF toolbox' is free software: you can redistribute it and/or modify it % under the terms of the GNU General Public License as published by the % the Free Software Foundation, either version 3 of the License, or (at % your option) any later version. % % 'NMF toolbox' is distributed in the hope that it will be useful, but % WITHOUT ANY WARRANTY; without even the implied warranty of % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General % Public License for more details. % % You should have received a copy of the GNU General Public License along % with 'NMF toolbox'. If not, see http://www.gnu.org/licenses/. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % check parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if nargin < 2 parameter = []; end if nargin < 1 error('Please specify input matrix V.'); end if ~isfield(parameter,'numComp') parameter.numComp = 3; end if ~isfield(parameter,'numIter') parameter.numIter = 30; end if ~isfield(parameter,'numTemplateFrames') parameter.numTemplateFrames = 8; end if ~isfield(parameter,'beta') parameter.beta = 0; end if ~isfield(parameter,'sparsityWeight') parameter.sparsityWeight = 0; end if ~isfield(parameter,'uncorrWeight') parameter.uncorrWeight = 0; end % use parameter nomenclature as in [2] [K,M] = size(V); T = parameter.numTemplateFrames; R = parameter.numComp; L = parameter.numIter; beta = parameter.beta; sparsityWeight = parameter.sparsityWeight; uncorrWeight = parameter.uncorrWeight; % use initial templates if ~isfield(parameter,'initW') initW = initTemplates(parameter,'random'); else initW = parameter.initW; end % stack the templates into a tensor for r = 1:R tensorW(:,r,:) = squeeze(initW{r}); end % use initial activations if ~isfield(parameter,'initH') initH = initActivations(parameter,'uniform'); else initH = parameter.initH; end % copy initial H = initH; % this is important to prevent initial jumps in the divergence measure V = V./(eps+sum(V(:))); for iter = 1:L % compute first approximation Lambda = convModel(tensorW,H); LambdaBeta = Lambda.^beta; Q = V.*LambdaBeta./Lambda; costMat = V.*(V.^beta - Lambda.^beta)/(eps+beta) - (V.^(beta+1) - Lambda.^(beta+1))/(eps+beta+1); costFunc(iter) = mean(costMat(:)); for t = 1:T % respect zero index tau = t-1; % use tau for shifting and t for indexing transpH = shiftOperator(H,tau)'; numeratorUpdateW = Q*transpH; denominatorUpdateW = eps + LambdaBeta*transpH + uncorrWeight*sum(squeeze(tensorW(:,:,setdiff(1:T,t))),3); tensorW(:,:,t) = tensorW(:,:,t) .* (numeratorUpdateW./denominatorUpdateW); end numeratorUpdateH = convModel(permute(tensorW, [2 1 3]),fliplr(Q)); denominatorUpdateH = convModel(permute(tensorW,[2 1 3]),fliplr(LambdaBeta)) + sparsityWeight + eps; H = H .* fliplr(numeratorUpdateH./denominatorUpdateH); % normalize templates to unit sum normVec = sum(sum(tensorW,1),3); tensorW = bsxfun(@times,tensorW,1./(eps+normVec)); % subplot(2,1,1) % imagesc(H); % axis xy % drawnow % % subplot(2,1,2) % imagesc(sqrt(reshape(permute(tensorW,[1 3 2]),K,T*R))); % axis xy % drawnow % 1; end % compute final output approximation for r = 1:R W{r} = squeeze(tensorW(:,r,:)); cnmfY{r} = convModel(tensorW(:,r,:), H(r,:)); end