환무의 블로그

data를 자석 높이별로 구분한 뒤, 분석하여 bead에 가해진 힘 결정

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%  Top level section 3:  FX3_ANALYZE %%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if 1

    %%% This is the third section in a series of routines to analyze

    %%% force-extension data taken with the new multi bead code for MT

    %%% (J. P. Cnossen, D. Dulin and N. H. Dekker Rev Sci Instrum 85, (2014).)

    %%%

    %%% This section assumes that you have read the (z-offset corrected) raw

    %%% data traces

    %%%

    %%% Author: Jan Lipfert

    %%% Date: 2013-09-16, updated 2014-07-17

    clc;

    plotflag = 1; %%% Whether or not to show the fits

    %%% Some variables ----------------------------------------------

    kT = 4.1;

    Rbead = 1.4;      % Bead radius in mum

    %% 1um bead 사용시 수정해야 함

    %%% ---------------------------------------------------------------

    %%% Figure out where the magnets are moving, from the motor file

    %%% This determines the "plateaus", where the nmagnet height is

    %%% constant and where we we want to analyze the forces

    %%% ---------------------------------------------------------------

    Nsmooth_zmag = 200;

    Nsmooth_dzmag = 200;

    small = 10^(-5); %%% Threshold to determine where it is moving

    %%% Some variables for error checking

    Nmin_points_plat = 100; % Minimum number of points in a plateau

    zmags = [];

    if 1

        zmag_smooth = smooth(zmag, Nsmooth_zmag, 'moving');

        diff_zmag = smooth(diff(zmag_smooth), Nsmooth_dzmag, 'lowess');

        %%% Find all the "plateaus", i.e. the sets of points where the magnets

        %%% are not moving up or down

        Nfirst =1; Nlast  =1;

        %%% Check whether we are starting in a plateau

        if (abs(diff_zmag(1)) < small & abs(diff_zmag(2)) < small)

            tplat(1).first = 1;

            Nfirst = Nfirst +1;

        end

        for i=2:length(diff_zmag)

            if (abs(diff_zmag(i)) < small & abs(diff_zmag(i-1)) > small)

                tplat(Nfirst).first = i;

                Nfirst = Nfirst +1;

            end

            if (abs(diff_zmag(i)) > small & abs(diff_zmag(i-1)) < small)

                tplat(Nlast).last = i;

                Nlast = Nlast +1;

            end

        end

        %%% Check whether we are ending in a plateau

        if (abs(diff_zmag(end)) < small & abs(diff_zmag(end-1)) < small)

            tplat(Nlast).last = length(diff_zmag);

            Nlast = Nlast +1;

        end

        Nplat = Nfirst-1;

        display(['Found ' num2str(Nplat) ' Zmag plateaus']);

        %display(['Nfirst ' num2str(Nfirst) ' Nlast ' num2str(Nlast)]);

        %%% Throw out plateaus that have fewer than Nmin_points_plat points

        %%% This is not quite it yet, since it keeps the plateaus separate that

        %%% are interupted by a "fake" (i.e. too short) plateau, but getting it

        %%% perfect is tricky

        if 1

            count = 1;

            Ngoodplat = 0;

            for j = 1:Nplat

                if length(tplat(j).first:tplat(j).last) > Nmin_points_plat

                    plat(count).first = tplat(j).first;

                    plat(count).last = tplat(j).last;

                    Ngoodplat = Ngoodplat + 1;

                    count = count + 1;

                else

                    display(['Plateau # ' num2str(j) ' has only ' num2str(length(Nmin_points_plat)) ' points!' ])

                end

            end

            Nplat = Ngoodplat;

        end

        %%% Plot the magnet height information, after smoothing, with

        %%% plateaus annotated

        if 1

            figure(1);clf; hold on; box on; %%% Magnet rot. vs. time

            platind = find(abs(diff_zmag) < small);

            faketime = 1:length(zmag_smooth);

            plot(faketime , zmag_smooth, 'b-')

            plot(faketime(platind), zmag_smooth(platind), 'r.')

            for j = 1:Nplat

                plot(faketime(plat(j).first) , zmag_smooth(plat(j).first), 'ko', 'markersize', 10)

                plot(faketime(plat(j).last) , zmag_smooth(plat(j).last), 'mo', 'markersize', 10)

                %%% Capture the output

                zmags = [zmags zmag_smooth(plat(j).last)];

            end

            set(gca, 'fontsize', 14, 'fontweight', 'bold', 'linewidth', 1,'TickLength',[0.02 0.02])

            xlabel('Time'); ylabel('Zmag')

            title(['Red = plateaus; Black / magneta circles = start / stop; Found ' num2str(Nplat) ' plateaus.' ])

            figure(2);clf; hold on; box on; %%% Derivative of magnet rot. vs. time

            plot(1:length(diff_zmag) , diff_zmag, 'b-')

            set(gca, 'fontsize', 14, 'fontweight', 'bold', 'linewidth', 1,'TickLength',[0.02 0.02])

            xlabel('Time'); ylabel('dzmag/dt')

        end

    end

    %%% ------- END PLATEAU FINDIND --------- %%%

     %pause;

    %%% ------------------------------------------------------- %%%

    %%% Look over the beads and determine the force for each plateau

    %%% ------------------------------------------------------- %%%

    if 1

        for i=1:Nbeads

            bead(i).ext = zeros(1,Nplat);

            bead(i).Fx_real = zeros(1,Nplat);

            bead(i).Fy_real = zeros(1,Nplat);

            bead(i).PSDforce = zeros(1,Nplat);

            bead(i).AVforce = zeros(1,Nplat);

            bead(i).fcorner = zeros(1,Nplat);

            bead(i).SAforce = zeros(1,Nplat);

            bead(i).fcorner2 = zeros(1,Nplat);

            %%% Use the script "analyze_one_trace" to determine the force for each trace

            for k=1:Nplat

                if length(plat(k).first:plat(k).last) > Nmin_points_plat

                    [bead(i).ext(k),...

                        bead(i).Fx_real(k),...

                        bead(i).Fy_real(k),...

                        bead(i).PSDforce(k),...

                        PSDfit,...

                        bead(i).AVforce(k),...

                        AVfit,...

                        bead(i).fcorner(k),...

                        SAfit,...

                        bead(i).SAforce(k),...

                        bead(i).fcorner2(k)] = ...

                        analyze_one_trace(bead(i).time(plat(k).first:plat(k).last),...

                        bead(i).x(plat(k).first:plat(k).last),...

                        bead(i).y(plat(k).first:plat(k).last),...

                        bead(i).z(plat(k).first:plat(k).last),...

                        i, F_IND, freq, plotflag); %For this example, we replace x dimension with y dimension to keep the size of data small

                    display(['Working on bead # ' num2str(i) ', plateau number ' num2str(k)])

                    if plotflag == 1

                                   %          pause;

                    end

                end

            end

            %% 각 높이에서 힘 및 코너 프리컨시 추출

            %%% Fit a two exponential model to the data

            %%% The model: F0 + Amp * exp(Zmag/l_decay) + Amp * exp(Zmag/l_decay)

            %%% Note that this is just an empirical fit, no real physics

            s=exp2fit(zmags, bead(i).PSDforce,2);

            bead(i).PSD_F0 = s(1);

            bead(i).PSD_Amp = s(2);

            bead(i).PSD_2_decay = s(3);

            bead(i).PSD_Amp_2 = s(4);

            bead(i).PSD_2_decay_2 = s(5);

            bead(i).PSD_Zfit = min(zmags):0.01:max(zmags);

            bead(i).PSD_Ffit = s(1)+s(2)*exp(-bead(i).PSD_Zfit/s(3))+s(4)*exp(-bead(i).PSD_Zfit/s(5));

            s=exp2fit(zmags, bead(i).AVforce,2);

            bead(i).AV_F0 = s(1);

            bead(i).AV_Amp = s(2);

            bead(i).AV_2_decay = s(3);

            bead(i).AV_Amp_2 = s(4);

            bead(i).AV_2_decay_2 = s(5);

            bead(i).AV_Zfit = min(zmags):0.01:max(zmags);

            bead(i).AV_Ffit = s(1)+s(2)*exp(-bead(i).AV_Zfit/s(3))+s(4)*exp(-bead(i).AV_Zfit/s(5));

            s=exp2fit(zmags, bead(i).SAforce,2);

            bead(i).SA_F0 = s(1);

            bead(i).SA_Amp = s(2);

            bead(i).SA_2_decay = s(3);

            bead(i).SA_Amp_2 = s(4);

            bead(i).SA_2_decay_2 = s(5);

            bead(i).SA_Zfit = min(zmags):0.01:max(zmags);

            bead(i).SA_Ffit = s(1)+s(2)*exp(-bead(i).SA_Zfit/s(3))+s(4)*exp(-bead(i).SA_Zfit/s(5));

            if 1

                %%% Fit the WLC chain model to the data

                %%% This is real physics and determines Lp (persistence length)

                %%% and Lc (the contour length)

                [bead(i).PSD_Lp, bead(i).PSD_Lc, bead(i).PSD_Forcefit, bead(i).PSD_Extfit] = inExtWLCfit(bead(i).ext, bead(i).PSDforce, 0.1);

                [bead(i).AV_Lp, bead(i).AV_Lc, bead(i).AV_Forcefit, bead(i).AV_Extfit] = inExtWLCfit(bead(i).ext, bead(i).AVforce, 0.1);

                [bead(i).SA_Lp, bead(i).SA_Lc, bead(i).SA_Forcefit, bead(i).SA_Extfit] = inExtWLCfit(bead(i).ext, bead(i).SAforce, 0.1);

            end

            %%% Show a Force vs. Zmag plot for this bead

            if 1

                %%% Plot F vs. Zmag, linear

                figure(3); clf; hold on; box on;

                plot(zmags, bead(i).Fx_real, 'ko', 'markersize', 15, 'linewidth', 1)

                plot(zmags, bead(i).PSDforce, 'bo', 'markersize', 15, 'linewidth', 2)

                plot(zmags, bead(i).AVforce, 'ro', 'markersize', 15, 'linewidth', 2)

                plot(zmags, bead(i).SAforce, 'go', 'markersize', 15, 'linewidth', 2)

                %%% Plot the fits

                plot(bead(i).PSD_Zfit,  bead(i).PSD_Ffit, 'b-', 'linewidth', 2)

                plot(bead(i).AV_Zfit,  bead(i).AV_Ffit, 'r-', 'linewidth', 2)

                plot(bead(i).SA_Zfit,  bead(i).SA_Ffit, 'g-', 'linewidth', 2)

                set(gca, 'fontsize', 18, 'fontweight', 'bold', 'linewidth', 1,'TickLength',[0.02 0.02])

                xlabel('Magnet position (mm)');

                ylabel('Force (pN)');

                legend('Real space', 'PSD', 'Allan', 'SA')

                %%% Show the fitted parameter from the Allan variance fit in the title

                title(['From AV fit: Amp=' num2str(bead(i).AV_Amp,3) ';  l(decay)=' num2str(bead(i).AV_2_decay,3) '; Amp2=' num2str(bead(i).AV_Amp_2,3) ';  l(decay)2=' num2str(bead(i).AV_2_decay_2,3) ';   F0=' num2str(bead(i).AV_F0,2) ])

                %%% Plot F vs. Zmag, log scaling

                figure(4); clf; hold on; box on;

                plot(zmags, bead(i).Fx_real, 'ko', 'markersize', 15, 'linewidth', 1)

                plot(zmags, bead(i).PSDforce, 'bo', 'markersize', 15, 'linewidth', 2)

                plot(zmags, bead(i).AVforce, 'ro', 'markersize', 15, 'linewidth', 2)

                plot(zmags, bead(i).SAforce, 'go', 'markersize', 15, 'linewidth', 2)

                %%% Plot the fits

                plot(bead(i).PSD_Zfit,  bead(i).PSD_Ffit, 'b-', 'linewidth', 2)

                plot(bead(i).AV_Zfit,  bead(i).AV_Ffit, 'r-', 'linewidth', 2)

                plot(bead(i).SA_Zfit,  bead(i).SA_Ffit, 'g-', 'linewidth', 2)

                set(gca, 'fontsize', 18, 'fontweight', 'bold', 'linewidth', 1,'TickLength',[0.02 0.02])

                xlabel('Magnet position (mm)');

                ylabel('Force (pN)');

                set(gca, 'yscale', 'log')

                legend('Real space', 'PSD', 'Allan', 'SA')

                title(['Bead # ' num2str(i)])

            end

            %%% Show a FX plot for this bead

            if 1

                %%% Plot force vs. DNA extension

                figure(5); clf; hold on; box on;

                plot(bead(i).ext, bead(i).Fx_real, 'ko', 'markersize', 15, 'linewidth', 1)

                plot(bead(i).ext, bead(i).PSDforce, 'bo', 'markersize', 15, 'linewidth', 2)

                plot(bead(i).ext, bead(i).AVforce, 'ro', 'markersize', 15, 'linewidth', 2)

                plot(bead(i).ext, bead(i).SAforce, 'go', 'markersize', 15, 'linewidth', 2)

                if 1

                    plot(bead(i).PSD_Extfit, bead(i).PSD_Forcefit, 'b-', 'markersize', 15, 'linewidth', 2)

                    plot(bead(i).AV_Extfit, bead(i).AV_Forcefit, 'r-', 'markersize', 15, 'linewidth', 2)

                    plot(bead(i).SA_Extfit, bead(i).SA_Forcefit, 'g-', 'markersize', 15, 'linewidth', 2)

                end

                set(gca, 'fontsize', 18, 'fontweight', 'bold', 'linewidth', 1,'TickLength',[0.02 0.02])

                xlabel('Extension (um)');

                ylabel('Force (pN)');

                legend('Real space', 'PSD', 'Allan', 'SA', 'location', 'northwest')

                title(['From AV fit: Lp=' num2str(bead(i).AV_Lp*1000,3) 'nm;  Lc=' num2str(bead(i).AV_Lc,3) 'um'  ])

                %%% Log plot

                figure(6); clf; hold on; box on;

                plot(bead(i).ext, bead(i).Fx_real, 'ko', 'markersize', 15, 'linewidth', 1)

                plot(bead(i).ext, bead(i).PSDforce, 'bo', 'markersize', 15, 'linewidth', 2)

                plot(bead(i).ext, bead(i).AVforce, 'ro', 'markersize', 15, 'linewidth', 2)

                plot(bead(i).ext, bead(i).SAforce, 'go', 'markersize', 15, 'linewidth', 2)

                if 1

                    plot(bead(i).PSD_Extfit, bead(i).PSD_Forcefit, 'b-', 'markersize', 15, 'linewidth', 2)

                    plot(bead(i).AV_Extfit, bead(i).AV_Forcefit, 'r-', 'markersize', 15, 'linewidth', 2)

                    plot(bead(i).SA_Extfit, bead(i).SA_Forcefit, 'g-', 'markersize', 15, 'linewidth', 2)

                end

                set(gca, 'fontsize', 18, 'fontweight', 'bold', 'linewidth', 1,'TickLength',[0.02 0.02])

                xlabel('Extension (um)');

                ylabel('Force (pN)');

                set(gca, 'yscale', 'log')

                legend('Real space', 'PSD', 'Allan', 'SA', 'location', 'northwest')

                title(['Bead # ' num2str(i)])

            end

            %%% Show a corner freq vs. Zmag plot for this bead

            if 1

                figure(7); clf; hold on; box on;

                plot(zmags, bead(i).fcorner2, 'ro', 'markersize', 15, 'linewidth', 2)

                set(gca, 'fontsize', 14, 'fontweight', 'bold', 'linewidth', 1,'TickLength',[0.02 0.02])

                xlabel('Magnet position (mm)');

                ylabel('f(corner, Hz)');

            end

            %         pause;

        end

    end

   % 2016.10.07 ??????

    %% 분석된 값들 저장

   outlength = length(zmags); 

   out = zeros(6,outlength);

   out(1,:)=zmags;

   out(2,:)=bead.Fx_real;

   out(3,:)=bead.Fy_real;

   out(4,:)=bead.PSDforce;

   out(5,:)=bead.AVforce;

   out(6,:)=bead.SAforce;

   out = transpose(out);

   out1 = zeros(3,2);

   out1(1,1) = bead.PSD_Lp;

   out1(2,1) = bead.AV_Lp;

   out1(3,1) = bead.SA_Lp;

   out1(1,2) = bead.PSD_Lc;

   out1(2,2) = bead.AV_Lc;

   out1(3,2) = bead.SA_Lc;

   save('ForceValues.txt','out','out1','-ascii')

