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
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 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]=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