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