超级会员免费看
特定Matlab脚本的使用与分析
在实验研究中，为了控制实验结果的不确定性，往往需要采集多组数据。下面将详细介绍如何使用不同的Matlab脚本来实现从示波器采集多组数据、进行数值积分以及蒙特卡罗模拟等操作。
1. 从两台示波器采集多组数据
为了控制实验结果的不确定性，需要采集多组数据。这里使用两台示波器分别采集来自两个相机和信号发生器的多组数据。由于相机长时间运行会发热，因此每个相机大约每隔30分钟进行一次数据采集，采集间隙相机处于关闭冷却状态。
1.1 DSOS604A高速示波器
DSOS604A是一款中速示波器，可实现信号的自动连续采集。以下是采集和存储数据的操作步骤：
1. 选择“Trigger”选项卡。
2. 点击“Trigger Action …”。
3. 勾选“Perform Multipurpose On Trigger”旁边的框。
4. 在“Max No Of Actions”下设置每次触发需要保存的次数。
5. 点击“Multipurpose On Trigger”部分的“Setup …”。
6. 在“Multipurpose”下拉列表中选择“QuickWaveform”。
7. 设置文件格式和保存位置，文件将以逗号分隔值（*.csv）格式保存。
1.2 HP Infiniium 54846B低速示波器
HP Infiniium 54846B是一款低速示波器，可通过USB/GBIP接口连接，并使用以下Matlab脚本连续采集多个信号：
% 03 May , 2017
% To collect the data from HP infinnium 54846B oscilloscope via USB/ GBIP interface:
% Capturing 65536 points for SR4000 / DS325 / SigGen 8648B using the oscilloscope
clc %Clear the MATLAB command window
clear %Clear the MATLAB variables
% Setup :
instrumentConnected = true;
% Ignore the instrument when doing a test:
if instrumentConnected
% Instrument
visaObjOscope = visa('agilent', 'GPIB0::7::INSTR');
% Set the buffer size
visaObjOscope.InputBufferSize = 100000;
% Set the timeout value
visaObjOscope.Timeout = 10;
% Set the Byte order
visaObjOscope.ByteOrder = 'littleEndian';
% Open the connection
fopen(visaObjOscope);
end
for k = 501:1050 % Use 500 acquisition slots (taken ~30 mins), since the camera is bit heating.
if instrumentConnected
fprintf(visaObjOscope, ':DISK:STOR CHAN1, "F:\\SR400\\SR30_8GSa%d.csv", TEXT, XYP, ON', k);
end
pause(4); % About 3.3 sec needs to store each file for SR4000 & DS325
end
if instrumentConnected
% Check for any instrument errors
instrumentError = query(visaObjOscope, ':SYST:ERR?');
% Close the instrument
fclose(visaObjOscope);
end
这个脚本的执行流程可以用以下mermaid流程图表示：
graph TD;
A[开始] --> B[清除命令窗口和变量];
B --> C{仪器是否连接};
C -- 是 --> D[初始化仪器连接];
C -- 否 --> E[跳过仪器操作];
D --> F[设置仪器参数];
F --> G[打开仪器连接];
G --> H[循环采集数据];
H --> I[保存数据];
I --> J[暂停4秒];
J --> H;
H --> K{采集完成};
K -- 是 --> L[检查仪器错误];
L --> M[关闭仪器连接];
M --> N[结束];
K -- 否 --> H;
E --> N;
2. Romberg积分算法
Romberg积分算法是一种数值积分方法，用于计算函数在给定区间内的积分。以下是实现该算法的Matlab函数：
function [R] = rombf(f, a, b, N)
% 18 June , 2019
% Romberg method is used for computing the integral of function f in the range [a, b] using
% at most k Richardson extrapolations. Stop when q extrapolations have been performed.
% Inputs:
% f - function,
% a, b - lower and upper bound of the integral
% N - Romberg table N by N with 2^(N−1) panels.
% Output:
% R - Romberg table with R[N, N] is the best estimate of the value of the integral.
h = b - a;
np = 1;
% First term R(1, 1)
R(1, 1) = h * (feval(f, a) + feval(f, b)) / 2;
for i = 2:N
h = h / 2;
np = 2 * np;
sum = 0; % Compute midpoint sum
for k = 1 : 2 : np - 1
sum = sum + feval(f, a + k * h);
end
% First column entry R(i, 1)
R(i, 1) = R(i - 1, 1) / 2 + h * sum;
m = 1;
% Other columns R(i, 2), ..., R(i, i)
for k = 2:i % Repeated Richardson extrapolation
m = 4 * m;
R(i, k) = R(i, k - 1) + (R(i, k - 1) - R(i - 1, k - 1)) / (m - 1);
end
end
end
该算法的主要步骤如下：
1. 初始化步长
h
和面板数
np
。
2. 计算第一项
R(1, 1)
。
3. 循环计算其他列的值，通过重复的Richardson外推法不断提高积分的精度。
3. 蒙特卡罗模拟
蒙特卡罗模拟用于评估随机抖动对距离测量的影响，该部分包含两个文件：一个用于使用蒙特卡罗方法生成数据并通过非参数估计获得概率密度函数（PDF），另一个用于生成相位误差结果并分析结果的不确定性。
3.1 生成数据的脚本：genMCsRJ.m
% July 04 , 2019
% This is for analyse the effect of the random jitter (RJ)
% (i.e., epsS & epsL (random variates)) on range measurements by frequency domain.
% Light function is funL = cos(2*pi*frqL*(t + epsL) - phiToF)
% Shutter function is funS = cos(2*pi*frqS*(t + epsS) + theta)
clear all;
nTrpz = 60e6; % number of steps in 'trapz' integration
nEval = 500; % The number of function evaluations
nSet = 30; % Repeating the number of set for uncertainty calculation
setNo = 1: nSet; % for x axis plot
T = 1e-3; % integration period in sec (FIXED)
phiToF = (0 : 6) * pi / 6; % ToF phase shift in rad
rad2deg(phiToF);
frqL = 600e6; % mod frequency of light signal --> [30 50 200 ... 1000] MHz
frqS = frqL; % mod frequency of shutter signal --> [30 50 200 ... 1000] MHz
% Random jitter parameters
muL = 0e-12; % mean in sec
sgmL = (0 : 50 : 1000) * 1e-12; %[0 50 100 200 400 ... 1000]*1e-12; % std dev in sec
muS = 0e-12; % mean in sec
sgmS = [0 sgmL(2 : end) + 1e-12]; %[0 51 101 201 401 ... 1001]*1e-12; % std dev in sec
tt = linspace(-T / 2, T / 2, nTrpz);
% Run the correlation function
ifRJ = 1; % 1 for RJ
for ii = 2:3 % 0 30 60 ... degrees
tic; % starting the time
for nL = 2:2 % 0 50 100 ... ps
for nS = nL %2:2:4% length(sgmS)
for nSt = 1: nSet
parfor nEvl = 1: nEval
esL = muL + randn(nTrpz, 1) * sgmL(nL);
esS = muS + randn(nTrpz, 1) * sgmS(nS);
if ifRJ == 0
funI = @(tt) exp(2i * pi * frqS * tt) .* cos(2 * pi * frqL * tt - phiToF(ii));
else
funI = exp(2i * pi * frqS * (tt + esS')) .* cos(2 * pi * frqL * (tt + esL') - phiToF(ii));
end
angIone = angle(pi * trapz(tt, funI'));
errAngPhi(nEvl) = abs(angIone) - abs(phiToF(ii));
end
stdErrAngPhi(nSt) = std(errAngPhi);
% Plotting the density estimations
optBWhist = 2 / nEval ^ (1 / 3) * iqr(errAngPhi);
if sgmL(nL) == 0
noBinsHist = 1;
else
noBinsHist = floor((max(errAngPhi) - min(errAngPhi)) / optBWhist);
end
optBWkde = 1.06 / nEval ^ (1 / 5) * min(stdErrAngPhi(nSt), iqr(errAngPhi) / 1.34);
xPos = min(errAngPhi); %min
yPosH = max(hist(errAngPhi, noBinsHist));
yPosK = max(kde(errAngPhi));
figure;
subplot(1, 2, 1);
hh = histfit(errAngPhi, noBinsHist);
hh(1).FaceColor = 'b';
text(xPos, yPosH, strcat('\sigma_{\phi_{err}} = ', num2str(stdErrAngPhi)), 'FontSize', 14);
title(['Histogram for the PDF of phase error \phi_{err} for \phi = ', num2str(phiToF(ii)), ' rad']);
xlabel('Phase error \phi_{err} (rad)');
ylabel('Frequency count');
set(gca, 'FontSize', 20);
[yVal, xVal] = kde(errAngPhi, 'width', optBWkde);
subplot(1, 2, 2);
plot(xVal, yVal, 'LineWidth', 1.5, 'Color', 'b');
text(xPos, yPosK, strcat('\sigma_{\phi_{err}} = ', num2str(stdErrAngPhi)), 'FontSize', 14);
title(['KDE for the PDF of phase error \phi_{err} for \phi = ', num2str(phiToF(ii)), ' rad']);
xlabel('Phase error \phi_{err} (rad)');
ylabel('Estimated density function for \phi_{err}');
set(gca, 'FontSize', 20);
% Calculate the uncertainty
muErrAngPhi(nSt) = mean(errAngPhi);
varErrAngPhi(nSt) = var(errAngPhi);
end
end
meanStdErrAngPhi1(nL) = mean(stdErrAngPhi);
uncStdErrAngPhi1(nL) = std(stdErrAngPhi) / sqrt(length(stdErrAngPhi));
meanStdErrAngPhi(nL) = sqrt(sum(varErrAngPhi) / length(varErrAngPhi));
uncStdErrAngPhi(nL) = sqrt(sum((stdErrAngPhi - sqrt(sum(varErrAngPhi) / length(varErrAngPhi))).^2)) / length(varErrAngPhi);
% Plotting std. dev. of the phase error by nSet wise
figure;
plot(setNo, stdErrAngPhi{nL}, '.-');
grid on;
title(strcat('Standard deviation of the phase error \phi_{err} for the number of sets: ', num2str(nSet), ' for \phi = ', num2str(phiToF(ii)), ' rad when f_{mod} = ', num2str(fl / 1e6), ' MHz, \sigma_l = ', num2str(sgmL(nL) * 1e12), ' ps and \sigma_s = ', num2str(sgmS(nS) * 1e12), ' ps'));
xlabel('Set number');
ylabel('Standard deviation of phase error \sigma_{\phi_{err}}');
text(max(setNo) / 2, max(stdErrAngPhi{nL}), strcat('Standard deviation = ', num2str(meanStdErrAngPhi(nL)), ' \pm ', num2str(uncStdErrAngPhi(nL))));
end
endT = toc; % completing time
strT = strcat(num2str(endT / 60), ' mins');
valPhiToF = phiToF(ii);
% Plotting uncertainty of the std. dev. of the phase error by sigL wise
figure;
plot(sgmL * 1e12, meanStdErrAngPhi, 'o', 'LineWidth', 1, 'MarkerSize', 10, 'color', 'b');
grid on;
title(strcat('Uncertainty of the standard deviation of the phase error \phi_{err} by sigL wise with sets: ', num2str(nSet), ' for \phi = ', num2str(valPhiToF), ' rad when f_{mod} = ', num2str(frqL / 1e6), ' MHz'));
xlabel('RMS of the injected random jitter \sigma_{RJ} (ps)');
ylabel('Typical phase error $\overline{\sigma_{\phi_{err}}}$ (rad)', 'interpreter', 'latex');
set(gca, 'FontSize', 20);
xlim([0 1050]);
text(0, 0, strT);
% Save the work space as a .mat file
filename = strcat('phaErrRJ', num2str(rad2deg(valPhiToF)), ' deg ', num2str(frqL / 1e6), 'MHz', num2str(nSet), ' sets.mat');
dataDir = 'E:\GA\Matlab\JitterAnalysis\wavesSR\jitsCalculation\SR4000\mcsRJ\';
cd(dataDir);
save(filename, 'nSet', 'frqL', 'sgmL', 'sgmS', 'valPhiToF', 'stdErrAngPhi', 'meanStdErrAngPhi', 'uncStdErrAngPhi', 'strT', '-v7.3');
end
该脚本的主要步骤如下：
1. 初始化参数，包括积分步数、函数评估次数、重复次数等。
2. 定义随机抖动参数，如均值和标准差。
3. 循环计算不同相位偏移和随机抖动下的相位误差。
4. 绘制相位误差的直方图和核密度估计图。
5. 计算相位误差的标准差和不确定性。
6. 绘制标准差和不确定性的图形。
7. 保存工作空间为.mat文件。
3.2 分析结果不确定性的脚本：rstMCsRJ.m
% Modified on September 12 , 2019
% This is for -->
% (1) to plot the phase error for all frequencies on same graph
% (2) to analyse the uncertainty of the resulted phase individually
clear all;
close all;
phaToF = [30 60]; % ToF phase shift
modFrq = [30 50 100 200 400 600 800 1000]; % Modulation frequencies in MHz
dataDir = 'E:\GA\Matlab\JitterAnalysis\wavesSR\jitsCalculation\SR4000\mcsRJ\';
cd(dataDir);
symMrk = {'∗−c', '+−k', '^−m', 'd−g', 's−b', 'o−y', '.−r', 'x−k'};
symMrkUnc = {'∗c', '+k', ' ^m', ' dg', ' sy', ' ob', ' .r', ' xk'};
for pp = 1: length(phaToF)
matName = strcat('phaErrRJ', num2str(phaToF(pp)), ' ∗30 sets.mat');
matFiles = dir(matName);
% (1) to plot the phase error
figure;
for qq = 1: length(matFiles)
matNameFrqWise = strcat('phaErrRJ', num2str(phaToF(pp)), ' deg ', num2str(modFrq(qq)), ' MHz30sets');
load(matNameFrqWise);
meanStdErr(qq, :) = meanStdErrAngPhi;
uncStdErr(qq, :) = uncStdErrAngPhi;
rmsRJ(qq, :) = sgmL;
modFrq(qq) = frqL;
semilogy(rmsRJ(qq, :) * 1e12, meanStdErr(qq, :), symMrk{qq}, 'LineWidth', 1, 'MarkerSize', 10, 'color', symMrk{qq}(end));
strLegNum{qq} = num2str(frqL / 1e6);
hold on;
end
hanLegNum = legend(strLegNum, 'Location', 'northwest', 'FontSize', 20, 'interpreter', 'latex');
legNumTtl = get(hanLegNum, 'Title');
set(legNumTtl, 'String', '{\it f} (MHz)', 'FontWeight', 'normal', 'FontSize', 20, 'interpreter', 'latex');
xlabel('RMS of the injected random jitter $\sigma_{RJ}$ (ps)', 'interpreter', 'latex');
ylabel('Typical phase error $\overline{\sigma_{\phi_{err}}}$ (rad)', 'interpreter', 'latex');
set(gca, 'FontSize', 20);
xlim([0 1050]);
hold off;
% (2) to analyse the uncertainty
for qq = 1: length(matFiles)
figure;
rmsVals = [rmsRJ(qq, :) * 1e12, rmsRJ(qq, :) * 1e12];
uncVals = [uncStdErr(qq, :), -uncStdErr(qq, :)];
stem(rmsVals, uncVals, symMrk{qq}, 'LineWidth', 1, 'MarkerSize', 10, 'color', symMrk{qq}(end));
set(gca, 'yscal', 'log');
grid on;
xlabel('RMS of the injected random jitter $\sigma_{RJ}$ (ps)', 'interpreter', 'latex');
ylabel('Error bars for phase error $\overline{\sigma_{\phi_{err}}}$ (rad)', 'interpreter', 'latex');
set(gca, 'FontSize', 25);
xlim([0 1050]);
title(strcat('Uncertainty of the mean of standard deviation of the phase error $\overline{\sigma_{\phi_{err}}}$ (rad)', num2str(nSet), ' sets for \phi = ', num2str(valPhiToF), ' rad when f_{mod} = ', strLegNum{qq}, ' MHz'), 'FontSize', 10);
end
% Plot the corresponding 0.1 mrad phase error for RJ against Mod freq.
modFrq = [30 50 100 200 400 600 800 1000]; %in MHz
sgmRJ = [530, 410, 288, 203, 142, 115, 97, 87]; % in ps
figure;
plot(modFrq, sgmRJ, ' ok', 'MarkerSize', 10);
hold on;
legend('RMS \sigma_{RJ}', 'Fitted curve');
xlabel('Modulation frequency $f$ (MHz)', 'interpreter', 'latex');
ylabel('Random jitter of RMS $\sigma_{RJ}$ (ps)', 'interpreter', 'latex');
set(gca, 'FontSize', 25, 'FontName', 'Times New Roman');
xlim([0 1050]);
ylim([0 550]);
end
该脚本的主要功能是：
1. 绘制所有频率下的相位误差在同一图形上。
2. 分析每个频率下相位误差结果的不确定性，并绘制相应的图形。
3. 绘制对应0.1 mrad相位误差的随机抖动与调制频率的关系图。
通过以上的Matlab脚本和操作步骤，可以实现从示波器采集多组数据、进行数值积分以及蒙特卡罗模拟等功能，从而更好地控制实验结果的不确定性，并对实验数据进行深入分析。
特定Matlab脚本的使用与分析
4. 脚本使用总结与注意事项
4.1 数据采集脚本总结
- DSOS604A示波器 ：操作相对直观，通过在示波器界面上进行一系列设置即可实现自动连续采集数据。在设置过程中，需要注意文件格式和保存位置的选择，确保数据能够正确保存。
- HP Infiniium 54846B示波器 ：需要使用Matlab脚本通过USB/GBIP接口进行连接和数据采集。在运行脚本前，要确保仪器连接正常，并且检查脚本中的仪器地址是否正确。同时，由于相机发热问题，设置了采集间隙，在实际应用中可以根据相机的散热情况调整采集间隔时间。
| 示波器型号 | 采集方式 | 注意事项 |
|---|---|---|
| DSOS604A | 示波器界面设置 | 文件格式和保存位置选择 |
| HP Infiniium 54846B | Matlab脚本连接 | 仪器连接和地址检查，采集间隔调整 |
4.2 Romberg积分算法脚本总结
Romberg积分算法通过不断细分区间和使用Richardson外推法来提高积分的精度。在使用该脚本时，需要注意输入参数的设置，包括被积函数、积分区间和Romberg表的大小。合理选择这些参数可以在保证计算精度的同时，提高计算效率。
4.3 蒙特卡罗模拟脚本总结
- genMCsRJ.m脚本 ：该脚本用于生成随机抖动数据并分析其对距离测量的影响。在运行脚本前，需要根据实际需求调整参数，如积分步数、函数评估次数、重复次数等。同时，要注意随机抖动参数的设置，不同的均值和标准差会对结果产生不同的影响。
- rstMCsRJ.m脚本 ：主要用于绘制相位误差图和分析结果的不确定性。在运行该脚本前，要确保之前生成的数据文件存在，并且文件命名符合脚本中的要求。
5. 实际应用案例分析
假设我们需要进行一个距离测量实验，为了控制实验结果的不确定性，我们可以使用上述的Matlab脚本进行数据采集和分析。
5.1 数据采集阶段
- 选择DSOS604A和HP Infiniium 54846B两台示波器分别采集来自两个相机和信号发生器的多组数据。按照前面介绍的步骤，对DSOS604A示波器进行界面设置，对HP Infiniium 54846B示波器运行相应的Matlab脚本进行数据采集。
- 在采集过程中，由于相机长时间运行会发热，每个相机每隔30分钟进行一次数据采集，采集间隙相机关闭冷却。
5.2 数据分析阶段
- Romberg积分算法应用 ：在某些情况下，我们可能需要对采集到的数据进行积分运算，以评估周期性抖动对距离测量的影响。此时，可以使用Romberg积分算法脚本，输入相应的被积函数、积分区间和Romberg表的大小，得到积分结果。
- 蒙特卡罗模拟应用 ：为了评估随机抖动对距离测量的影响，我们可以运行genMCsRJ.m脚本生成随机抖动数据，并分析其对相位误差的影响。然后，使用rstMCsRJ.m脚本绘制相位误差图和分析结果的不确定性，从而更好地了解随机抖动对实验结果的影响程度。
以下是整个实验过程的mermaid流程图：
graph TD;
A[开始实验] --> B[数据采集];
B --> C{选择示波器};
C -- DSOS604A --> D[示波器界面设置采集];
C -- HP Infiniium 54846B --> E[运行Matlab脚本采集];
D --> F[数据保存];
E --> F;
F --> G[数据分析];
G --> H{分析需求};
H -- 周期性抖动分析 --> I[使用Romberg积分算法];
H -- 随机抖动分析 --> J[运行genMCsRJ.m脚本];
J --> K[运行rstMCsRJ.m脚本];
I --> L[结果输出];
K --> L;
L --> M[实验结束];
6. 未来扩展与优化方向
6.1 数据采集方面
- 可以考虑增加更多的示波器或传感器，进一步丰富采集的数据类型，以提高实验结果的准确性。
- 优化采集间隔时间的设置，根据相机的实时温度监测结果自动调整采集间隔，实现更智能的采集策略。
6.2 数据分析方面
- 对Romberg积分算法进行优化，减少计算时间和内存占用。例如，可以采用并行计算的方法，加快积分计算的速度。
- 改进蒙特卡罗模拟的算法，提高随机抖动数据的生成效率和准确性。同时，可以尝试使用更复杂的数据分析方法，如机器学习算法，对实验结果进行更深入的分析和预测。
6.3 脚本通用性方面
- 将现有的脚本进行封装，开发成更通用的函数库，方便在不同的实验场景中重复使用。同时，提供详细的文档和示例代码，降低用户的使用门槛。
通过以上的总结、案例分析和未来扩展方向的探讨，我们可以更好地理解和使用这些Matlab脚本，为实际的实验研究提供有力的支持。在实际应用中，我们可以根据具体的需求对脚本进行适当的修改和优化，以满足不同的实验要求。
扫一扫
3718
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
