# v_snrseg

## PURPOSE

V_SNRSEG Measure segmental and global SNR [SEG,GLO]=(S,R,FS,M,TF)

## SYNOPSIS

function [seg,glo,tc,snf,vf]=v_snrseg(s,r,fs,m,tf)

## DESCRIPTION

    V_SNRSEG Measure segmental and global SNR [SEG,GLO]=(S,R,FS,M,TF)

    Usage: (1) seg=v_snrseg(s,r,fs);                  % s & r are noisy and clean signal
           (2) seg=v_snrseg(s,r,fs,'wz');             % no VAD or inerpolation used ['Vq' is default]
           (3) [seg,glo]=v_snrseg(s,r,fs,'Vq',0.03);  % 30 ms frames

     Inputs:    s  test signal
                r  reference signal
               fs  sample frequency (Hz)
                m  mode [default = 'Vq']
                     w = No VAD - use whole file
                     v = use sohn VAD to discard silent portions
                     V = use P.56-based VAD to discard silent portions [default]
                     a = A-weight the signals
                     b = weight signals by BS-468
                     q = use linear interpolation to remove delays +- 1 sample [default]
                     z = do not do any alignment
                     p = plot results
               tf  frame increment [0.01]

     Outputs: seg = Segmental SNR in dB
              glo = Global SNR in dB (typically 7 dB greater than SNR-seg)
              tc  = time at centre of each frame (seconds)
              snf = Segmental SNR in dB in each frame
              vf  = Boolean mask indicating frames that valid (from VAD)

     This function compares a noisy signal, S, with a clean reference, R, and
     computes the segemntal signal-to-noise ratio (SNR) in dB. The signals,
     which must be of the same length, are split into non-overlapping frames
     of length TF (default 10 ms) and the SNR of each frame in dB is calculated.
     The segmental SNR is the average of these values, i.e.
             SEG = mean(10*log10(sum(Ri^2)/sum((Si-Ri)^2))
     where the mean is over frames and the sum runs over one particular frame.
     Two optional modifications can be made to this basic formula:

        (a) Frames are excluded if there is no significant energy in the R
            signal. The idea is to limit the calculation to frames in which
            speech is active. By default, the v_voicebox function "v_activlev" is
            used to detect the inactive frames (the 'V' mode option).

        (b) In each frame independently, the reference signal is shifted by up
            to +- 1 sample to find the alignment than minimizes the noise
            component (S-R)^2. This shifting accounts for small misalignments
            and/or sample frequency differences between the two signals. For
            larger shifts, you can use the v_voicebox function "v_sigalign".
            Accurate alignemnt is especially important at high SNR values.

     If no M argument is specified, both these modifications will be applied;
     this is equivalent to specifying M='Vq'.

## CROSS-REFERENCE INFORMATION

This function calls:

- v_activlev V_ACTIVLEV Measure active speech level as in ITU-T P.56 [LEV,AF,FSO]=(sp,FS,MODE)
- v_stdspectrum V_STDSPECTRUM Generate standard acoustic/speech spectra in s- or z-domain [B,A,SI,SN]=(S,M,F,N,ZI,BS,AS)
- v_vadsohn V_VADSOHN implements a voice activity detector [VS,ZO]=(S,FSZ,M,P)

This function is called by:

## SOURCE CODE

    0001 function [seg,glo,tc,snf,vf]=v_snrseg(s,r,fs,m,tf)
    0002 %V_SNRSEG Measure segmental and global SNR [SEG,GLO]=(S,R,FS,M,TF)
    0003 %
    0004 %Usage: (1) seg=v_snrseg(s,r,fs);                  % s & r are noisy and clean signal
    0005 %       (2) seg=v_snrseg(s,r,fs,'wz');             % no VAD or inerpolation used ['Vq' is default]
    0006 %       (3) [seg,glo]=v_snrseg(s,r,fs,'Vq',0.03);  % 30 ms frames
    0007 %
    0008 % Inputs:    s  test signal
    0009 %            r  reference signal
    0010 %           fs  sample frequency (Hz)
    0011 %            m  mode [default = 'Vq']
    0012 %                 w = No VAD - use whole file
    0013 %                 v = use sohn VAD to discard silent portions
    0014 %                 V = use P.56-based VAD to discard silent portions [default]
    0015 %                 a = A-weight the signals
    0016 %                 b = weight signals by BS-468
    0017 %                 q = use linear interpolation to remove delays +- 1 sample [default]
    0018 %                 z = do not do any alignment
    0019 %                 p = plot results
    0020 %           tf  frame increment [0.01]
    0021 %
    0022 % Outputs: seg = Segmental SNR in dB
    0023 %          glo = Global SNR in dB (typically 7 dB greater than SNR-seg)
    0024 %          tc  = time at centre of each frame (seconds)
    0025 %          snf = Segmental SNR in dB in each frame
    0026 %          vf  = Boolean mask indicating frames that valid (from VAD)
    0027 %
    0028 % This function compares a noisy signal, S, with a clean reference, R, and
    0029 % computes the segemntal signal-to-noise ratio (SNR) in dB. The signals,
    0030 % which must be of the same length, are split into non-overlapping frames
    0031 % of length TF (default 10 ms) and the SNR of each frame in dB is calculated.
    0032 % The segmental SNR is the average of these values, i.e.
    0033 %         SEG = mean(10*log10(sum(Ri^2)/sum((Si-Ri)^2))
    0034 % where the mean is over frames and the sum runs over one particular frame.
    0035 % Two optional modifications can be made to this basic formula:
    0036 %
    0037 %    (a) Frames are excluded if there is no significant energy in the R
    0038 %        signal. The idea is to limit the calculation to frames in which
    0039 %        speech is active. By default, the v_voicebox function "v_activlev" is
    0040 %        used to detect the inactive frames (the 'V' mode option).
    0041 %
    0042 %    (b) In each frame independently, the reference signal is shifted by up
    0043 %        to +- 1 sample to find the alignment than minimizes the noise
    0044 %        component (S-R)^2. This shifting accounts for small misalignments
    0045 %        and/or sample frequency differences between the two signals. For
    0046 %        larger shifts, you can use the v_voicebox function "v_sigalign".
    0047 %        Accurate alignemnt is especially important at high SNR values.
    0048 %
    0049 % If no M argument is specified, both these modifications will be applied;
    0050 % this is equivalent to specifying M='Vq'.
    0051
    0052 % Bugs/suggestions
    0053 % (1) Optionally restrict the bandwidth to the smaller of the two
    0054 %     bandwidths either with an extra parameter or automatically determined
    0055 % (2) Optionally apply a gain to the s signal to maximize the SNR and
    0056 %     output the gain as an additional parameter
    0057
    0058 %      Copyright (C) Mike Brookes 2011
    0059 %      Version: $Id: v_snrseg.m 10865 2018-09-21 17:22:45Z dmb $
    0060 %
    0061 %   VOICEBOX is a MATLAB toolbox for speech processing.
    0062 %   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
    0063 %
    0064 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0065 %   This program is free software; you can redistribute it and/or modify
    0066 %   it under the terms of the GNU Lesser General Public License as published by
    0067 %   the Free Software Foundation; either version 3 of the License, or
    0068 %   (at your option) any later version.
    0069 %
    0070 %   This program is distributed in the hope that it will be useful,
    0071 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
    0072 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    0073 %   GNU Lesser General Public License for more details.
    0074 %
    0075 %   You can obtain a copy of the GNU Lesser General Public License from
    0076 %   https://www.gnu.org/licenses/ .
    0077 %
    0078 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0079
    0080 if nargin<4 || ~ischar(m)
    0081     m='Vq';
    0082 end
    0083 if nargin<5 || ~numel(tf)
    0084     tf=0.01;                                            % default frame length is 10 ms
    0085 end
    0086 snmax=100;                                              % clipping limit for SNR
    0087
    0088 % filter the input signals if required
    0089
    0090 if any(m=='a')                                          % A-weighting
    0091     [b,a]=v_stdspectrum(2,'z',fs);                      % create bandpass filter
    0092     s=filter(b,a,s);                                    % filter test signal
    0093     r=filter(b,a,r);                                    % filter reference signal
    0094 elseif any(m=='b')                                      %  BS-468 weighting
    0095     [b,a]=v_stdspectrum(8,'z',fs);                      % create bandpass filter
    0096     s=filter(b,a,s);                                    % filter test signal
    0097     r=filter(b,a,r);                                    % filter reference signal
    0098 end
    0099
    0100 mq=~any(m=='z');                                        % flag indicating we need to do alignment
    0101 nr=min(length(r), length(s));                           % signal length to align
    0102 kf=round(tf*fs);                                        % length of frame in samples
    0103 ifr=kf+mq:kf:nr-mq;                                     % ending sample of each frame (allowing +-mq at each end)
    0104 ifl=ifr(end);                                           % end sample of last frame
    0105 nf=numel(ifr);                                          % number of frames
    0106 if mq                                                   % perform linear interpolation
    0107     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0108     % Consider the linearly interpolated signal sm(i)=(1-am)*s(i)+am*s(i+1) where 0<=am<=1.     %
    0109     % For each frame we calculate the value of am that minimizes the error sum((sm(i)-r(i))^2). %
    0110     % The solution is am = sum((s(i)-r(i)).*(s(i)-s(i+1))) / sum((s(i)-s(i+1))^2).              %
    0111     % Similarly, we calculate the optimum ap in sp(i)=(1-ap)*s(i)+ap*s(i-1) and then            %
    0112     % finally, we select whichever of sm(i) and sp(i) with the lowest error.                    %
    0113     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0114     ssm=reshape(s(2:ifl)-s(3:ifl+1),kf,nf);
    0115     ssp=reshape(s(2:ifl)-s(1:ifl-1),kf,nf);
    0116     sr=reshape(s(2:ifl)-r(2:ifl),kf,nf);                % enframed error: test-ref (one frame per column)
    0117     am=min(max(sum(sr.*ssm,1)./sum(ssm.^2,1),0),1);     % optimum dist between s(2:ifl) and s(3:ifl+1)
    0118     ap=min(max(sum(sr.*ssp,1)./sum(ssp.^2,1),0),1);     % optimum dist between s(2:ifl) and s(1:ifl-1)
    0119     ef=min(sum((sr-repmat(am,kf,1).*ssm).^2,1),sum((sr-repmat(ap,kf,1).*ssp).^2,1)); % select the best for each frame
    0120 else                                                    % no interpolation
    0121     ef=sum(reshape((s(1:ifl)-r(1:ifl)).^2,kf,nf),1);    % sum of squared errors in each frame
    0122 end
    0123 rf=sum(reshape(r(mq+1:ifl).^2,kf,nf),1);                % sum of sqared reference signal in each frame
    0124 em=ef==0;                                               % mask for zero noise frames
    0125 rm=rf==0;                                               % mask for zero reference frames
    0126 snf=10*log10((rf+rm)./(ef+em));
    0127 snf(rm)=-snmax;
    0128 snf(em)=snmax;
    0129
    0130 % select the frames to include
    0131
    0132 if any(m=='w')                                          % 'w' option: do not discard silent frames
    0133     vf=true(1,nf);                                      % include all frames
    0134 elseif any(m=='v')                                      % 'v' option: use Sohn VAD to discard silent frames
    0135     vs=v_vadsohn(r,fs,'na');
    0136     nvs=length(vs);
    0137     [vss,vix]=sort([ifr'; vs(:,2)]);
    0138     vjx=zeros(nvs+nf,5);
    0139     vjx(vix,1)=(1:nvs+nf)';                             % sorted position
    0140     vjx(1:nf,2)=vjx(1:nf,1)-(1:nf)';                    % prev VAD frame end (or 0 or nvs+1 if none)
    0141     vjx(nf+1:end,2)=vjx(nf+1:end,1)-(1:nvs)';           % prev snr frame end (or 0 or nvs+1 if none)
    0142     dvs=[vss(1)-mq; vss(2:end)-vss(1:end-1)];           % number of samples from previous frame boundary
    0143     vjx(:,3)=dvs(vjx(:,1));                             % number of samples from previous frame boundary
    0144     vjx(1:nf,4)=vs(min(1+vjx(1:nf,2),nvs),3);           % VAD result for samples between prev frame boundary and this one
    0145     vjx(nf+1:end,4)=vs(:,3);                            % VAD result for samples between prev frame boundary and this one
    0146     vjx(1:nf,5)=1:nf;                                   % SNR frame to accumulate into
    0147     vjx(vjx(nf+1:end,2)>=nf,3)=0;                       % zap any VAD frame beyond the last snr frame
    0148     vjx(nf+1:end,5)=min(vjx(nf+1:end,2)+1,nf);          % SNR frame to accumulate into
    0149     vf=full(sparse(1,vjx(:,5),vjx(:,3).*vjx(:,4),1,nf))>kf/2; % accumulate into SNR frames and compare with threshold
    0150 else                                                    % default is 'V': use P.56 VAD to discard silent frames
    0151     [lev,af,fso,vad]=v_activlev(r,fs);                  % do VAD on reference signal
    0152     vf=sum(reshape(vad(mq+1:ifl),kf,nf),1)>kf/2;        % find frames that are mostly active
    0153     tc=((1:nf)*kf+(1-kf)/2)/fs;                         % time at centre of frame
    0154 end
    0155 seg=mean(snf(vf));
    0156 glo=10*log10(sum(rf(vf))/sum(ef(vf)));
    0157 if ~nargout || any (m=='p')
    0158     subplot(312);                                       % plot reference signal
    0159     plot((1:length(r))/fs,r);
    0160     ylabel('Reference');
    0161     axh(2)=gca;
    0162     subplot(313);                                       % plot SNR + global and segmental means
    0163     snv=snf;
    0164     snv(~vf)=NaN;                                       % plot valid frames in blue
    0165     snu=snf;
    0166     snu(vf>0)=NaN;                                      % plot invalid frames in red
    0167     plot([1 nr]/fs,[glo seg; glo seg],':k',tc,snv,'-b',tc,snu,'-r');
    0168     ylabel('Frame SNR');
    0169     xlabel('Time (s)');
    0170     set(gca,'ylim',min([snv(:);snu(:)])+[-3 60]);       % restrict displayed SNR range
    0171     axh(3)=gca;
    0172     subplot(311);                                       % plot test signal
    0173     plot((1:length(s))/fs,s);
    0174     ylabel('Signal');
    0175     title(sprintf('SNR = %.1f dB, SNR_{seg} = %.1f dB',glo,seg));
    0176     axh(1)=gca;
    0177     linkaxes(axh,'x');                                  % link time axes
    0178     linkaxes(axh(1:2),'y');                             % link signal amplitude axes
    0179     set(gca,'xlim',[1 nr]/fs);
    0180 end