V_LOGSUM v_logsum(x,d,k)=log(sum(k.*exp(x),d)) Usage: (1) y=v_logsum(x) % log(sum(exp(x))) (2) f=0.1*log(10); y=logsm(x*f)/f; % add powers expressed in dB Identities: (1) v_logsum(x+c)=v_logsum(x)+c % where c is a constant (2) sum(exp(x-logsum(x)))=1 % normalize to make exp(x) sum to 1 Inputs: X(M,N,...) data matrix to sum D optional dimension to sum along [1st non-singular dimension] K(M,N,...) optional scaling matrix. It must either be idential in size to X, or else be a vector whose length is equal to the size of dimension D of X Outputs: Y(1,N,...) = log(sum(exp(X).*K,D)) This routine evaluates the given expression for Y but takes care to avoid overflow or underflow.
0001 function y=v_logsum(x,d,k) 0002 %V_LOGSUM v_logsum(x,d,k)=log(sum(k.*exp(x),d)) 0003 % 0004 % Usage: (1) y=v_logsum(x) % log(sum(exp(x))) 0005 % (2) f=0.1*log(10); y=logsm(x*f)/f; % add powers expressed in dB 0006 % 0007 % Identities: (1) v_logsum(x+c)=v_logsum(x)+c % where c is a constant 0008 % (2) sum(exp(x-logsum(x)))=1 % normalize to make exp(x) sum to 1 0009 % 0010 % Inputs: X(M,N,...) data matrix to sum 0011 % D optional dimension to sum along [1st non-singular dimension] 0012 % K(M,N,...) optional scaling matrix. It must either be idential 0013 % in size to X, or else be a vector whose length is 0014 % equal to the size of dimension D of X 0015 % 0016 % Outputs: Y(1,N,...) = log(sum(exp(X).*K,D)) 0017 % 0018 % This routine evaluates the given expression for Y but takes care to avoid 0019 % overflow or underflow. 0020 0021 % Copyright (C) Mike Brookes 1998 0022 % Version: $Id: v_logsum.m 10865 2018-09-21 17:22:45Z dmb $ 0023 % 0024 % VOICEBOX is a MATLAB toolbox for speech processing. 0025 % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 0026 % 0027 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0028 % This program is free software; you can redistribute it and/or modify 0029 % it under the terms of the GNU General Public License as published by 0030 % the Free Software Foundation; either version 2 of the License, or 0031 % (at your option) any later version. 0032 % 0033 % This program is distributed in the hope that it will be useful, 0034 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0035 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0036 % GNU General Public License for more details. 0037 % 0038 % You can obtain a copy of the GNU General Public License from 0039 % http://www.gnu.org/copyleft/gpl.html or by writing to 0040 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0041 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0042 0043 if nargin==1 || ~numel(d) 0044 d=[find(size(x)-1) 1]; 0045 d=d(1); 0046 end 0047 n=size(x,d); 0048 if n<=1, % use efficient computation if only one term in the sum 0049 if nargin<3 0050 y=x; 0051 else 0052 y=x+log(k); 0053 end 0054 return; 0055 end 0056 s=size(x); 0057 p=[d:ndims(x) 1:d-1]; 0058 z=reshape(permute(x,p),n,prod(s)/n); 0059 q=max(z,[],1); % we subtract y from each row to avoid dynamic range problems 0060 a=(q==Inf)|(q==-Inf); % check for infinities 0061 if nargin<3 0062 y=q+log(sum(exp(z-q(ones(n,1),:)),1)); 0063 elseif numel(k)==n 0064 y=q+log(sum(exp(z-q(ones(n,1),:)).*repmat(k(:),1,prod(s)/n),1)); 0065 else 0066 y=q+log(sum(exp(z-q(ones(n,1),:)).*reshape(permute(k,p),n,prod(s)/n),1)); 0067 end 0068 y(a)=q(a); % correct any column whose max is +-Inf 0069 s(d)=1; % update the dimension of the summed component 0070 y=ipermute(reshape(y,s(p)),p); 0071