TailCorr/script_m/z_dsp.m

classdef z_dsp < handle
    properties
        %input
        fs_L;
        fs_H;
        TargetFrequency;
        G;
        simulink_time;
        intp_mode;
        dac_mode_sel;
        route_num;
        env_num;
        
        %output
        Ideal2Low;
        Ideal2Target;
        wave_pre;
        wave_preL;
        amp_real;
        amp_imag;
        time_real;
        time_imag;
        name;
        wave_revised;
        wave_revisedL;
        DownsamplingBy12GDataAlign;
        HardwareMeanIntpDataAlign;
        Delay;
        Delay_mode;
        pause_time;
        filename;
        rpt_num;
        FallingEdge;
        Amp;
        itv_time;    %信号具有周期性时的间隔
    end

    methods

        function obj = z_dsp(fs_L,fs_H,TargetFrequency,G,simulink_time,intp_mode,dac_mode_sel)
            obj.fs_L            = fs_L;
            obj.fs_H            = fs_H;
            obj.TargetFrequency = TargetFrequency;
            obj.G               = G;
            obj.simulink_time   = simulink_time;
            obj.intp_mode       = intp_mode;
            obj.dac_mode_sel    = dac_mode_sel;
            obj.Ideal2Low       = fs_H/(fs_L/2);
            obj.Ideal2Target    = fs_H/TargetFrequency;
            obj.name = [
                "第一组S21参数_flattop_上升沿2ns_持续时间30ns_下降沿",...    
                "第一组S21参数_flattop_上升沿4ns_持续时间30ns_下降沿",...
                "第一组S21参数_flattop_上升沿4ns_持续时间50ns_下降沿",...
                "第一组S21参数_flattop_上升沿4ns_持续时间1000ns_下降沿",...
                "第一组S21参数_flattop_上升沿100ns_持续时间10000ns_下降沿",...
                "第一组S21参数_acz_持续时间30ns_下降沿",...
                "第一组S21参数_acz_持续时间50ns_下降沿";
                "第二组S21参数_flattop_上升沿2ns_持续时间30ns_下降沿",...    
                "第二组S21参数_flattop_上升沿4ns_持续时间30ns_下降沿",...
                "第二组S21参数_flattop_上升沿4ns_持续时间50ns_下降沿",...
                "第二组S21参数_flattop_上升沿4ns_持续时间1000ns_下降沿",...
                "第二组S21参数_flattop_上升沿100ns_持续时间10000ns_下降沿",...
                "第二组S21参数_acz_持续时间30ns_下降沿",...
                "第二组S21参数_acz_持续时间50ns_下降沿";
                "第三组S21参数_flattop_上升沿2ns_持续时间30ns_下降沿",...    
                "第三组S21参数_flattop_上升沿4ns_持续时间30ns_下降沿",...
                "第三组S21参数_flattop_上升沿4ns_持续时间50ns_下降沿",...
                "第三组S21参数_flattop_上升沿4ns_持续时间1000ns_下降沿",...
                "第三组S21参数_flattop_上升沿100ns_持续时间10000ns_下降沿",...
                "第三组S21参数_acz_持续时间30ns_下降沿",...
                "第三组S21参数_acz_持续时间50ns_下降沿";
                "第四组S21参数_flattop_上升沿2ns_持续时间30ns_下降沿",...    
                "第四组S21参数_flattop_上升沿4ns_持续时间30ns_下降沿",...
                "第四组S21参数_flattop_上升沿4ns_持续时间50ns_下降沿",...
                "第四组S21参数_flattop_上升沿4ns_持续时间1000ns_下降沿",...
                "第四组S21参数_flattop_上升沿100ns_持续时间10000ns_下降沿",...
                "第四组S21参数_acz_持续时间30ns_下降沿",...
                "第四组S21参数_acz_持续时间50ns_下降沿";
                "第五组S21参数_flattop_上升沿2ns_持续时间30ns_下降沿",...    
                "第五组S21参数_flattop_上升沿4ns_持续时间30ns_下降沿",...
                "第五组S21参数_flattop_上升沿4ns_持续时间50ns_下降沿",...
                "第五组S21参数_flattop_上升沿4ns_持续时间1000ns_下降沿",...
                "第五组S21参数_flattop_上升沿100ns_持续时间10000ns_下降沿",...
                "第五组S21参数_acz_持续时间30ns_下降沿",...
                "第五组S21参数_acz_持续时间50ns_下降沿";
                "第一组S21参数_flattop_上升沿2ns_持续时间30ns_下降沿后",...    
                "第一组S21参数_flattop_上升沿4ns_持续时间30ns_下降沿后",...
                "第一组S21参数_flattop_上升沿4ns_持续时间50ns_下降沿后",...
                "第一组S21参数_flattop_上升沿4ns_持续时间1000ns_下降沿后",...
                "第一组S21参数_flattop_上升沿100ns_持续时间10000ns_下降沿后",...
                "第一组S21参数_acz_持续时间30ns_下降沿后",...
                "第一组S21参数_acz_持续时间50ns_下降沿后";
                "第二组S21参数_flattop_上升沿2ns_持续时间30ns_下降沿后",...    
                "第二组S21参数_flattop_上升沿4ns_持续时间30ns_下降沿后",...
                "第二组S21参数_flattop_上升沿4ns_持续时间50ns_下降沿后",...
                "第二组S21参数_flattop_上升沿4ns_持续时间1000ns_下降沿后",...
                "第二组S21参数_flattop_上升沿100ns_持续时间10000ns_下降沿后",...
                "第二组S21参数_acz_持续时间30ns_下降沿后",...
                "第二组S21参数_acz_持续时间50ns_下降沿后";
                "第三组S21参数_flattop_上升沿2ns_持续时间30ns_下降沿后",...    
                "第三组S21参数_flattop_上升沿4ns_持续时间30ns_下降沿后",...
                "第三组S21参数_flattop_上升沿4ns_持续时间50ns_下降沿后",...
                "第三组S21参数_flattop_上升沿4ns_持续时间1000ns_下降沿后",...
                "第三组S21参数_flattop_上升沿100ns_持续时间10000ns_下降沿后",...
                "第三组S21参数_acz_持续时间30ns_下降沿后",...
                "第三组S21参数_acz_持续时间50ns_下降沿后";
                "第四组S21参数_flattop_上升沿2ns_持续时间30ns_下降沿后",...    
                "第四组S21参数_flattop_上升沿4ns_持续时间30ns_下降沿后",...
                "第四组S21参数_flattop_上升沿4ns_持续时间50ns_下降沿后",...
                "第四组S21参数_flattop_上升沿4ns_持续时间1000ns_下降沿后",...
                "第四组S21参数_flattop_上升沿100ns_持续时间10000ns_下降沿后",...
                "第四组S21参数_acz_持续时间30ns_下降沿后",...
                "第四组S21参数_acz_持续时间50ns_下降沿后";
                "第五组S21参数_flattop_上升沿2ns_持续时间30ns_下降沿后",...    
                "第五组S21参数_flattop_上升沿4ns_持续时间30ns_下降沿后",...
                "第五组S21参数_flattop_上升沿4ns_持续时间50ns_下降沿后",...
                "第五组S21参数_flattop_上升沿4ns_持续时间1000ns_下降沿后",...
                "第五组S21参数_flattop_上升沿100ns_持续时间10000ns_下降沿后",...
                "第五组S21参数_acz_持续时间30ns_下降沿后",...
                "第五组S21参数_acz_持续时间50ns_下降沿后";
            ];
            obj.pause_time = 0.5;
            obj.Amp = 1.5e4;
        end

        function env(obj)

            cd("D:\Work\EnvData\acz");
            obj1 = py.importlib.import_module('acz');
            py.importlib.reload(obj1);

            %按点数产生理想方波
            % amp_rect = 1.5e4;
            % %单位是ns front是到达时间，flat是持续时间，lagging是后边还有多少个0，会影响脚本的修正时间
            % [front(1), flat(1), lagging(1)] = deal(50,100,7400);% 50,100,7400;100ns方波     
            % [front(2), flat(2), lagging(2)] = deal(50,4000,11500);% 50,4000,11500；4us方波     
            % 
            % for i = 1:2
            %     front_H(i) = front(i)*fs_H/1e9; flat_H(i) = flat(i)*fs_H/1e9; lagging_H(i) = lagging(i)*fs_H/1e9;
            %     wave_pre{i} = amp_rect*cat(2,zeros(1,front_H(i)),ones(1,flat_H(i)),zeros(1,lagging_H(i)));%脚本的单位是点数
            % end
            %flattop波
            A = 1.5e4;
            [edge(1), length_flattop(1)] = deal(2,30);%ns，在fsn_L取1时是参数里的length
            [edge(2), length_flattop(2)] = deal(4,30);
            [edge(3), length_flattop(3)] = deal(4,50);
            [edge(4), length_flattop(4)] = deal(4,1000);
            [edge(5), length_flattop(5)] = deal(100,10000);

            for i = 1:length(length_flattop)
                [edge_H(i), length_H(i)] = deal(edge(i)*obj.fs_H/1e9,length_flattop(i)*obj.fs_H/1e9);
                obj.wave_pre{i} = flattop(A, edge_H(i), length_H(i), 1);
            end
            
            %acz波
            amplitude = 1.5e4;
            
            carrierFreq = 0.000000;
            carrierPhase = 0.000000;
            dragAlpha = 0.000000;
            thf = 0.864;
            thi = 0.05;
            lam2 = -0.18;
            lam3 = 0.04;
            
            length_acz(1) = 30;
            length_acz(2) = 50;
            
            for i = 1:length(length_acz)
                length_acz_H(i) = int32(length_acz(i)*obj.fs_H/1e9);
                obj.wave_pre{i+length(length_flattop)} = real(double(py.acz.aczwave(amplitude, length_acz_H(i), carrierFreq,carrierPhase, dragAlpha,thf, thi, lam2, lam3)));
            end
            
            obj.env_num = length(length_flattop) + length(length_acz);

            for i = 1:obj.env_num
                obj.wave_pre{i} = cat(2,repmat(cat(2,obj.wave_pre{i},zeros(1,round(30e-9*obj.fs_H))),1,obj.rpt_num),zeros(1,floor(obj.simulink_time*obj.fs_H)));    %校正前的高频信号
                obj.wave_preL{i} = obj.wave_pre{i}(1:obj.Ideal2Low:end);    %校正前的低频信号
            end

            assignin("base",'wave_preL',obj.wave_preL);
            obj.FallingEdge = [30e-9,30e-9,50e-9,1000e-9,10000e-9,30e-9,50e-9];
            
        end

        function route(obj)

            obj.amp_real{1}= [0.025 0.015 0.0002 0.2 0 0];
            obj.amp_imag{1}= [0 0 0 0 0 0];
            obj.time_real{1} = [-1/250, -1/650, -1/1600 -1/20 0 0];
            obj.time_imag{1} = [0 0 0 0 0 0];
            
            obj.amp_real{2}= [0.025 0.015 0.0002 0.2 0 0];
            obj.amp_imag{2}= [0 0 0 0 0 0];
            obj.time_real{2} = [-1/250, -1/650, -1/1600 -1/20 0 0];
            obj.time_imag{2} = [0 -1/300 -1/500 0 0 0];
            
            obj.amp_real{3}= [0.025 0.009 0.0002 0.2 0 0];
            obj.amp_imag{3}= [0 0.012 0 0 0 0];
            obj.time_real{3} = [-1/250, -1/650, -1/1600 -1/20 0 0];
            obj.time_imag{3} = [0 -1/300 -1/500 0 0 0];
            
            obj.amp_real{4}= [0.025 0.015 0.0002 0.2 0 0];
            obj.amp_imag{4}= [0 0 0 0 0 0];
            obj.time_real{4} = [-1/250, -1/2000, -1/1600 -1/20 0 0];
            obj.time_imag{4} = [0 -1/15 -1/50 0 0 0];
            
            obj.amp_real{5}= [0.025 0.009 0.0002 0.2 0 0];
            obj.amp_imag{5}= [0 0.012 0 0 0 0];
            obj.time_real{5} = [-1/250, -1/2000, -1/1600 -1/20 0 0];
            obj.time_imag{5} = [0 -1/15 -1/50 0 0 0];
            
            [m,n] = size(obj.amp_real);
            obj.route_num = n;
        end

        function py_cal(obj)

            cd("D:\Work\TailCorr_20241008_NoGit");
            obj2 = py.importlib.import_module('wave_calculation');
            py.importlib.reload(obj2);
            cd("D:\Work\TailCorr");    

            convolve_bound = int8(3);
            calibration_time = int32(20e3);
            cal_method = int8(1);
            sampling_rateL = int64(obj.fs_L/2);
            sampling_rate = int64(obj.fs_H);
            
            %校正后的高频信号
            for m = 1:obj.route_num
                for n = 1:obj.env_num
                    wave_cal = cell(py.wave_calculation.wave_cal(obj.wave_pre{1,n}, obj.amp_real{1,m}, obj.amp_imag{1,m}, obj.time_real{1,m}, obj.time_imag{1,m}, convolve_bound, calibration_time, cal_method, sampling_rate));
                    obj.wave_revised{m,n} = double(wave_cal{1,1});
                    wave_calL = cell(py.wave_calculation.wave_cal(obj.wave_preL{1,n}, obj.amp_real{1,m}, obj.amp_imag{1,m}, obj.time_real{1,m}, obj.time_imag{1,m}, convolve_bound, calibration_time, cal_method, sampling_rateL));
                    obj.wave_revisedL{m,n} = double(wave_calL{1,1});        
                end
                alpha{m} = double(wave_calL{1,2});
                beta{m} =  double(wave_calL{1,3});   
            end            
            alpha_wideth=32; 
            beta_width=32;
            %定点化系数
            for i = 1:obj.route_num
                alphaFixRe{i} = ceil((2^(alpha_wideth-1))*real(alpha{i}));
                alphaFixIm{i} = ceil((2^(alpha_wideth-1))*imag(alpha{i}));
                betaFixRe{i}  = ceil((2^(beta_width-1))*real(beta{i}));
                betaFixIm{i}  = ceil((2^(beta_width-1))*imag(beta{i}));
            end

            assignin('base', 'alphaFixRe', alphaFixRe);
            assignin('base', 'alphaFixIm', alphaFixIm);
            assignin('base', 'betaFixRe' , betaFixRe);
            assignin('base', 'betaFixIm' , betaFixIm);
        end

        function FIL(obj)
            
            for m = 1:obj.route_num
                assignin('base', 'm', m);
                for n = 1:obj.env_num
                    assignin('base', 'n', n);
                    optnons=simset('SrcWorkspace','current');
                    sim('z_dsp_FIL',[0,obj.simulink_time]);
                    sim2m = @(x)reshape(logsout.get(x).Values.Data,[],1);
                    dout0{m,n} = sim2m("dout0");
                    dout1{m,n} = sim2m("dout1");
                    dout2{m,n} = sim2m("dout2");
                    dout3{m,n} = sim2m("dout3");
                
                    N = length(dout0{m,n});
                    cs_wave{m,n} = zeros(4*N,1);            
                    
                    cs_wave{m,n}(1:4:4*N) = dout0{m,n};
                    cs_wave{m,n}(2:4:4*N) = dout1{m,n};
                    cs_wave{m,n}(3:4:4*N) = dout2{m,n};
                    cs_wave{m,n}(4:4:4*N) = dout3{m,n};
                
                    HardwareMeanIntpData{m,n} = cs_wave{m,n};%硬件校正后内插
                    DownsamplingBy12GData{m,n} = obj.wave_revised{m,n}(1:obj.Ideal2Target:end);
                    [obj.DownsamplingBy12GDataAlign{m,n},obj.HardwareMeanIntpDataAlign{m,n},obj.Delay(m,n)] = ...
                    alignsignals(DownsamplingBy12GData{m,n}(1:round(obj.TargetFrequency*20e-6)),HardwareMeanIntpData{m,n}(1:round(obj.TargetFrequency*20e-6)),"Method","xcorr");
                end
            end
            obj.Delay_mode = mode(obj.Delay,'all');
            fprintf('Delay_mode = %d\n',obj.Delay_mode);
        end

        function DataShow(obj,save)

            close all;
                      
            fileID = fopen(obj.filename, 'w');
            if fileID == -1
                disp('文件打开失败');
            else
                disp('文件打开成功');            
            end

            start_time = abs(obj.Delay_mode)/(obj.TargetFrequency/1e9)*1e-9;%由于相位修正后会有偏移的点数，所以需要考虑上这个偏移的时间，采样率为3GHz，3个点对应1ns

            if(obj.rpt_num == 1)
                for m = 1:obj.route_num
                    for n = 1:obj.env_num
                        edge_Align(n) = obj.FallingEdge(n) + start_time;
                        tmp(n) = edge_Align(n) + 10e-9;
                        a{n} = [start_time-5e-9 tmp(n)];%[1/obj.fs_H 50e-9];[50e-9 1.5e-6],[500e-9+10e-9 tmp-20e-9]
                        b{n} = [tmp(n) 20e-6];

                        figure('Units','normalized','Position',[0.0004    0.5174    0.4992    0.4229]);
                        obj.diff_plot_py(obj.TargetFrequency,obj.HardwareMeanIntpDataAlign{m,n}', obj.DownsamplingBy12GDataAlign{m,n}(1:floor(obj.TargetFrequency*20e-6)),obj.name(m,n),'硬件与脚本的差值',a{n},obj.Amp,edge_Align(n),fileID);
                        if(save == "save")
                            savefig(obj.name(m,n));
                        end
                        
                        figure('Units','normalized','Position',[0.0004    0.0340    0.4992    0.4229]);
                        obj.diff_plot_py(obj.TargetFrequency,obj.HardwareMeanIntpDataAlign{m,n}', obj.DownsamplingBy12GDataAlign{m,n}(1:floor(obj.TargetFrequency*20e-6)),obj.name(m+5,n),'硬件与脚本的差值',b{n},obj.Amp,edge_Align(n),fileID);
                        if(save == "save")
                            savefig(obj.name(m+5,n));
                        end
                    end
                end
            else
                for m = 1:obj.route_num
                    for n = 1:obj.env_num
                        figure('Units','normalized','Position',[0    0.0333    1.0000    0.9125]);
                        title(obj.name(m,n),Interpreter="none");
                        tiledlayout('vertical','TileSpacing','tight')
                        obj.diff_plot_py(obj.TargetFrequency,obj.HardwareMeanIntpDataAlign{m,n}', obj.DownsamplingBy12GDataAlign{m,n}(1:floor(obj.TargetFrequency*20e-6)),obj.name(m,n),'硬件与脚本的差值',obj.FallingEdge(n)+obj.itv_time,obj.Amp,start_time,fileID);
                        if(save == "save")
                            savefig(obj.name(m,n));
                        end                    
                    end
                end
            end

            fclose(fileID);

        end
    
        function RouteShow(obj,save)

            t = 0:1/(1e2):10000;
            for i = 1:5
                amp_routing{i}  = obj.amp_real{1,i} + 1j*obj.amp_imag{1,i};
                time_routing{i} = obj.time_real{1,i} + 1j*obj.time_imag{1,i};
                tau{i} = -1./time_routing{i};
            end            

            figure()
            set(gcf,"Position",[1          49        2560        1314])
            tiledlayout('flow','TileSpacing','tight');
            title_name = ["第一组S_{21}参数","第二组S_{21}参数","第三组S_{21}参数","第四组S_{21}参数","第五组S_{21}参数"];
            for m = 1:obj.route_num
                for n = 1:1:length(amp_routing{1,m})
                    S21_time{m}(:,n) = amp_routing{1,m}(n)*exp(time_routing{1,m}(n)*t);
                end
                nexttile
                plot(t*1e-9,real(sum(S21_time{m},2)));
                grid on
                title(title_name(m));
            end

            if(save == "save")
                savefig("S21线路参数");
            end

        end
        
        function FigDisplay(obj)
            if(obj.rpt_num == 1)
                for m = 1:obj.route_num*obj.env_num
                    figure(2*m-1)
                    figure(2*m)
                    pause(obj.pause_time);
                end

            else
                for m = 1:obj.route_num*obj.env_num
                    figure(m)
                    pause(obj.pause_time);
                end                
            end
        end

        function LoadFigAndDisplay(obj)
            for n = 1:obj.route_num
                for m = 1:obj.env_num
                open(strcat(obj.name(n,m),'.fig'));
                open(strcat(obj.name(n+5,m),'.fig'));
                pause(obj.pause_time);
                end
            end                        
        end

        function ErrAny(obj,save)
            fid = fopen(obj.filename,'r');
            if(obj.rpt_num == 1)
                data = textscan(fid,'Falling edge of 20ns~40ns mean :%s	 std :%s	Falling edge of 1us~1.1us mean :%s	 std :%s	The mean and std stably less than 1e-4 is :%s s');
                fclose(fid);
                data{1} = cellfun(@str2num,data{1});
                data{2} = cellfun(@str2num,data{2});
                data{3} = cellfun(@str2num,data{3});
                data{4} = cellfun(@str2num,data{4});
                data{5} = cellfun(@str2num,data{5});
                title_name = ["下降沿后20ns~40ns误差的平均值","下降沿后20ns~40ns误差的标准差","下降沿后1us~1.1us误差的平均值","下降沿后1us~1.1us误差的标准差","加窗参数"];    
                err_threshold = [1e-3 1e-3 1e-4 3e-4 5e-6];
            else
                data = textscan(fid,'每个周期拖尾误差均值的标准差 = %s s');
                fclose(fid);
                data{1} = cellfun(@str2num,data{1});
                title_name = ["多周期误差平均值的标准差"];
                err_threshold = [0.5e-3];                
            end
            [h,v] = size(data);
            figure()
            tiledlayout('flow','TileSpacing','tight')
            colors = lines(obj.route_num);   
            set(gcf,'Position', [1          49        2560        1314]);            
            for m = 1:v
                nexttile
                hold on
                for i = 1:(obj.route_num)
                    idx = (i-1)*(length(data{m})/obj.route_num) + 1 : i*(length(data{m})/obj.route_num);
                    plot(idx,abs(data{m}(idx)),'-o','Color', colors(i, :));
                end
                yline(err_threshold(m),'--r');
                title(title_name(m));
                set(gca,'YScale','log');
                legend("第一组线路","第二组线路","第三组线路","第四组线路","第五组线路",'Location','northwest');
            end
            if(obj.rpt_num == 1)
                if(save == "save")
                    savefig("单周期误差分析")
                end
            else
                if(save == "save")
                    savefig("多周期误差分析")
                end
            end

        end

        %compare FIL with python script
        function diff_plot_py(obj,fs,iir_out, Script_out,title1,title2,a,amp,edge,fileID)
            
            %输入数据长度不等时取其公共部分
            N = min(length(iir_out),length(Script_out));
            iir_out = iir_out(1:N);
            Script_out = Script_out(1:N);
            
            diff = (iir_out - Script_out)/amp;%求差，并归一化
            
            n = (0:1:N-1)/fs;
            %找出关心的数据点
            if(obj.rpt_num == 1)
                n_edge = find(n>=edge-1e-12);%edge代表下降沿
                n50 = find(n>=edge+20e-9-1e-12);%下降沿后20ns
                n20_40 = find((n>=edge+20e-9-1e-12) & (n<=edge+40e-9+1e-12));%下降沿后20ns到40ns
                n1000 = find(n>=edge+1000e-9-1e-12);%下降沿后1us
                n1000_1100 = find((n>=edge+1000e-9-1e-12) & (n<=edge+1100e-9+1e-12));%下降沿后1us到1.1us
                
                ne = find((abs(diff)>=1e-4) & (abs(diff)<1));%误差小于万分之一的点
                ne(1) = 1;
                
                window_length = 100e-9*fs;
                diff_mean_window = movmean(diff,window_length);
                diff_std_window = movstd(diff,window_length);
                n_mean_window = find((abs(diff_mean_window)>=1e-4) );%100ns窗，误差均值小于万分之一点
                n_std_window = find((abs(diff_std_window)>=1e-4) ); %100ns窗，误差方差小于万分之一点
                n_common = max(n_mean_window(end),n_std_window(end));
                %原始数据作图
                tiledlayout(2,1)
                ax1 = nexttile;
                plot(n,iir_out,n,Script_out)
                legend('硬件','软件');
                xlabel('t/s')
                xlim(a);
                title(title1,Interpreter="none");
                grid on
                hold on
                
                %差值做图
                ax2 = nexttile;
                plot(n,diff)
                xlabel('t/s')
                title('diff')
                grid on
                hold on
                xlim(a)
                title('硬件与脚本的差值',Interpreter="none");
                linkaxes([ax1,ax2],'x');
                
                plot_p = @(x)[
                    plot(n(x),diff(x),'r*');
                    text(n(x), diff(x)+diff(x)*0.1, ['(',num2str(n(x)),',',num2str(diff(x)),')'],'color','k');
                    ];
                
                ne(1) = 1;
                
                % [diff_max,R_mpos] = max(abs(diff));%误差最大值
                % plot_p(R_mpos);
                
                if a(2) <= 5e-6
                    plot_p(n_edge(1));%下降沿
                    % plot_p(R_mpos);
                elseif a(2) == 20e-6
                    plot_p(n50(1));   %下降沿20ns
                    plot_p(n1000(1)); %下降沿1us
                    plot_p(ne(end));  %误差小于万分之一
                    fprintf(fileID,"Falling edge of 20ns~40ns mean :%.4e\t std :%.4e\t",mean(diff(n20_40)),std(diff(n20_40)));
                    fprintf(fileID,"Falling edge of 1us~1.1us mean :%.4e\t std :%.4e\t",mean(diff(n1000_1100)),std(diff(n1000_1100)));
                    % fprintf("The error after falling edge of 1us is:%.4e\t",diff(n1000(1)));
                    % fprintf("The time of erroe less than 1e-4 is :%.4e us\n",(n(ne(end))-n(n_edge(1))));
                    fprintf(fileID,"The mean and std stably less than 1e-4 is :%.4e s\n",(n(n_common)-n(n_edge(1))));
                end 
            else
                n_start = find(n>=edge-1e-12);%edge代表下降沿
                
                % 确定周期长度对应的采样点数量
                T = a;    %在这种情况下，a这个参数用不到了，使用其传递周期，也就是说a这个参数有两种不同的涵义
                samples_per_period = round(T * fs);  % 每个周期采样点数
                num_periods = obj.rpt_num; % 总周期数
                period_means = zeros(1, num_periods); % 存储每周期均值

                    for i = 1:num_periods
                        % 提取当前周期的起止索引
                        start_idx(i) = n_start(1) + (i - 1) * samples_per_period;
                        end_idx(i) = n_start(1) + i * samples_per_period;
        
                        % 提取当前周期的数据
                        period_data = diff(start_idx(i):end_idx(i));
        
                        % 计算当前周期的均值
                        period_means(i) = mean(period_data);
                    end
                fprintf(fileID,"每个周期拖尾误差均值的标准差 = %.4e s\n",std(period_means));            
                ax1 = nexttile;
                plot(n,iir_out,n,Script_out);
                hold on
                plot(n(start_idx), Script_out(start_idx), 'r*'); % 标记每个周期的起始点            
                plot(n(end_idx), Script_out(end_idx), 'g*'); % 标记每个周期的起始点      
                legend('硬件','软件');
                xlabel('t/s');
                title(title1,Interpreter="none");
                ax2 = nexttile;
                hold on
                plot(n, diff); hold on;                     % 原始信号
                plot(n(end_idx), diff(end_idx), 'g*'); % 标记每个周期的起始点     
                xlabel('t/s');
                ylabel('归一化误差');
                linkaxes([ax1,ax2],'x');
                xlim([0,n(end_idx(end)) + 5e-7]);
                title(title2,Interpreter="none");
                
            end
        end

    
        
    end

end