findtune

FINDTUNE   get the tune value from turn by turn positions

TUNE=FINDTUNE(POS,METHOD)

POS:       Tune-by-turn particle position
METHOD:    Method for tune determination:
               1: Highest peak in fft
               2: Interpolation on fft results
               3: Windowing + interpolation

[TUNE,SPECTRUM]=FINDTUNE(...) Also returns the fft

0001 function [tune,spectrum]=findtune(pos,method)
0002 %FINDTUNE   get the tune value from turn by turn positions
0003 %
0004 %TUNE=FINDTUNE(POS,METHOD)
0005 %
0006 %POS:       Tune-by-turn particle position
0007 %METHOD:    Method for tune determination:
0008 %               1: Highest peak in fft
0009 %               2: Interpolation on fft results
0010 %               3: Windowing + interpolation
0011 %
0012 %[TUNE,SPECTRUM]=FINDTUNE(...) Also returns the fft
0013 
0014 
0015 
0016 if nargin < 2, method=3; end
0017 
0018 nturns=size(pos,1);
0019 nparts=size(pos,2);
0020 nt2=fix(nturns/2);
0021 
0022 posm=mean(pos);
0023 wrong=~isfinite(posm);
0024 
0025 switch method
0026 case 1
0027 methname='highest peak';
0028 pos2=pos-posm(ones(nturns,1),:);
0029 spectrum=fft(pos2);
0030 [vmax,rmax]=max(abs(spectrum(1:nt2,:))); %#ok<ASGLU>
0031 tune=(rmax-1)/nturns;
0032 
0033 case 2
0034 methname='interpolation';
0035 pos2=pos-posm(ones(nturns,1),:);
0036 spectrum=fft(pos2);
0037 [vmax,rmax]=max(abs(spectrum(1:nt2,:))); %#ok<ASGLU>
0038 rmax=rmax+(rmax==1);
0039 kmax=sub2ind([nturns nparts],rmax,1:nparts);
0040 back=(spectrum(kmax-1) > spectrum(kmax+1));
0041 k1=kmax-back;
0042 k2=k1+1;
0043 v1=abs(spectrum(k1));
0044 v2=abs(spectrum(k2));
0045 tune=(rmax-back-1 +(v2./(v1+v2)))/nturns;
0046 
0047 case 3
0048 methname='window + interp.';
0049 w=hann_window(nturns);
0050 pos2=(pos-posm(ones(nturns,1),:)).*w(:,ones(1,nparts));
0051 spectrum=fft(pos2);
0052 [vmax,rmax]=max(abs(spectrum(1:nt2,:))); %#ok<ASGLU>
0053 rmax=rmax+(rmax==1);
0054 kmax=sub2ind([nturns nparts],rmax,1:nparts);
0055 back=(spectrum(kmax-1) > spectrum(kmax+1));
0056 k1=kmax-back;
0057 k2=k1+1;
0058 v1=abs(spectrum(k1));
0059 v2=abs(spectrum(k2));
0060 tune=(rmax-back-1 +((2*v2-v1)./(v1+v2)))/nturns;
0061 %tune2=(rmax-back-1)/nturns + asin(phi(v1,v2,cos(2*pi/nturns))*sin(2*pi/nturns))/2/pi;
0062 %disp(['method 4 tune: ' num2str(mean2(tune2')) ' (rms: ' num2str(std(tune2')) ')']);
0063 
0064 end
0065 tune(wrong)=NaN;
0066 errmax=2.5*std(tune');
0067 keep=(abs(tune-mean(tune'))<=errmax);
0068 reject=find(~(keep | wrong));
0069 for bpm=reject
0070     [bname,kdx]=sr.bpmname(bpm);
0071     disp(['reject ' bname ' (' num2str(kdx) ')']);
0072 end
0073 fprintf('%20s tune:%g (rms:%g)\n',methname, mean(tune(keep)'),std(tune(keep)'));
0074 
0075 function vv=phi(a,b,c) %#ok<DEFNU>
0076 d1=c*(a+b);
0077 delt=d1.*d1 - 2*a.*b.*(2*c*c-c-1).*a.*a - b.*b - 2*a.*b*c;
0078 vv=(-(a+b*c).*(a-b) + b.*sqrt(delt))./(a.*a + b.*b  +2*a.*b*c);
0079 
0080 
0081 function w=hann_window(n)
0082 w=0.5*(1-cos(2*pi*(0:n-1)'/n));