# v_ppmvu

## PURPOSE

V_PPMVU calculate PPM, VU or EBU level of an audio signal [V,FX,FX1]=(X,FSX,M)

## SYNOPSIS

function [v,fx,fx1]=v_ppmvu(x,fsx,m)

## DESCRIPTION

    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)

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

## SOURCE CODE

    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 Lesser General Public License as published by
    0052 %   the Free Software Foundation; either version 3 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 Lesser General Public License for more details.
    0059 %
    0060 %   You can obtain a copy of the GNU Lesser General Public License from
    0061 %   https://www.gnu.org/licenses/ .
    0062 %
    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