Acquisition procedure

function master(experimenttime, recordingtime, duration, rate, savedrate, plotresolution,savepath)
    variable experimenttime, recordingtime, duration, rate, savedrate, plotresolution
    string savepath
    variable psdflag=1 // first time, make the periodogram
    NewPath/O Path1,savepath//defines where to save the file
       
        // rate= sampling rate for data acquiring, like 40000 Hz
    // savedrate= rate for high time-resolution data recording, like 200 Hz
    // plotresolution= interval between each point on strip graph
    // recordingtime= interval for saving on disk, like 10s or 30s
    // experimenttime= global experiment time, after which acquisition stops, like 3600 s (break stops before and saves)

<!--break-->
    string timenow=date()+" "+time() // date and time now
    string filename
    variable t0 = stopMSTimer(-2)/10^6 // consults the microseconds timer, from the moment of PC startup
    variable t1=0, dt1=0, dt2=0 // time init
    if (plotresolution<duration)
        plotresolution=duration // cannot update faster than acquire
    endif
    make/O/N=0 timestrip //store the time
    make/O/N=0 deflstrip  //where to store the deflection every plotresolution>duration
    make /O/N=0 variancestrip //where to store the variance every plotresolution>duration
    setscale/P x, t1, (plotresolution), "s", timestrip, deflstrip, variancestrip
    setscale/P y, 0,1,"V" deflstrip, variancestrip
    setscale/P y,0,1,"s",timestrip
    //display
    display/W=(1,1,800,150) deflstrip vs timestrip
    AppendToGraph/L=L2 variancestrip vs timestrip
    ModifyGraph axisEnab(left)={0,0.45}
    ModifyGraph axisEnab(L2)={0.55,1}
    ModifyGraph freePos(L2)=0
    ModifyGraph lblPosMode(L2)=2
    ModifyGraph rgb(deflstrip)=(0,0,0)
 
 
    make /O/N=0 storeaverage, storevariance, storetime, waveaverage, wavevariance, wavetime
    display/W=(8, 310, 400,550) storeaverage vs storetime
    variable acq_cycle=0 // counter for acquisition cycle
    do //big loop for acquisition and plot
       t1=stopMSTimer(-2)/10^6-t0 // evaluate the experimental time
       //print "t1=",t1
       waveacquire(t1, duration, rate, savedrate, psdflag) // DATA ACQUISITION
        dt1+=stopMSTimer(-2)/10^6-t0-t1 // increment the delta time for plotting
       dt2+=stopMSTimer(-2)/10^6-t0-t1 // increment the delta time for saving
 
       insertpoints numpnts(storeaverage), numpnts(waveaverage), storeaverage // extends storeaverage
       storeaverage[numpnts(storeaverage)-numpnts(waveaverage),
           &&linewrap numpnts(storeaverage)-1]=waveaverage[p-numpnts(storeaverage)+numpnts(waveaverage)-1]
       insertpoints numpnts(storevariance), numpnts(wavevariance), storevariance // extends storevariance
       storevariance[numpnts(storevariance)-numpnts(wavevariance),
           &&linewrap numpnts(storevariance)-1]=wavevariance[p-numpnts(storeaverage)+numpnts(waveaverage)-1]
       insertpoints numpnts(storetime), numpnts(wavetime), storetime // extends storetime
       storetime[numpnts(storetime)-numpnts(wavetime),
           &&linewrap numpnts(storetime)-1]=wavetime[p-numpnts(storeaverage)+numpnts(waveaverage)-1]
 
       acq_cycle=+1 //increment the counter for acquisition cycle
       if (psdflag==1)
         psdflag=0
       endif // switches off psd
       //print dt1," ",dt2
       if (dt1>=plotresolution) // after exiting this loops, it adds a plot point ...
         // insertpoints numpnts(deflstrip), 1, deflstrip, variancestrip // makes room
         timestrip[numpnts(timestrip)]={stopMSTimer(-2)/10^6-t0} // to store the time too
         deflstrip[numpnts(deflstrip)]={mean(storeaverage, pnt2x(storeaverage,
           &&linewrap numpnts(storeaverage)-acq_cycle*rate/savedrate), pnt2x(storeaverage, numpnts(storeaverage)))}
         variancestrip[numpnts(variancestrip)]={mean(storevariance, pnt2x(storevariance,
           &&linewrap numpnts(storevariance)-acq_cycle*rate/savedrate), pnt2x(storevariance,
           &&linewrap numpnts(storevariance)))/sqrt(acq_cycle*rate/savedrate)}
         doupdate // explicitly ask to update the plots (in hope it has not done it before at insertpoint)
        // print "plotplot"
         dt1=0 //resets delta time for plotting
         acq_cycle=0 // resets acquisition cycle, so it later knows what to use for plotting
       endif
       if  (dt2>=recordingtime)
         // saves storeaverage and storevariance on file with automatic names
         // empty storeaverage and storevariance for next acquisitions
         timenow=date()+" "+time() // updates time
         filename="DeflectionStrip "+num2str(savedrate)+"Hz on "+timenow+"-"+num2istr(dt2)+".txt"
         save /G/P=Path1 storetime storeaverage storevariance as filename  
         redimension/N=0 storetime, storeaverage, storevariance //clean up
         dt2=0 //resets delta time for recording
         psdflag=1
         endif
        while (t1<experimenttime)
    // now save storeaverage and storevariance,
    timenow=date()+" "+time() // updates time
    filename="DeflectionStrip "+num2str(savedrate)+"Hz on "+timenow+"-"+num2istr(dt2)+".txt" // for saving
    save /G/P=Path1 storetime storeaverage storevariance as filename       
    // and also save deflstrip and variancestrip
    timenow=date()+" "+time() // updates time
    filename="ExperimentStrip "+"finished on "+timenow+".txt" // for saving
    save /G/P=Path1 timestrip deflstrip variancestrip as filename  
    //print "I want to save deflstrip and variancestrip"
    // now acquire last PSD
    waveacquire(t1, duration, rate, savedrate, 1) // DATA ACQUISITION (and saving of PSD)
