findtunetom2.m 2.94 KB
Newer Older
Thomas Perron's avatar
Thomas Perron committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108
function [tunef,numdc,peakss,phasef,camp,uvect]=findtunetom2(M,divider,window)

%[tunef,numdc]=findtunetom2(M,divider)
%finds thew tune for real or complex sets of orbits. ALl orbits should have
%the same length. one orbit is one collumn, and there ashould be as many
%collums
%as there are different orbits to treat.
% the size of M should be orbit_lengthxnumber_of _orbits. 
% WINDOW should be equal to 1 if a hanex window should be aplied.

%*******************************************************************
%THE VALUE OF THE TUNE IS GIVEN CONSIDERED THAT THE ORBITS ARE RED ON A
%SINGLE BPM (SAMPELING AT REVOLUTION FREQUENCY). IF THE ORBIT IS READ ON
%MORE THAN ONE BPM THE GIVEN TUNE VALUE SHOULD BE MULTIPLIED BY THE NUMBER
%OF BPMS CONSIDERED.
%EG: ALL 224 BPM ON A FIRST TURN: MULTIPLY THE COMPUTED TUNE
%VALUE BY 224
%*******************************************************************


% tune is then a line vector containing the tunes. The tune is calculated
% doing a hanning window and making an interpolation to find the right peak
% value on the spectrum.
%M is the orbit matrix, DIVIDER is the step to resolve the main peak (Hz). step=
%FFTstep/DIVIDER
phases=NaN;
tune=NaN;
pointnotok=~isfinite(M);
size(M);
indbaddata=sum(pointnotok,1)>size(M,1)*0.15;
M(pointnotok)=0;
numdc=0;
c=size(M,2);
n=size(M,1);
nu=[0:1:n-1];
han=(0.53836-0.46164*cos(2*pi*nu./(n-1)))';
hanext=han(:,ones(1,c));
if nargin>2
    if window==1
M=M.*hanext;   
    end
end

res=abs(fft(M(:,~indbaddata)));
Mok=M(:,~indbaddata);
res([1:3 end-3:end],:)=0;
matexp=(i*2*pi*nu./n)';

[bid,ind]=max(res);
if isreal(M) res(round(n/2):n,:)=0; end

for k=1:size(res,2)
    serie=res(:,k);
          
    if ind(k)==1
     subind=[0:1/(divider-1):1];
   col=Mok(:,k);
   
   matcomp=sum(exp(matexp(:,ones(1,divider)).*subind(ones(n,1),:)).*col(:,ones(1,divider)));
    prodscal=abs((matcomp));
    p=phase((matcomp));  
    [bid,indn]=sort(prodscal,2,'descend');
   peakss(k)=prodscal(indn(1)); 
   tune(k)=(subind(indn(1)))/n;
   phases(k)=p(indn(1));
   camp(k)=matcomp(indn(1))/n;
   size(matexp(:,ones(1,divider)))
   size(subind(ones(n,1),:))
   
   uvect(:,k)=exp(matexp.*subind(indn(1)));
   %exp(matexp(:,ones(1,divider).*subind(ones(n,1),:)));
   %peak(k)=0
    %    tune(k)=0;
     %  numdc=numdc+1;
    %elseif ind(k)==n
    %tune(k)=1-1/n;
    
   
    
    else 
    a=ind(k);
    col=Mok(:,k);
    
    %subind=[ind(k)-1.5:1/(divider-1):ind(k)-0.5];
   subind=[ind(k)-1.5:1/(divider-1):ind(k)-0.5];
   
   
   
   matcomp=sum(exp(matexp(:,ones(1,divider)).*subind(ones(n,1),:)).*col(:,ones(1,divider)));
   p=phase((matcomp));  
   prodscal=abs((matcomp));
    [bid,indn]=sort(prodscal,2,'descend');
  
    peakss(k)=prodscal(indn(1)); 
   tune(k)=subind(indn(1))/(n);
    phases(k)=p(indn(1));
     camp(k)=matcomp((indn(1)))/(n);  
   
     uvect(:,k)=exp(matexp.*subind(indn(1)));
    end
  
    
end
    
tunef(indbaddata)=NaN;
   tunef(~indbaddata)=tune;
phasef(indbaddata)=NaN;
phasef(~indbaddata)=phases;