1 function [delayStructNew] = align(xf, yf, delayStruct, i, trueDelay); 2 3 %%%%%%% 4 % Bastiaan's algorithm copied 5 %%%%%%% 6 Ap500 = [1.00, -4.95, 9.801, -9.70299, 4.80298005, -0.9509900499]; 7 Bp500 = [0.662743088639636, -2.5841655608125, 3.77668102146288, -2.45182477425154, 0.596566274575251, 0.0]; 8 Ap200 = [1.00, -4.875, 9.50625, -9.26859375, 4.518439453125, -0.881095693359375]; 9 Bp200 = [0.862545460994275, -3.2832804496114, 4.67892032308828, -2.95798023879133, 0.699796870041299, 0.0]; 10 11 oldMethod = 1; % Turn on or off the old method. The new one is Bastiaan's August 2008 updates 12 THReSHoLD = 2.0; % ADJUSTABLE threshold factor; 4.0 seems good 13 %%%%%%%%%%%%%%%%%%% 14 % use log domain (showed improved performance) 15 xxf = sqrt(real(xf.*conj(xf))+1e-20); 16 yyf = sqrt(real(yf.*conj(yf))+1e-20); 17 delayStruct.sxAll2(:,i) = 20*log10(xxf); 18 delayStruct.syAll2(:,i) = 20*log10(yyf); 19 20 mD = min(i-1,delayStruct.maxDelayb); 21 if oldMethod 22 factor = 1.0; 23 histLenb = 250; 24 xthreshold = factor*median(delayStruct.sxAll2(:,i-mD:i),2); 25 ythreshold = factor*median(delayStruct.syAll2(:,i-mD:i),2); 26 else 27 xthreshold = sum(delayStruct.sxAll2(:,i-mD:i),2)/(delayStruct.maxDelayb+1); 28 29 [yout, delayStruct.z200] = filter(Bp200, Ap200, delayStruct.syAll2(:,i), delayStruct.z200, 2); 30 yout = yout/(delayStruct.maxDelayb+1); 31 ythreshold = mean(delayStruct.syAll2(:,i-mD:i),2); 32 ythreshold = yout; 33 end 34 35 delayStruct.bxspectrum(i) = getBspectrum(delayStruct.sxAll2(:,i), xthreshold, delayStruct.bandfirst, delayStruct.bandlast); 36 delayStruct.byspectrum(i) = getBspectrum(delayStruct.syAll2(:,i), ythreshold, delayStruct.bandfirst, delayStruct.bandlast); 37 38 delayStruct.bxhist(end-mD:end) = delayStruct.bxspectrum(i-mD:i); 39 40 delayStruct.bcount(:,i) = hisser2(delayStruct.byspectrum(i), flipud(delayStruct.bxhist), delayStruct.bandfirst, delayStruct.bandlast); 41 [delayStruct.fout(:,i), delayStruct.z500] = filter(Bp500, Ap500, delayStruct.bcount(:,i), delayStruct.z500, 2); 42 if oldMethod 43 %delayStruct.new(:,i) = sum(delayStruct.bcount(:,max(1,i-histLenb+1):i),2); % using the history range 44 tmpVec = [delayStruct.fout(1,i)*ones(2,1); delayStruct.fout(:,i); delayStruct.fout(end,i)*ones(2,1)]; % using the history range 45 tmpVec = filter(ones(1,5), 1, tmpVec); 46 delayStruct.new(:,i) = tmpVec(5:end); 47 %delayStruct.new(:,i) = delayStruct.fout(:,i); % using the history range 48 else 49 [delayStruct.fout(:,i), delayStruct.z500] = filter(Bp500, Ap500, delayStruct.bcount(:,i), delayStruct.z500, 2); 50 % NEW CODE 51 delayStruct.new(:,i) = filter([-1,-2,1,4,1,-2,-1], 1, delayStruct.fout(:,i)); %remv smth component 52 delayStruct.new(1:end-3,i) = delayStruct.new(1+3:end,i); 53 delayStruct.new(1:6,i) = 0.0; 54 delayStruct.new(end-6:end,i) = 0.0; % ends are no good 55 end 56 [valuen, tempdelay] = min(delayStruct.new(:,i)); % find minimum 57 if oldMethod 58 threshold = valuen + (max(delayStruct.new(:,i)) - valuen)/4; 59 thIndex = find(delayStruct.new(:,i) <= threshold); 60 if (i > 1) 61 delayDiff = abs(delayStruct.delay(i-1)-tempdelay+1); 62 if (delayStruct.oneGoodEstimate & (max(diff(thIndex)) > 1) & (delayDiff < 10)) 63 % We consider this minimum to be significant, hence update the delay 64 delayStruct.delay(i) = tempdelay; 65 elseif (~delayStruct.oneGoodEstimate & (max(diff(thIndex)) > 1)) 66 delayStruct.delay(i) = tempdelay; 67 if (i > histLenb) 68 delayStruct.oneGoodEstimate = 1; 69 end 70 else 71 delayStruct.delay(i) = delayStruct.delay(i-1); 72 end 73 else 74 delayStruct.delay(i) = tempdelay; 75 end 76 else 77 threshold = THReSHoLD*std(delayStruct.new(:,i)); % set updata threshold 78 if ((-valuen > threshold) | (i < delayStruct.smlength)) % see if you want to update delay 79 delayStruct.delay(i) = tempdelay; 80 else 81 delayStruct.delay(i) = delayStruct.delay(i-1); 82 end 83 % END NEW CODE 84 end 85 delayStructNew = delayStruct; 86 87 % administrative and plotting stuff 88 if( 0) 89 figure(10); 90 plot([1:length(delayStructNew.new(:,i))],delayStructNew.new(:,i),trueDelay*[1 1],[min(delayStructNew.new(:,i)),max(delayStructNew.new(:,i))],'r',[1 length(delayStructNew.new(:,i))],threshold*[1 1],'r:', 'LineWidth',2); 91 %plot([1:length(delayStructNew.bcount(:,i))],delayStructNew.bcount(:,i),trueDelay*[1 1],[min(delayStructNew.bcount(:,i)),max(delayStructNew.bcount(:,i))],'r','LineWidth',2); 92 %plot([thedelay,thedelay],[min(fcount(:,i)),max(fcount(:,i))],'r'); 93 %title(sprintf('bin count and known delay at time %5.1f s\n',(i-1)*(support/(fs*oversampling)))); 94 title(delayStructNew.oneGoodEstimate) 95 xlabel('delay in frames'); 96 %hold off; 97 drawnow 98 end 99