end // of function

function waveacquire(t1, duration, rate, savedrate, psdflag)
    variable t1, duration, rate, savedrate, psdflag
    string timenow, filename
    variable i=0
    variable blocksize=rate/savedrate // blocksize for block averaging
    variable segmentsize=0
    if (duration*rate>8192*4)
            segmentsize=8192// segment size for Periodogram calculation
    else
            segmentsize=(duration*rate)/5
    endif
    //print segmentsize
    make/O/N=(rate*duration) wave0
    make/O/N=(rate*duration/blocksize) wavevariance // wave variance decimated
    make/O/N=(rate*duration/blocksize) waveaverage // wave after decimation
    make/O/N=(rate*duration/blocksize) wavetime // to store the time points
    make/O/N=0 media //to store the recordings
    make/O/N=0 timepoints //to store the time of recordings
    make/O/N=(segmentsize/2+1) periodogramdomain

    // variable timezero=(stopMSTimer(-2)-t0)/10^6 // resets the time
    //print timezero
    setscale/P x, t1, (1/rate), "s", wave0 // sets the spacing of the wave, so the sampling rate
    // next is commented so this is a SIMULATION for acquisition
    DAQmx_scan/DEV="dev4" WAVES="wave0, 2;" //acquires from DEV device on channel 2 and puts in wave0
    // must create a rate*duration size wave0
    // for adding a frequency, use sin(freq*x*2*Pi) x is already scaled in 1/rate seconds
    // variable phase=Pi*gnoise(1) // to start the noise at different phases
    // wave0=gnoise(1)+10*sin(3000*x*2*Pi)+sin(0.1*x*2*Pi+phase)// 3000 Hz vibration+ 0.1Hz vibration + gaussian noise
    // sleep/S duration // pause to simulate acquisition and do not mess the timings
    if (psdflag==1)
            DSPPeriodogram/Q/DLSG/NOR=(segmentsize*segmentsize/2)/SEGN={(segmentsize),(segmentsize/2)}/Win=hanning wave0
            timenow=date()+" "+time() // updates time
            filename="PSD "+num2str(savedrate)+"Hz on "+timenow+".txt" // for saving
            periodogramdomain=pnt2x(W_periodogram, p)
            save /G/P=Path1 periodogramdomain W_periodogram as filename
    endif

    wavevariance=variance(wave0, pnt2x(wave0, blocksize*p), pnt2x(wave0, (blocksize*(p+1))-1))/blocksize
    setscale/P x, t1, blocksize*(1/rate), "s", wavevariance // sets the spacing of the wave, so the sampling rate
    waveaverage=sum(wave0, pnt2x(wave0, blocksize*p), pnt2x(wave0, (blocksize*(p+1))-1))/blocksize
    setscale/P x, t1, blocksize*(1/rate), "s", waveaverage // sets the spacing of the wave, so the sampling rate
    wavetime=pnt2x(wavevariance, p) // timings from the scale (?)
    variable stoptimer=stopMSTimer(-2) // stops the microseconds timer
end
<!--break-->