clc;clear;close all

% hdlsetuptoolpath('ToolName','Xilinx Vivado','ToolPath','D:\SoftWare\Xilinx\Vivado\2019.2\bin\vivado.bat');	

addpath(genpath('D:\Work\EnvData'));	
addpath(genpath('D:\Work\EnvData\data-v2'));	
% addpath(genpath('D:\Work\TailCorr_20241008_NoGit'));	

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

fs_L = 0.75e9;    %硬件频率
fs_H = 12e9;      %以高频近似理想信号
TargetFrequency = 3e9;
Ideal2Low = fs_H/(fs_L/2);
Ideal2Target = fs_H/TargetFrequency;
G = 1;
DownSample = 2;
simulink_time = 20e-6;  %1.5*16e-6;1.5e-3
intp_mode = 3;     %0不内插，1内插2倍，2内插4倍,3内插8倍
dac_mode_sel = 0;  %选择DAC模式，0出八路，1邻近插值，2邻近插值

%按点数产生理想方波
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(6,50);

for i = 1:4
    [edge_H(i), length_H(i)] = deal(edge(i)*fs_H/1e9,length_flattop(i)*fs_H/1e9);
    wave_pre{i+2} = 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:2
    length_acz_H(i) = int32(length_acz(i)*fs_H/1e9);
    wave_pre{i+6} = real(double(py.acz.aczwave(amplitude, length_acz_H(i), carrierFreq,carrierPhase, dragAlpha,thf, thi, lam2, lam3)));
end
% signalAnalyzer(wave_pre{2},'SampleRate',fs_H);

for i = 1:8
    wave_pre{i} = cat(2,wave_pre{i},zeros(1,floor(simulink_time*fs_H)));    %校正前的高频信号
    wave_preL{i} = wave_pre{i}(1:Ideal2Low:end);    %校正前的低频信号
end

% signalAnalyzer(HardwareMeanIntpDataAlign{1},'SampleRate',3e9);

%%%python脚本校正结果
%S21参数
amp_real = [0.025 0.015 0.0002 0 0 0];
amp_imag = [0 0 0 0 0 0];
time_real = [-1/250, -1/650, -1/1600 0 0 0];
time_imag = [0 0 0 0 0 0];

% amp_real = [0.0539981,-0.0319766,0.084015161,0.0048,0,0];
% amp_imag = [0,-0.041014189,-0.052936266,0,0,0];
% time_real = [-0.0024820146,-0.0080529118,-0.006728925,-0.0001,0,0];
% time_imag = [0,-0.008137675,-0.0033212836,0,0,0];
% 
% amp_real = [0.025 0.015 0.0002 0.2 0 0];
% amp_imag = [0 0 0 0 0 0];
% time_real = [-1/250, -1/650, -1/1600 -1/20 0 0];
% time_imag = [0 0 0 0 0 0];

amp_routing = amp_real + 1j*amp_imag;
time_routing = time_real + 1j*time_imag;
tau = -1./time_routing;

convolve_bound = int8(3);
calibration_time = int32(20e3);
cal_method = int8(1);
sampling_rateL = int64(fs_L/2);
sampling_rate = int64(fs_H);

%校正后的高频信号
for i = 1:8
    wave_cal = cell(py.wave_calculation.wave_cal(wave_pre{i}, amp_real, amp_imag, time_real, time_imag, convolve_bound, calibration_time, cal_method, sampling_rate));
    wave_revised{i} = double(wave_cal{1,1});
    wave_calL = cell(py.wave_calculation.wave_cal(wave_preL{i}, amp_real, amp_imag, time_real, time_imag, convolve_bound, calibration_time, cal_method, sampling_rateL));
    wave_revisedL{i} = double(wave_calL{1,1});
end
%校正后的低频信号
alpha = double(wave_calL{1,2});
beta =  double(wave_calL{1,3});
beta(5:6) = 0;
alpha_wideth=32; 
beta_width=32;
alphaFixRe = ceil((2^(alpha_wideth-1))*real(alpha));
alphaFixIm = ceil((2^(alpha_wideth-1))*imag(alpha));
betaFixRe  = ceil((2^(beta_width-1))*real(beta));
betaFixIm  = ceil((2^(beta_width-1))*imag(beta));

%%%仿真
for i = 1:8
    options=simset('SrcWorkspace','current');
    sim('z_dsp',[0,simulink_time]);
    sim2m = @(x)reshape(logsout.get(x).Values.Data,[],1);
    dout0{i} = sim2m("dout0");
    dout1{i} = sim2m("dout1");
    dout2{i} = sim2m("dout2");
    dout3{i} = sim2m("dout3");

    N(i) = length(dout0{i});
    cs_wave{i} = zeros(4*N(i),1);            
    
    cs_wave{i}(1:4:4*N) = dout0{i};
    cs_wave{i}(2:4:4*N) = dout1{i};
    cs_wave{i}(3:4:4*N) = dout2{i};
    cs_wave{i}(4:4:4*N) = dout3{i};

    HardwareMeanIntpData{i} = cs_wave{i};%硬件校正后内插
    DownsamplingBy12GData{i} = wave_revised{i}(1:Ideal2Target:end);
    [DownsamplingBy12GDataAlign{i},HardwareMeanIntpDataAlign{i},Delay(i)] = ...
        alignsignals(DownsamplingBy12GData{i}(1:round(TargetFrequency*20e-6)),HardwareMeanIntpData{i}(1:round(TargetFrequency*20e-6)),"Method","xcorr");
