V_PPMVU calculate PPM, VU or EBU level of an audio signal [V,FX,FX1]=(X,FSX,M) Usage: (1) v=v_ppmvu(x,fs,'a') % calculate PPM level of signal x with sampling freq fs (2) [v,f]=v_ppmvu(x,fs,'aw') % calculate PPM + fast version as well (3) v=v_ppmvu(x,fs,'au') % calculate VU level in linear units rather than dB (4) v=v_ppmvu(x,fs,'c') % calculate VU level (5) v=v_ppmvu(x,fs,'e') % calculate EBU loudness level (6) [v1,fx]=v_ppmvu(x1,fs,'a') v2=v_ppmvu(x2,fx) % process in chunks: same as v_ppmvu([x1; x2],fs,'a') Inputs: x = input signal, one **column** per channel fsx = sample frequency or fx output from a previous call m = either a structure with algorithm parameters (see below) or an attack/decay time constsant or a character string: 'a' UK PPM characteristic [default] 'b' DIN PPM 'c' VU American 'd' VU French 'e' EBU short [default toggle: 'q'] 'f' EBU medium [default toggle: 'q'] followed by any combination of modifier toggles: 'z' remove mean (not yet implemented) 'p' preemphasis (not yet implemented) 's' set rise time to zero 'o' oversample x 4 'w' give fast output (with 0 decay time) as well as slow output 'v' plot graph 'q' average squared signal 'u' output magnitude (or mean square value if 'q') instead of dB Outputs: y = selected output (same size as x) fx = cell array holding algorithm state or, if 'w' option specified, the fast version of the output fx1 = cell array holding algorithm state (only if 'w' specified) Algorithm Parameters: mm text string with options ta attach time constant or zero if no attack smoothing (seconds) tm averaging filter duration or zero (seconds) td decay time constant or zero if no decay smoothing (seconds)
0001 function [v,fx,fx1]=v_ppmvu(x,fsx,m) 0002 %V_PPMVU calculate PPM, VU or EBU level of an audio signal [V,FX,FX1]=(X,FSX,M) 0003 % 0004 % Usage: (1) v=v_ppmvu(x,fs,'a') % calculate PPM level of signal x with sampling freq fs 0005 % (2) [v,f]=v_ppmvu(x,fs,'aw') % calculate PPM + fast version as well 0006 % (3) v=v_ppmvu(x,fs,'au') % calculate VU level in linear units rather than dB 0007 % (4) v=v_ppmvu(x,fs,'c') % calculate VU level 0008 % (5) v=v_ppmvu(x,fs,'e') % calculate EBU loudness level 0009 % (6) [v1,fx]=v_ppmvu(x1,fs,'a') 0010 % v2=v_ppmvu(x2,fx) % process in chunks: same as v_ppmvu([x1; x2],fs,'a') 0011 % 0012 % Inputs: x = input signal, one **column** per channel 0013 % fsx = sample frequency or fx output from a previous call 0014 % m = either a structure with algorithm parameters (see below) 0015 % or an attack/decay time constsant or a character string: 0016 % 'a' UK PPM characteristic [default] 0017 % 'b' DIN PPM 0018 % 'c' VU American 0019 % 'd' VU French 0020 % 'e' EBU short [default toggle: 'q'] 0021 % 'f' EBU medium [default toggle: 'q'] 0022 % followed by any combination of modifier toggles: 0023 % 'z' remove mean (not yet implemented) 0024 % 'p' preemphasis (not yet implemented) 0025 % 's' set rise time to zero 0026 % 'o' oversample x 4 0027 % 'w' give fast output (with 0 decay time) as well as slow output 0028 % 'v' plot graph 0029 % 'q' average squared signal 0030 % 'u' output magnitude (or mean square value if 'q') instead of dB 0031 % 0032 % Outputs: y = selected output (same size as x) 0033 % fx = cell array holding algorithm state 0034 % or, if 'w' option specified, the fast version of the output 0035 % fx1 = cell array holding algorithm state (only if 'w' specified) 0036 % 0037 % Algorithm Parameters: 0038 % mm text string with options 0039 % ta attach time constant or zero if no attack smoothing (seconds) 0040 % tm averaging filter duration or zero (seconds) 0041 % td decay time constant or zero if no decay smoothing (seconds) 0042 0043 % Copyright (C) Mike Brookes 2013 0044 % Version: $Id: v_ppmvu.m 3387 2013-08-23 12:32:47Z dmb $ 0045 % 0046 % VOICEBOX is a MATLAB toolbox for speech processing. 0047 % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 0048 % 0049 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0050 % This program is free software; you can redistribute it and/or modify 0051 % it under the terms of the GNU General Public License as published by 0052 % the Free Software Foundation; either version 2 of the License, or 0053 % (at your option) any later version. 0054 % 0055 % This program is distributed in the hope that it will be useful, 0056 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0057 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0058 % GNU General Public License for more details. 0059 % 0060 % You can obtain a copy of the GNU General Public License from 0061 % http://www.gnu.org/copyleft/gpl.html or by writing to 0062 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0063 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0064 0065 % define persistent constants 0066 persistent m4 k4 h1 g1 tamd deft p0 0067 if isempty(m4) 0068 % upsample by 4 using two cascaded half-band filters 0069 m4=5; % one quarter of the first filter order 0070 h1f=kaiser(4*m4+1,5.3)'.*sinc(-m4:0.5:m4); 0071 h1=h1f(2:2:4*m4); 0072 k4=3; 0073 g1f=kaiser(4*k4+1,5.3)'.*sinc(-k4:0.5:k4); 0074 g1=g1f(2:2:4*k4); 0075 tamd=[0.00632 0 1.01335; % 'a' UK PPM characteristic [default] 0076 0.00316 0 0.73830; % 'b' DIN PPM 0077 0.10433 0 0; % 'c' VU American 0078 0.13089 0 0; % 'd' VU French 0079 0 0.4 0; % 'e' EBU momentary 0080 0 3 0]; % 'f' EBU short 0081 deft={'','','','','q','q'}; 0082 p0.mm='a'; 0083 p0.ta=0; 0084 p0.tm=0; 0085 p0.td=0; 0086 end 0087 % decode input arguments 0088 if iscell(fsx) 0089 fx=fsx; 0090 [nx,nc]=size(x); 0091 fsp=fx{1}; 0092 fsx=fsp(1); 0093 fy=fsp(2); 0094 nt=fsp(3); 0095 nm=fsp(4); 0096 p=fx{2}; 0097 mop=fx{3}; 0098 else 0099 p=p0; % default values 0100 if nargin<3 || isempty(m) 0101 m='a'; 0102 end 0103 if ischar(m) 0104 p.mm=m; 0105 ix=m(1)-'a'+1; 0106 if ix<=0 || ix>6 0107 ix=1; 0108 end 0109 p.ta=tamd(ix,1); 0110 p.tm=tamd(ix,2); 0111 p.td=tamd(ix,3); 0112 p.mm=[p.mm deft{ix}]; % add in default toggles 0113 elseif isstruct(m) 0114 p=m; 0115 else 0116 p.mm=' '; 0117 p.ta=m; 0118 p.tm=0; 0119 p.td=m; 0120 end 0121 [nx,nc]=size(x); 0122 mop=rem(sum(repmat('zpsowvqu',length(p.mm),1)-p.mm(ones(8,1),:)'==0,1),2); % set option toggles 0123 fy=fsx*(1+3*mop(4)); % effective sample rate 0124 nt=round(fy*0.1); % 0.1 s resolution in samples fro MA filter 0125 nm=round(fy*p.tm/nt); % number of blocks for MA filter 0126 fsp=[fsx fy nt nm]; 0127 fx={fsp p mop zeros(m4,nc),zeros(2*m4-1,nc),zeros(k4,nc),zeros(2*k4-1,nc), zeros(1,nc), 0, zeros(nm+1,nc), zeros(0,nc)}; 0128 end 0129 0130 % Stage 1: oversampling we use two cascaded half-band filters to upsample by 4 0131 0132 if mop(4) 0133 nx2=2*nx; 0134 v=zeros(nx2,nc); 0135 v(1:2:2*m4,:)=fx{4}; 0136 v(2*m4+1:2:nx2,:)=x(1:nx-m4,:); 0137 fx{4}=x(nx-m4+1:nx,:); 0138 [v(2:2:nx2,:),fx{5}]=filter(h1,1,x,fx{5}); % delayed by 2*m4 samples 0139 nx4=4*nx; 0140 y=zeros(nx4,nc); 0141 y(1:2:2*k4,:)=fx{6}; 0142 y(2*k4+1:2:nx4-1,:)=v(1:nx2-k4,:); 0143 fx{6}=v(nx2-k4+1:nx2,:); 0144 [y(2:2:nx4,nc),fx{7}]=filter(g1,1,v,fx{7}); % delayed by 4*m4+2*k4 samples 0145 else 0146 y=x; 0147 end 0148 ty=1/fy; 0149 % Stage 2 make positive 0150 if mop(7) 0151 y=y.^2; 0152 else 0153 y=abs(y); 0154 end 0155 % Stage 2: attack filter 0156 if p.ta>0 && ~mop(3) 0157 za=exp(-ty/p.ta); 0158 [y,fx{8}]=filter(1-za,[1 -za],y,fx{8}); 0159 end 0160 0161 % Stage 3: moving average filter 0162 if nm>0 0163 ny=size(y,1); 0164 y=cumsum(y,1); 0165 nr=fx{9}; 0166 jj=nt-nr:nt:ny; % end of frame indices 0167 kk=length(jj); 0168 if kk>0 0169 nmkk=nm+kk; 0170 v=zeros(nmkk,nc); 0171 v(nm+2:nmkk,:)=y(jj(2:kk),:)-y(jj(1:kk-1),:); % v(:,nc) contains the sum of 0.1 second blocks 0172 v(nm+1,:)=y(jj(1),:)+fx{10}(2,:); 0173 fx{10}(2,:)=y(ny,:)-y(jj(kk),:); 0174 v(2:nm,:)=fx{10}(3:nm+1,:); % saved values 0175 fx{10}(3:nm+1,:)=v(2+kk:nmkk,:); 0176 v=cumsum(v,1); 0177 v(nm+1:nmkk,:)=v(nm+1:nmkk,:)-v(1:kk,:); % perform MA filter 0178 v(nm,:)=fx{10}(1,:); % final MA output from previous chunk 0179 y=v(nm+floor((nr+1:nr+ny)/nt),:)/(nm*nt); % copy MA ouptuts into output buffer 0180 fx{10}(1,:)=v(nmkk,:); % save the final MA filter output for the next chunk 0181 fx{9}=ny-jj(kk); % update the length of the tail 0182 else % no completed 100 ms blocks in this chunk 0183 fx{9}(1)=ny+nr; % update the length of the tail 0184 fx{10}(2,:)=y(ny,:)+fx{10}(2,:); % update the sum of the tail 0185 y=fx{10}(ones(ny,1),:); % set outputs to the final MA output from previous chunk 0186 end 0187 end 0188 % Stage 4: decay filter 0189 if p.td>0 0190 [v,vk,fx{11}]=maxfilt(y,exp(-ty/p.td),Inf,1,fx{11}); 0191 else 0192 v=y; 0193 end 0194 0195 % Stage 5: decimate and ouput 0196 if mop(5) % if 'w' option specified, output y as well 0197 if mop(4) % if oversampling, we need to decimate 0198 y=y(1:4:end,:); 0199 v=v(1:4:end,:); 0200 end 0201 if ~mop(8) 0202 y=(20-10*mop(7))*log10(y); 0203 v=(20-10*mop(7))*log10(v); 0204 end 0205 fx1=fx; 0206 fx=y; 0207 else 0208 if mop(4) % if oversampling, we need to decimate 0209 v=v(1:4:end,:); 0210 end 0211 if ~mop(8) 0212 v=(20-10*mop(7))*log10(v); 0213 end 0214 end 0215 if mop(6) || nargout==0 0216 t=(1:length(x))/fsx; 0217 ax(1)=subplot(211); 0218 plot(t,x,'-b'); 0219 if mop(7)>0 && mop(8)>0 0220 v=sqrt(v); 0221 if mop(5) 0222 y=sqrt(y); 0223 end 0224 end 0225 ax(2)=subplot(212); 0226 if mop(5) % if 'w' option specified, output y as well 0227 plot(t,v,'-b',t,y,'-r'); 0228 else 0229 plot(t,v,'-b'); 0230 end 0231 if ~mop(8) 0232 set(gca,'ylim',[max(min(v)-1,max(v)-40) max(v)+1]); 0233 ylabel('dB'); 0234 elseif mop(7) 0235 ylabel('x_{RMS}'); 0236 else 0237 ylabel('|x|'); 0238 end 0239 linkaxes(ax,'x'); 0240 xlabel(['Time (' xticksi 's)']); 0241 end