TailCorr/script_m/z_dsp.m

565 lines
27 KiB
Matlab
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

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,115004us方波
%
% 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;%由于相位修正后会有偏移的点数所以需要考虑上这个偏移的时间采样率为3GHz3个点对应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