# v_sigalign

## PURPOSE

V_SIGALIGN align a clean reference with a noisy signal [d,g,rr,ss]=(s,r,maxd,m,fs)

## SYNOPSIS

function [d,g,rr,ss]=v_sigalign(s,r,maxd,m,fs)

## DESCRIPTION

    V_SIGALIGN align a clean reference with a noisy signal [d,g,rr,ss]=(s,r,maxd,m,fs)
     Inputs:
                m  mode
                     u = unity gain
                     g = find optimal gain [default]
                     a = A-weight the signals
                     b = weight signals by BS-468
                     s = find delay to maximize the correlation coefficient between r and s [default]
                     S = find delay to maximize the energy of the component of r in s
                     p = plot result
                s  test signal
                r  reference signal
             maxd  [+-max] or [min max] delay allowed in samples or fractions of length(r)
                   default is maximum that ensures at least 50% of r or s in the overlap
               fs  sample frequency (only used for filtering and plotting)

     Outputs:
                d = optimum delay to apply to r
                g = optimal gain to apply to r
               rr = g*r(* -d)  [zero padded to match s if ss output is not given]
               ss = s truncated if necessary to match the length of rr

## CROSS-REFERENCE INFORMATION

This function calls:

- v_irfft V_IRFFT Inverse fft of a conjugate symmetric spectrum X=(Y,N,D)
- v_rfft V_RFFT Calculate the DFT of real data Y=(X,N,D)
- 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)

This function is called by:

