Changeset 643 in ETALON for SPESO


Ignore:
Timestamp:
Mar 31, 2017, 5:22:59 PM (7 years ago)
Author:
hodnevuc
Message:
 
Location:
SPESO/ana2015
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • SPESO/ana2015/analyse_channel.m

    r634 r643  
    11function channelInfo=analyse_channel(fnamedir,fnamebase,channelLoop, channelParameters,plotShot,plots_dir)
     2
    23global basedir;
    34global imageSizes;
    … …  
    89channelInfo=[];
    910channelContent='';
     11
    1012channelFileNameBase=[ basedir fnamedir fnamebase '_ch' num2str(channelLoop)  ];
    1113imgName=[ plots_dir fnamebase '_ch' num2str(channelLoop)  ];
    1214imgTitle=[ strrep(fnamebase,'_',' ') ' - ch' num2str(channelLoop) ];
     15
    1316if ((exist([ channelFileNameBase  '.csv.gz' ],'file'))||(exist([ channelFileNameBase  '.csv' ],'file')))
    1417    if (exist([ channelFileNameBase  '.csv.gz' ],'file'))
    15       disp(['Reading ' channelFileNameBase  '.csv.gz' ]);
    16       cmdGzip = ['gzip -c -d ' channelFileNameBase  '.csv.gz' ];
    17       [~, channelContent] = system(cmdGzip);
     18        disp(['Reading ' channelFileNameBase  '.csv.gz' ]);     
     19        gunzip([channelFileNameBase '.csv.gz']);
     20        data_num=csvread([channelFileNameBase '.csv'],2);
     21        delete([channelFileNameBase '.csv'])
    1822    elseif exist([ channelFileNameBase  '.csv' ],'file')
    19       disp(['Reading ' channelFileNameBase  '.csv' ]);
    20         fid=fopen([ channelFileNameBase  '.csv' ]);
    21         channelContent=fread(fid);
    22         fclose(fid);
     23        disp(['Reading ' channelFileNameBase  '.csv' ]);
     24        data_num=csvread([ channelFileNameBase  '.csv' ],2);
    2325    end %read file
    24     beginTxt = strfind(char(channelContent), '[V]')+5;
    25     data = char(channelContent(beginTxt:end)');
    26     data2 = strrep(data', ',', ' ');
    27     data_num=str2num(data2);
    28     if (length(channelContent) < 5000)
     26   
     27if (length(data_num) < 100)
    2928        warning_sp(['Channel ' num2str(channelLoop) ' has a short/empty content in ' channelFileNameBase ])
    3029    return
    3130    end % empty channel
    32     elseif exist([ basedir fnamedir fnamebase  '_image.png' ],'file')
    33           fname=[ basedir fnamedir fnamebase  '_image.png' ]
    34           disp(['Reading ' fname  ' - ch ' num2str(channelLoop) ]);
    35 fileinfo=dir(fname);
    36 if(fileinfo.bytes==0) 
    37   warning_sp( [ 'File ' fname ' is empty!']);
    38 data_num=[];
    39 else
    40       try
    41       imgdata=imread([ basedir fnamedir fnamebase  '_image.png' ]);
    42       catch
    43       [status,res]=system( [ ' ls -lart '  basedir fnamedir  fnamebase  '_image.png ' ] )
    44       warning_sp([ 'Unable to read '  basedir fnamedir  fnamebase  '_image.png ; ls result: ' res  ]);
    45       imgdata=[]
    46       end
     31elseif exist([ basedir fnamedir fnamebase  '_image.png' ],'file')
     32   
     33        fname=[ basedir fnamedir fnamebase  '_image.png' ]
     34        disp(['Reading ' fname  ' - ch ' num2str(channelLoop) ]);
     35    fileinfo=dir(fname);
     36   
     37    if(fileinfo.bytes==0) 
     38        warning_sp( [ 'File ' fname ' is empty!']);
     39        data_num=[];
     40    else
     41        try
     42        imgdata=imread([ basedir fnamedir fnamebase  '_image.png' ]);
     43        catch
     44        [status,res]=system( [ ' ls -lart '  basedir fnamedir  fnamebase  '_image.png ' ] )
     45        warning_sp([ 'Unable to read '  basedir fnamedir  fnamebase  '_image.png ; ls result: ' res  ]);
     46        imgdata=[]
     47        end
     48   
    4749        if (~(isempty(imgdata)))
    48       chdata=imgdata(:,:,2)-imgdata(:,:,1);
    49       nHorizontalPoints=950;
    50       croppedchdata=chdata(52:652,41:41-1+nHorizontalPoints);
    51       [val,pos]=max(croppedchdata);
    52       idxnonzero=find(pos>10);
    53       data_num=zeros(nHorizontalPoints,2);
    54       data_num_tmp(:,1)=(1:nHorizontalPoints)*channelParameters.scaleX;
    55       data_num_tmp(:,2)=(255-pos)*channelParameters.scaleY;
    56       data_num=data_num_tmp(idxnonzero,1:2);
    57 showScopeImg=0;
    58 if (showScopeImg==1)
    59       figure(201)
    60       clf
    61       image(imgdata)
    62       figure(202)
    63       clf
    64       image(chdata)
    65       figure(203)
    66       clf
    67       image(croppedchdata)
    68       figure(205)
    69       clf
    70       plot(data_num(:,1),data_num(:,2))
    71 end % showScopeImg
    72       else
    73         data_num=[]
    74       end % imgdata not empty
     50        chdata=imgdata(:,:,2)-imgdata(:,:,1);
     51        nHorizontalPoints=950;
     52        croppedchdata=chdata(52:652,41:41-1+nHorizontalPoints);
     53        [val,pos]=max(croppedchdata);
     54        idxnonzero=find(pos>10);
     55        data_num=zeros(nHorizontalPoints,2);
     56        data_num_tmp(:,1)=(1:nHorizontalPoints)*channelParameters.scaleX;
     57        data_num_tmp(:,2)=(255-pos)*channelParameters.scaleY;
     58        data_num=data_num_tmp(idxnonzero,1:2);
     59        showScopeImg=0;
     60        if (showScopeImg==1)
     61            figure(201)
     62            clf
     63            image(imgdata)
     64            figure(202)
     65            clf
     66            image(chdata)
     67            figure(203)
     68            clf
     69            image(croppedchdata)
     70            figure(205)
     71            clf
     72            plot(data_num(:,1),data_num(:,2))
     73        end % showScopeImg
     74        else
     75            data_num=[];
     76    end % imgdata not empty
    7577end %file is empty
    7678else
    … …  
    8991   
    9092end
    91 
     93       %% check if signal is inverted
     94       
     95    if  ( isfield(channelParameters, 'Inverted'))
     96    if (channelParameters.Inverted==1)
     97        data_num(:,2)=-data_num(:,2);% invert back
     98    end
     99   
     100    end
     101   
    92102for iPulseMode=1:2
    93103    if (iPulseMode==1)
    … …  
    102112        error_sp('Wrong pulse mode');
    103113    end
    104    
     114    %% FFT filtering
     115        if  ( isfield(channelParameters, 'Filtering'))
     116     FFT=fft(data_num(:,2));
     117 Am=abs(FFT);
     118 Ph=angle(FFT);
     119 d=channelParameters.Filtering;%filter limit
     120 Am((d+1):(end-d))=0;
     121 S=Am.*exp(1i.*Ph);
     122 SS=real(ifft(S));
     123 data_num(:,2)=SS(:);   
     124        end
     125   
     126    %%
    105127   
    106128    sigDuration=sigEnd-sigStart;
    … …  
    125147    %%% signal and background calculation
    126148   
    127     sigstartidx=find(data_num(:,1)>sigStart, 1 );
    128     sigstopidx=find(data_num(:,1)>sigEnd, 1 );
    129    
    130     bgBfStartIdx=find(data_num(:,1)>bgBfStart, 1 );
    131     bgAfStartIdx=find(data_num(:,1)>bgAfStart, 1 );
    132    
    133     bgBfEndIdx=find(data_num(:,1)>bgBfEnd, 1 );
    134     bgAfEndIdx=find(data_num(:,1)>bgAfEnd, 1 );
     149    sigstartidx=find(data_num(:,1)>=sigStart, 1 );
     150    sigstopidx=find(data_num(:,1)>=sigEnd, 1 );
     151   
     152    bgBfStartIdx=find(data_num(:,1)>=bgBfStart, 1 );
     153    bgAfStartIdx=find(data_num(:,1)>=bgAfStart, 1 );
     154   
     155    bgBfEndIdx=find(data_num(:,1)>=bgBfEnd, 1 );
     156    bgAfEndIdx=find(data_num(:,1)>=bgAfEnd, 1 );
    135157   
    136158    %%% signal
    137     signalMean=sum(data_num(sigstartidx:sigstopidx,2))/(sigstopidx-sigstartidx+1);
    138     signalRMS=rms(data_num(sigstartidx:sigstopidx,2)-signalMean);
     159 
     160       
     161     SigInWindow=  data_num(sigstartidx:sigstopidx,2);
     162    [~,pos]=max(SigInWindow);%find position of maximum
     163    diap=0;
     164    if ((pos- diap)<1)
     165        lowBound=1;
     166    else
     167        lowBound=pos- diap;
     168    end
     169    if ((pos+ diap)>length(SigInWindow))
     170        highBound=length(SigInWindow);
     171    else
     172                highBound=pos+ diap;
     173    end
     174    signalMean=mean(SigInWindow(lowBound:highBound));
     175    signalRMS=rms(SigInWindow(lowBound:highBound)-signalMean);
     176
     177    % signalMean=mean(data_num(sigstartidx:sigstopidx,2));
     178   % signalRMS=rms(data_num(sigstartidx:sigstopidx,2)-signalMean);
    139179    signalNegPeakRaw=min(data_num(sigstartidx:sigstopidx,2));
    140180   
    141     bgBfMean=sum(data_num(bgBfStartIdx:bgBfEndIdx,2))/(bgBfStartIdx-bgBfEndIdx+1);
     181    bgBfMean=mean(data_num(bgBfStartIdx:bgBfEndIdx,2));
    142182    bgBfRMS=rms(data_num(bgBfStartIdx:bgBfEndIdx,2)-bgBfMean);
    143     bgAfMean=sum(data_num(bgAfStartIdx:bgAfEndIdx,2))/(bgAfStartIdx-bgAfEndIdx+1);
     183    bgAfMean=mean(data_num(bgAfStartIdx:bgAfEndIdx,2));
    144184    bgAfRMS=rms(data_num(bgAfStartIdx:bgAfEndIdx,2)-bgAfMean);
    145185   

  • SPESO/ana2015/get_parameters.m

    r634 r643  
    55if (dateNum>=20180101)
    66  error_sp([ 'NO parameters defined on ' dateStr ]);
     7  elseif (dateNum>=20170311)
     8  myParameters.ch(1).on=1;
     9 
     10  myParameters.ch(1).SPM.sigStart=-7e-8;
     11  myParameters.ch(1).SPM.sigEnd=-5.5e-8;
     12  myParameters.ch(1).SPM.riseTime=1.0e-8;
     13
     14  myParameters.ch(1).LPM.sigStart=-15e-8;
     15  myParameters.ch(1).LPM.sigEnd=15e-8;
     16  myParameters.ch(1).LPM.riseTime=1e-8;
     17  myParameters.ch(1).Filtering=50;
     18  myParameters.ch(1).Inverted=0; 
     19  myParameters.ch(2).on=1;
     20 
     21  myParameters.ch(2).SPM.sigStart=-8e-8;
     22  myParameters.ch(2).SPM.sigEnd=-6.5e-8;
     23  myParameters.ch(2).SPM.riseTime=1.0e-8;
     24
     25  myParameters.ch(2).LPM.sigStart=-15e-8;
     26  myParameters.ch(2).LPM.sigEnd=15e-8;
     27  myParameters.ch(2).LPM.riseTime=1e-8;
     28   myParameters.ch(2).Inverted=0; 
     29    myParameters.ch(2).Filtering=50;
     30
     31  myParameters.ch(3).on=1;
     32 
     33  myParameters.ch(3).SPM.sigStart=-0.73e-7;
     34  myParameters.ch(3).SPM.sigEnd=-0.61e-7;
     35  myParameters.ch(3).SPM.riseTime=2.0e-8;
     36
     37  myParameters.ch(3).LPM.sigStart=-15e-8;
     38  myParameters.ch(3).LPM.sigEnd=15e-8;
     39  myParameters.ch(3).LPM.riseTime=1e-8;
     40   myParameters.ch(3).Inverted=1; 
     41    myParameters.ch(3).Filtering=200;
     42
     43  myParameters.ch(4).on=0;
     44elseif (dateNum>=20161208)
     45  myParameters.ch(1).on=1;
     46 
     47  myParameters.ch(1).SPM.sigStart=2.65e-8;
     48  myParameters.ch(1).SPM.sigEnd=3.05e-8;
     49  myParameters.ch(1).SPM.riseTime=1.0e-8;
     50
     51  myParameters.ch(1).LPM.sigStart=-1e-8;
     52  myParameters.ch(1).LPM.sigEnd=1e-8;
     53  myParameters.ch(1).LPM.riseTime=1e-8;
     54  myParameters.ch(1).Filtering=50;
     55  myParameters.ch(1).Inverted=0; 
     56  myParameters.ch(2).on=1;
     57 
     58  myParameters.ch(2).SPM.sigStart=1.90e-8;
     59  myParameters.ch(2).SPM.sigEnd=2.20e-8;
     60  myParameters.ch(2).SPM.riseTime=1.0e-8;
     61
     62  myParameters.ch(2).LPM.sigStart=-1e-8;
     63  myParameters.ch(2).LPM.sigEnd=1e-8;
     64  myParameters.ch(2).LPM.riseTime=1e-8;
     65   myParameters.ch(2).Inverted=0; 
     66    myParameters.ch(2).Filtering=50;
     67
     68  myParameters.ch(3).on=1;
     69 
     70  myParameters.ch(3).SPM.sigStart=-0.73e-7;
     71  myParameters.ch(3).SPM.sigEnd=-0.61e-7;
     72  myParameters.ch(3).SPM.riseTime=2.0e-8;
     73
     74  myParameters.ch(3).LPM.sigStart=-1e-8;
     75  myParameters.ch(3).LPM.sigEnd=1e-8;
     76  myParameters.ch(3).LPM.riseTime=1e-8;
     77       myParameters.ch(3).Inverted=1; 
     78    myParameters.ch(3).Filtering=200;
     79
     80  myParameters.ch(4).on=0;
    781elseif (dateNum>=20161122)
    882  myParameters.ch(1).on=1;
Note: See TracChangeset for help on using the changeset viewer.