end

%% data 읽어서 자석의 높이별로 구분

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%  Top level section 2: FX2_LOAD_RAW  %%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if 1

    %%% This is the second section in a series of routines to analyze

    %%% force-extension data taken with the new multi bead code for MT

    %%% (J. P. Cnossen, D. Dulin and N. H. Dekker Rev Sci Instrum 85, (2014).)

    %%%

    %%% This section is used to load the data for the force-extension

    %%% measurements, namely the large matrix of position traces for different

    %%% magnet heights.

    %%%

    %%% The script assumes that you have determined the z-offsets previously and

    %%% loads them z-offsets from a file

    %%%

    %%% Author: Jan Lipfert

    %%% Date: 2013-09-16

    clear all; clc; close all;

    %%% Index for each data set

    DATA = 1;

    %%% Write here some descriptive text about your data

    %%% ---

    %%% Example of force extension data

    %%% 21 kbp DNA, M270 beads, 1 mm gap verticaly oriented magnets

    if DATA ==1

        traces_file = 'bead.txt';

        motors_file = 'bead_motors.txt';

        zoffsets_file = 'FX_offsets.txt';

        Nref = 1;

        F_IND = 1; %%% whether to use x or y (1=x, 2=y)

        freq = 60; %%% Acquisition frequency in Hz

    end

    %% 읽을 파일 및 자석의 정렬방향, 카메라 속도 설정 

   %% (우리 셋업은 x축 방향으로 자석 정렬, 60Hz로 이미지 취득)

    %%%--- Read in the previously determined z-offsets from file

    zoff_data = load(zoffsets_file);

    zoffsets = zoff_data(:,2);

    %% 1번 함수에서 만들어 둔 바닥면에 대한 정보 취득

    %%%--- Read in data ---

    data = load(traces_file);

    zmag = load(motors_file);

    Nbeads = 1;

    %% data 및 자석에 대한 정보 loading

    %%%--- Parse the bead data ---

    for i=1:Nbeads

        display(['Parsing data for bead ' num2str(i) ' of ' num2str(Nbeads)])

        bead(i).time = 1:length(data(:,1));

        bead(i).time = bead(i).time/freq;% ????????

        bead(i).x = data(:,1);

        bead(i).y = data(:,2);

        bead(i).z = data(:,3);

    end

    %% 시간 설정 및 data 할당

    %%% Subtract previously determined z-offsets

    if 1

        if length(zoffsets) == Nbeads

            for i=1:Nbeads

                %bead(i).z = bead(i).z - zoffsets(i);

                bead(i).z = zoffsets(i)- bead(i).z ; % max(z) -Z ; our setup

                %bead(i).z = (zoffsets(i)- bead(i).z)*0.878 ; % max(z) -Z ; our setup

            end

            display('Subtracted pre-determined z-offsets.')

        else

            display(['Number of beads in the big data file (' num2str(Nbeads) ') does not agree with the number of beads in the zoffset file (' num2str(length(zoffsets)) ')'])

        end

    end

    %% 읽어둔 바닥면 정보 이용하여 z축 길이 결정