## SOURCE CODE

    0001 function [d,g,rr,ss]=v_sigalign(s,r,maxd,m,fs)
    0002 %V_SIGALIGN align a clean reference with a noisy signal [d,g,rr,ss]=(s,r,maxd,m,fs)
    0003 % Inputs:
    0004 %            m  mode
    0005 %                 u = unity gain
    0006 %                 g = find optimal gain [default]
    0007 %                 a = A-weight the signals
    0008 %                 b = weight signals by BS-468
    0009 %                 s = find delay to maximize the correlation coefficient between r and s [default]
    0010 %                 S = find delay to maximize the energy of the component of r in s
    0011 %                 p = plot result
    0012 %            s  test signal
    0013 %            r  reference signal
    0014 %         maxd  [+-max] or [min max] delay allowed in samples or fractions of length(r)
    0015 %               default is maximum that ensures at least 50% of r or s in the overlap
    0016 %           fs  sample frequency (only used for filtering and plotting)
    0017 %
    0018 % Outputs:
    0019 %            d = optimum delay to apply to r
    0020 %            g = optimal gain to apply to r
    0021 %           rr = g*r(* -d)  [zero padded to match s if ss output is not given]
    0022 %           ss = s truncated if necessary to match the length of rr
    0023
    0024
    0025 %      Copyright (C) Mike Brookes 2011
    0026 %      Version: $Id: v_sigalign.m 10865 2018-09-21 17:22:45Z dmb $
    0027 %
    0028 %   VOICEBOX is a MATLAB toolbox for speech processing.
    0029 %   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
    0030 %
    0031 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0032 %   This program is free software; you can redistribute it and/or modify
    0033 %   it under the terms of the GNU Lesser General Public License as published by
    0034 %   the Free Software Foundation; either version 3 of the License, or
    0035 %   (at your option) any later version.
    0036 %
    0037 %   This program is distributed in the hope that it will be useful,
    0038 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
    0039 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    0040 %   GNU Lesser General Public License for more details.
    0041 %
    0042 %   You can obtain a copy of the GNU Lesser General Public License from
    0043 %   https://www.gnu.org/licenses/ .
    0044 %
    0045 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0046
    0047 % Bugs/Suggestions
    0048 % 1. add option to calculate a DC offset
    0049 % 2. optionally find optimal fractional time shift
    0050 % 3. split long signals into chunks to reduce memory requirements
    0051
    0052 ns=length(s);
    0053 nr=length(r);
    0054 if numel(s)~=ns || numel(r)~=nr
    0055     error('Inputs cannot be matrices');
    0056 end
    0057 s=s(:);
    0058 r=r(:);
    0059 if nargin<3
    0060     maxd=[];
    0061 end
    0062 switch numel(maxd)
    0063     case 0
    0064         if nr<ns
    0065             lmm=[-0.25*nr ns-0.75*nr];
    0066         else
    0067             lmm=[-0.25*ns nr-0.75*ns];
    0068         end
    0069     case 1
    0070         lmm=[-maxd maxd];
    0071     otherwise
    0072         lmm=maxd(1:2);
    0073 end
    0074 lmm=round(lmm.*(1+(nr-1)*(abs(lmm)<1)));  % convert fractions of nr to samples
    0075 lmin=lmm(1);
    0076 lmax=lmm(2);
    0077 lags=lmax-lmin+1;
    0078 if lags<=0
    0079     error('Invalid lag limits');
    0080 end
    0081 if nargin<4 || ~numel(m)
    0082     m='gs';
    0083 end
    0084 if nargin<5 || ~numel(fs)
    0085     fs=[];
    0086 else
    0087     if any(m=='a')
    0088         [b,a]=v_stdspectrum(2,'z',fs);
    0089         s=filter(b,a,s);
    0090         r=filter(b,a,r);
    0091     elseif any(m=='b')
    0092         [b,a]=v_stdspectrum(8,'z',fs);
    0093         s=filter(b,a,s);
    0094         r=filter(b,a,r);
    0095     end
    0096 end
    0097
    0098 % now do cross correlation
    0099
    0100 rxi=max(1,1-lmin);   % first reference sample needed
    0101 rxj=min(nr,ns-lmax); % last reference sample needed
    0102 nrx=rxj-rxi+1;          % length of reference segment
    0103 if nrx<1
    0104     error('Reference signal too short');
    0105 end
    0106 fl=2^nextpow2(lmax-lmin+nrx);
    0107 sxi=max(1,rxi+lmin); % first signal sample needed
    0108 sxj=min(ns,rxj+lmax); % last signal sample needed
    0109 rs=v_irfft(v_rfft([s(sxi:sxj); zeros(fl-sxj+sxi-1,1)]).*conj(v_rfft([r(rxi:rxj); zeros(fl-rxj+rxi-1,1)])));
    0110 rsu=rs(1:lags);
    0111 ssq=cumsum(s(sxi:sxj).^2);
    0112 ssqd=[ssq(nrx); ssq(nrx+1:lmax-lmin+nrx)-ssq(1:lmax-lmin)];
    0113 if any (m=='S') % maximize energy of common component
    0114     [cmx,icx]=max(abs(rsu)); % maximize cross correlation
    0115 else
    0116     [cmx,icx]=max(rsu.^2./ssqd); % maximize correlation coefficient
    0117 end
    0118 d=icx-1+lmin;
    0119 ia=max(1,d+1); % first sample of s in common region
    0120 ja=min(ns,d+nr); % last sample of s in common region
    0121 ija=ia:ja;
    0122 ijad=ija-d;
    0123 rr=r(ijad);
    0124 ss=s(ija);
    0125 if any (m=='u')
    0126     g=1;
    0127 else
    0128 g=sum(rr.*ss)/sum(rr.^2);   % gain to apply to r
    0129 end
    0130 rr=rr*g;
    0131 if ~nargout || any(m=='p')
    0132     xco=sum(rr.*ss)/sqrt(sum(rr.^2)*sum(ss.^2));
    0133     snr=sum(rr.^2)/sum((rr-ss).^2);
    0134     if numel(fs)==1
    0135         tun='s';
    0136     else
    0137         tun='samples';
    0138         fs=1;
    0139     end
    0140     subplot(311);
    0141     plot(ija/fs,rr);
    0142     pm='+-';
    0143     title(sprintf('Ref delay = %.2g %s, %cGain = %.2g dB, Xcorr = %.2g, SNR = %.2g dB',d/fs,tun,pm(1+(g<0)),20*log10(g),xco,10*log10(snr)));
    0144     ylabel('Reference');
    0145     set(gca,'XLim',ija([1 end])/fs);
    0146     axh(2)=gca;
    0147     subplot(312);
    0148     plot(ija/fs,ss);
    0149     ylabel('Signal');
    0150     set(gca,'XLim',ija([1 end])/fs);
    0151     axh(1)=gca;
    0152     subplot(313);
    0153     plot(ija/fs,ss-rr);
    0154     ylabel('Residual');
    0155     xlabel(sprintf('Time (%s)',tun));
    0156     set(gca,'XLim',ija([1 end])/fs);
    0157     axh(3)=gca;
    0158     linkaxes(axh(1:3),'x');
    0159 end
    0160 if nargout==3
    0161     rr=[zeros(ia-1,1); rr; zeros(ns-ja,1)]; % force to be the size of s
    0162 end
    0163
    0164
    0165