Home > voicebox > sigalign.m

sigalign

PURPOSE ^

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

SYNOPSIS ^

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

DESCRIPTION ^

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 martch to the length of rr

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [d,g,rr,ss]=sigalign(s,r,maxd,m,fs)
0002 %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 martch to the length of rr
0023 
0024 
0025 %      Copyright (C) Mike Brookes 2011
0026 %      Version: $Id: sigalign.m 713 2011-10-16 14:45:43Z 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 General Public License as published by
0034 %   the Free Software Foundation; either version 2 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 General Public License for more details.
0041 %
0042 %   You can obtain a copy of the GNU General Public License from
0043 %   http://www.gnu.org/copyleft/gpl.html or by writing to
0044 %   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
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]=stdspectrum(2,'z',fs);
0089         s=filter(b,a,s);
0090         r=filter(b,a,r);
0091     elseif any(m=='b')
0092         [b,a]=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=irfft(rfft([s(sxi:sxj); zeros(fl-sxj+sxi-1,1)]).*conj(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

Generated on Tue 10-Oct-2017 08:30:10 by m2html © 2003