end

% signalAnalyzer(DownsamplingBy12GDataAlign{1},HardwareMeanIntpDataAlign{1},'SampleRate',3e9);

%% 绘图并保存
close all;
Amp = 1.5e4;
FallingEdge = [
    150e-9,4050e-9,...%矩形波
    30e-9,30e-9,50e-9,50e-9,...%flattop
    30e-9,50e-9%acz
    ];
name = [
    "rect_100ns_校正后的波形_下降沿后10ns.fig","rect_100ns_校正后的波形_下降沿后1us.fig";
    "rect_4us_校正后的波形_下降沿后10ns.fig","rect_4us_校正后的波形_下降沿后1us.fig";   
    "flattop_上升沿2ns_持续时间30ns_校正后的波形_下降沿后10ns.fig","flattop_上升沿2ns_持续时间30ns_校正后的波形_下降沿后1us.fig";   
    "flattop_上升沿4ns_持续时间30ns_校正后的波形_下降沿后10ns.fig","flattop_上升沿4ns_持续时间30ns_校正后的波形_下降沿后1us.fig";   
    "flattop_上升沿4ns_持续时间50ns_校正后的波形_下降沿后10ns.fig","flattop_上升沿4ns_持续时间50ns_校正后的波形_下降沿后1us.fig";  
    "flattop_上升沿6ns_持续时间50ns_校正后的波形_下降沿后10ns.fig","flattop_上升沿6ns_持续时间50ns_校正后的波形_下降沿后1us.fig";
    "acz_持续时间30ns_校正后的波形_下降沿后10ns.fig","acz_持续时间30ns_校正后的波形_下降沿后1us.fig";
    "acz_持续时间50ns_校正后的波形_下降沿后10ns.fig","acz_持续时间50ns_校正后的波形_下降沿后1us.fig";
];
Delay_mode = mode(Delay);
for i = 1:8
    start_time(i) = abs(Delay_mode)/(TargetFrequency/1e9)*1e-9;%由于相位修正后会有偏移的点数，所以需要考虑上这个偏移的时间，采样率为3GHz，3个点对应1ns
    edge_Align(i) = FallingEdge(i) + start_time(i);
    tmp(i) = edge_Align(i) + 10e-9;
    a{i} = [start_time(i)-5e-9 tmp(i)];%[1/fs_H 50e-9];[50e-9 1.5e-6],[500e-9+10e-9 tmp-20e-9]
    b{i} = [tmp(i) 10e-6];
    fig1 = figure('Units','normalized','Position',[0.000390625,0.517361111111111,0.49921875,0.422916666666667]);
    diff_plot_py(TargetFrequency,HardwareMeanIntpDataAlign{i}', DownsamplingBy12GDataAlign{i}(1:floor(TargetFrequency*20e-6)),'HardwareRevised','ScriptRevised',a{i},Amp,edge_Align(i));
    title(name(i,1),Interpreter="none");
    savefig(name(i,1));
    fig2 = figure('Units','normalized','Position',[0.000390625,0.034027777777778,0.49921875,0.422916666666667]);
    diff_plot_py(TargetFrequency,HardwareMeanIntpDataAlign{i}', DownsamplingBy12GDataAlign{i}(1:floor(TargetFrequency*20e-6)),'HardwareRevised','ScriptRevised',b{i},Amp,edge_Align(i));
    title(name(i,2),Interpreter="none");
    savefig(name(i,2));
end

%% 可视化S21参数
t = 0:1/(1e2):10000;

for i = 1:1:length(amp_routing)
    S21_time(:,i) = amp_routing(i)*exp(time_routing(i)*t);
end

figure
plot(t*1e-9,real(sum(S21_time,2)));
grid on
title("s(t)");
savefig("S21参数");

% signalAnalyzer(real(sum(S21_time,2)),'SampleRate',1e11);%时间是1ns，还得加上采样率

% rmpath(genpath('D:\Work\EnvData'));
% rmpath(genpath('D:\Work\EnvData\data-v2'));
% rmpath(genpath('D:\Work\TailCorr_20241008_NoGit'));
%% 图像可视化
cd("D:\Work\TailCorr\仿真结果\20241101_125M八倍内插至1G_第1组S21参数")
for i = 1:8
    close all
    open(name(i,1));
    open(name(i,2));
    pause()
end