end

offset data 읽어서 바닥면 정의

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%  Top level section 1: FX1_DETERMINE_ZOFFSET  %%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if 1

    %%% This is the first section in a series of routines to analyze

    %%% force-extension data taken with the new multi bead code for MT

    %%% (J. P. Cnossen, D. Dulin and N. H. Dekker Rev Sci Instrum 85, (2014).)

    %%%

    %%% This section is used to determine the z-offset for each bead

    %%% ---

    %%% To determine the z-offsets, record bead position traces under conditions

    %%% where the beads are essentially at the surface, either by going to very

    %%% low (essentially zero) stretching force or by going to a large number

    %%% of supercoils /turns under low forces.

    %%%

    %%% - Set the file names and path informaiton in the "DATA" block below.

    %%% The variable "DATA" controls which data sets are read and analyzed.

    %%%

    %%% - The code determines the offset by smoothing the position traces with a

    %%% window of Nsmooth points abd takes the minimum of the smoothed traces

    %%% as the surface position = z-offset

    %%%

    %%% - The results get saved into a file of which the name is set through the

    %%% "output_name" variable. Format is first colum = bead number, second

    %%% column = offset in mum

    %%%

    %%% Author: Jan Lipfert

    %%% Date: 2013-09-16

    clear all; clc; close all;

    %%% Index to select which data set to run

    DATA = 1;

    %%% Write here some descriptive text about your data

    %%% ---

    %%% Example data set on FX measurements; M270 beads, 21 kbp DNA

    if DATA ==1

        traces_file = 'offset.txt';

        motors_file = 'offset_motors.txt';

        output_name = 'FX_offsets.txt';

        Nref = 1;

    end

   %% 읽을 파일 및 출력할 파일에 대한 이름 정의    

    Nsmooth = 100;

    %%%%%%%%%%%%%%%%%%%%%%%

    %%%--- Read in data ---

    %%%%%%%%%%%%%%%%%%%%%%%

    data = load(traces_file);

    zmag = load(motors_file); % for this example, motor file contains just the data of zmag

   %% 바닥면에 대한 bead data와 자석의 높이에 대한 정보값 loading    

    %%%--- Parse the bead data ---

    % For this example, there is just one bead, which is already corrected by reference

    % bead. The dimensions of y and z are the only columns saved for this example.

    for i=1

        display(['Parsing data for bead ' num2str(i) ' of ' num2str(1)])

        bead(i).time = 1:length(data(:,1));

        bead(i).x = data(:,1);

        bead(i).y = data(:,2);

        bead(i).z = data(:,3);

    end

       %% bead 구조체에 loading 한 data 할당

    %%% --- Plot the motor information --- %%%

    if 1

        %%% Plot the magnet height information

        figure(100);clf; hold on; box on; %%% Magnet rot. vs. time

        faketime = 1:length(zmag);

        plot(faketime , zmag, 'r-')

        set(gca, 'fontsize', 16, 'linewidth', 1, 'fontweight', 'bold','TickLength',[0.02 0.02]);

        xlabel('Time (frames)'); ylabel('Zmag (mm)')

        title('Motor height')

    end

    z_offsets = [];

    %%% Loop over beads and show the height information

    if 1

        for i=1 % for this example, there is just one bead

            figure(1); clf; hold on; box on;

            plot(bead(i).time, bead(i).z, 'k-', 'linewidth', 1)

            %%% Smooth and find minimum

            smooth_z = smooth(bead(i).z, Nsmooth, 'moving');

            plot(bead(i).time, smooth_z, 'r-', 'linewidth', 2)

            %[min_z ind] = min(smooth_z);

            [max_z ind] = max(smooth_z);

            plot(bead(i).time(ind), smooth_z(ind), 'bx', 'linewidth', 2, 'markersize', 20)

            %z_offsets = [z_offsets min_z];

            z_offsets = [z_offsets max_z];

         %%  우리 장치에서는 가장 큰 값이 바닥면에 대한 값으로 추출

            set(gca, 'fontsize', 16, 'linewidth', 1, 'fontweight', 'bold','TickLength',[0.02 0.02]);

            xlabel('Time (s)'); ylabel('z (um)')

            title(['Bead # ' num2str(i)])

            pause;

        end

    end

    %%% Plot the offsets

    if 1

        figure(2); clf; hold on; box on;

        plot(1, z_offsets, 'bo', 'linewidth', 2, 'markersize', 20)

        set(gca, 'fontsize', 16, 'linewidth', 1, 'fontweight', 'bold','TickLength',[0.02 0.02]);

        xlabel('Bead #'); ylabel('z-offset (um)')

    end

    %%% Save the data

    if 1

        foo = [(1)' z_offsets'];

        save(output_name, 'foo', '-ascii')

    end

       %% 바닥면에 대한 정보 추출 하여 저장

end