超级会员免费看
数学推导、统计定理与Matlab脚本应用
1. 数学推导
1.1 误差函数
误差函数在对正态分布进行积分时会遇到，正态分布是高斯函数的一种归一化形式。对于实或复自变量 $z$，误差函数 $\text{erf}(z)$ 相对于 $\eta$ 的定义为：
[
\text{erf}(z) \equiv \frac{2}{\sqrt{\pi}} \int_{0}^{z} e^{-\eta^2} d\eta
]
误差函数具有以下有用的性质：
- $\text{erf}(0) = 0$
- $\text{erf}(\infty) = 1$
- $\text{erf}(-z) = -\text{erf}(z)$
1.2 单调函数
在区间 $(a, b)$ 上的可逆函数 $y = g(x)$ 有以下几种类型：
1.
递减
：对于任意的 $x_1, x_2 \in (a, b)$，若 $x_1 < x_2$，则 $g(x_1) \geq g(x_2)$。
2.
严格递减
：对于任意的 $x_1, x_2 \in (a, b)$，若 $x_1 < x_2$，则 $g(x_1) > g(x_2)$。
3.
递增
：对于任意的 $x_1, x_2 \in (a, b)$，若 $x_1 < x_2$，则 $g(x_1) \leq g(x_2)$。
4.
严格递增
：对于任意的 $x_1, x_2 \in (a, b)$，若 $x_1 < x_2$，则 $g(x_1) < g(x_2)$。
如果函数 $y = g(x)$ 在区间 $(a, b)$ 上可逆，并且属于上述类型之一，则称该函数在给定区间上是单调的。
下面用 mermaid 流程图展示单调函数的判断流程：
graph TD;
A[开始] --> B{判断x1 < x2};
B -- 是 --> C{判断g(x1)与g(x2)关系};
C -- g(x1) ≥ g(x2) --> D[递减函数];
C -- g(x1) > g(x2) --> E[严格递减函数];
C -- g(x1) ≤ g(x2) --> F[递增函数];
C -- g(x1) < g(x2) --> G[严格递增函数];
B -- 否 --> H[结束];
2. 统计定理
2.1 定义
2.1.1 无意识统计学家定律
当连续随机变量 $X$ 的概率分布已知，但函数 $g(X)$ 的概率分布未知时，可以使用无意识统计学家定律（LOTUS）来求 $g(X)$ 的期望值。变量 $X$ 的任何函数 $g: \mathbb{R} \to \mathbb{R}$ 的期望值为：
[
E[g(X)] = \int_{-\infty}^{\infty} g(x) f_X(x) dx
]
2.1.2 随机变量的累积分布函数
连续随机变量 $X$ 的累积分布函数（CDF）是描述随机变量分布的一种方法，定义为：
[
F_X(x) = P(X \leq x), \quad \text{对于所有的 } x \in \mathbb{R}
]
其中 $P(X \leq x)$ 是 $X$ 小于或等于元素 $x$ 的概率。
2.1.3 高斯（正态）分布
正态分布是一种将随机变量 $X$ 与累积概率相关联的概率分布，定义为：
[
f_X(x) = \frac{1}{\sigma \sqrt{2\pi}} \exp\left(-\frac{1}{2} \left(\frac{x - \mu}{\sigma}\right)^2\right)
]
其中 $x \in X$，$\mu$ 和 $\sigma^2$ 分别是变量 $X$ 的均值和方差。
2.2 证明
2.2.1 LOTUS 定理
如果 $X$ 是具有概率密度函数（PDF）$f_X(x)$ 的连续随机变量，则任何函数 $g: \mathbb{R} \to \mathbb{R}$ 的均值为：
[
E[g(X)] = \int_{-\infty}^{\infty} g(x) f_X(x) dx
]
设 $Y = g(X)$ 是严格单调函数（即 $g(x)$ 可微，其反函数 $g^{-1}(y)$ 也是单调的），其导数为：
[
\frac{d}{dy} \left( g^{-1}(y) \right) = \frac{1}{g’(g^{-1}(y))}
]
通过变量替换法，可得：
[
\int_{-\infty}^{\infty} g(x) f_X(x) dx = \int_{-\infty}^{\infty} y f_X(g^{-1}(y)) \frac{1}{g’(g^{-1}(y))} dy
]
又因为 $x = g^{-1}(y)$，则 $dx = \frac{1}{g’(g^{-1}(y))} dy$。从累积分布函数可知：
[
F_Y(y) = P(Y \leq y) = P(g(X) \leq y) = P(X \leq g^{-1}(y)) = F_X(g^{-1}(y))
]
根据链式法则，$Y$ 的 PDF 为：
[
f_Y(y) = f_X(g^{-1}(y)) \frac{1}{g’(g^{-1}(y))}
]
所以：
[
\int_{-\infty}^{\infty} g(x) f_X(x) dx = \int_{-\infty}^{\infty} y f_Y(y) dy = E[Y]
]
即：
[
\int_{-\infty}^{\infty} g(x) f_X(x) dx = E[g(X)]
]
2.2.2 高斯积分
高斯积分（也称为概率积分）与误差函数 $\text{erf}(.)$ 相关，是一维高斯函数在 $(-\infty, \infty)$ 上的积分：
[
I = \int_{-\infty}^{\infty} e^{-x^2} dx
]
可以通过组合两个一维高斯变量来计算。则：
[
I^2 = \int_{-\infty}^{\infty} e^{-x^2} dx \int_{-\infty}^{\infty} e^{-y^2} dy = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} e^{-(x^2 + y^2)} dx dy
]
转换为极坐标，$x^2 + y^2 = r^2$，$dx dy = r dr d\theta$，可得：
[
I = \sqrt{\int_{0}^{2\pi} \int_{0}^{\infty} e^{-r^2} r dr d\theta} = \sqrt{2\pi \left[ \frac{e^{-r^2}}{-2} \right]_{0}^{\infty}} = \sqrt{\pi}
]
这个积分可用于求正态分布在给定区间（如 $\mu \pm \sigma$，$\mu \pm 2\sigma$，$\mu \pm 3\sigma$）内的相应面积，结果分别约为分布的 68.27%、95.45% 和 99.73%。
2.2.3 一些参数的均值和方差
2.2.3.1 连续随机变量
设 $X$ 是连续随机变量，$X$ 的期望值定义为：
[
E[X] = \int_{-\infty}^{\infty} x f_X(x) dx
]
方差为：
[
V[X] = E[(X - E[X])^2] = \int_{-\infty}^{\infty} (X - E[X])^2 f_X(x) dx = E[X^2] - (E[X])^2
]
其中 $E[X^2]$ 是变量 $X^2$ 的期望值，根据 LOTUS 定理：
[
E[X^2] = \int_{-\infty}^{\infty} x^2 f_X(x) dx
]
2.2.3.2 随机变量的线性组合
设 $Y = \sum_{k = 1}^{K} \alpha_k X_k$ 是 $K$ 个独立随机变量 $X_k$ 的线性组合，其中 $\alpha_k$ 为实常数，$\mu_k$ 和 $\sigma_k^2$ 分别是随机变量 $X_k$ 的均值和方差。则 $Y$ 的均值 $\mu_Y$ 和方差 $\sigma_Y^2$ 分别为：
[
\mu_Y = E[Y] = \sum_{k = 1}^{K} \alpha_k E[X_k] = \sum_{k = 1}^{K} \alpha_k \mu_k
]
[
\sigma_Y^2 = V[Y] = \sum_{k = 1}^{K} \alpha_k^2 \sigma_k^2
]
2.2.3.3 高斯分布
根据均值的定义，将高斯分布的概率密度函数代入可得：
[
E[X] = \int_{-\infty}^{\infty} x \frac{1}{\sigma \sqrt{2\pi}} \exp\left(-\frac{1}{2} \left(\frac{x - \mu}{\sigma}\right)^2\right) dx = \mu
]
方差为：
[
V[X] = \int_{-\infty}^{\infty} x^2 \frac{1}{\sigma \sqrt{2\pi}} \exp\left(-\frac{1}{2} \left(\frac{x - \mu}{\sigma}\right)^2\right) dx - \mu^2 = \sigma^2
]
下面用表格总结一些统计定理中的关键公式：
| 定理或参数 | 公式 |
| — | — |
| LOTUS 定理 | $E[g(X)] = \int_{-\infty}^{\infty} g(x) f_X(x) dx$ |
| 累积分布函数 | $F_X(x) = P(X \leq x)$ |
| 高斯分布概率密度函数 | $f_X(x) = \frac{1}{\sigma \sqrt{2\pi}} \exp\left(-\frac{1}{2} \left(\frac{x - \mu}{\sigma}\right)^2\right)$ |
| 连续随机变量均值 | $E[X] = \int_{-\infty}^{\infty} x f_X(x) dx$ |
| 连续随机变量方差 | $V[X] = E[X^2] - (E[X])^2$ |
| 随机变量线性组合均值 | $\mu_Y = \sum_{k = 1}^{K} \alpha_k \mu_k$ |
| 随机变量线性组合方差 | $\sigma_Y^2 = \sum_{k = 1}^{K} \alpha_k^2 \sigma_k^2$ |
| 高斯分布均值 | $E[X] = \mu$ |
| 高斯分布方差 | $V[X] = \sigma^2$ |
3. 特定的 Matlab 脚本
3.1 模拟数据的抖动测量
这部分包括两个文件：一个用于生成模拟数据，另一个用于提取模拟信号中的抖动。
3.1.1 生成模拟数据：genSG10to50.m
% % Generating simulation data for check the propsed algorithm:
% % File name : genSG10to50.m --> includes for RJ, PJ and RJ + PJ;
% % Prepared by Gehan Anthonys -on September, 2017;
% % Parameters as for oscilloscope:
sampRate = 8e9; % 8 GSa/s
Ts = 1/sampRate; % 125 ps
tt = (0:1:2^14 - 1) * Ts; % 65536 = 2^16 data points
lenTT = length(tt);
amplSig = 1;
modFreq = (10:10:50) * 1e+6;
% % 1 --> RJ only, 2 --> PJ only, 3 --> both.
boolPjRjBoth = 3;
switch boolPjRjBoth
case 1
foldName = 'withRJ\20 sets\';
case 2
foldName = 'withPJ\20 sets\';
case 3
foldName = 'withRJ+PJ\20 sets\';
end
% for RJ
if boolPjRjBoth == 1 || boolPjRjBoth == 3
lenPts = lenTT / 16 + 1; % no of random numbers : 4097 for ../16
rmsRJ = 5e-12; % rms: 5, 10, 20, 30, 40, 50 ps
valRJ = rmsRJ * randn(1, lenPts); % with 5 ps (rms) RJ
for k = 1:lenPts - 1
coeff = polyfit([tt(16*(k - 1) + 1), tt(16*k)], [valRJ(k), valRJ(k + 1)], 1);
interPts{k} = coeff(1) * tt(16*(k - 1) + 1:16*k) + coeff(2);
end
horzIntRJ = interPts{1};
for k = 2:length(interPts)
horzIntRJ = horzcat(horzIntRJ, interPts{k});
end
intpolRJ = horzIntRJ; % rms 5 ps of RJ
% % To check the RJ:
figure;
plot(tt, intpolRJ, '.-');
grid on;
% XOR --
consRJ = repmat(rmsRJ, 1, lenTT); % const rms 5ps of RJ, but the jitter
% will be cancelled out even the 1st element set to 0,
consRJ(1) = 0; % since ideal and data signals aligned with 1st element,
% THEN, nearly ALL POINTS ARE EXACTLY ALINGED. Therfore, this is USELESS.
% XOR --
halfConsRJ = horzcat(zeros(1, lenTT / 2), repmat(rmsRJ, 1, lenTT / 2));
% half zero + half const rms 5ps of RJ
end
% for PJs
if boolPjRjBoth == 2 || boolPjRjBoth == 3
amplPJ1 = 5e-12; % ampl: 5, 10, 20, 30, 40, 50 ps
freqPJ1 = 4.37e+6;
amplPJ2 = 0;
freqPJ2 = 10e+6;
end
nSet = 20; % no of data sets
for nn = 1:nSet
for ii = 3:3
switch boolPjRjBoth
case 1 % with RJ only
vol = amplSig * sin(2 * pi * modFreq(ii) * (tt + intpolRJ));
% with rms value of RJ, XOR --
% vol = amplSig * sin(2 * pi * modFreq(ii) * (tt + consRJ));
% with const rms 5ps of RJ, XOR --
% vol = amplSig * sin(2 * pi * modFreq(ii) * (tt + halfConsRJ));
% with half const rms 5 ps of RJ
case 2 % with PJs only
vol = amplSig * sin(2 * pi * modFreq(ii) * (tt + amplPJ1 * cos(2 * pi * freqPJ1 * tt) + amplPJ2 * cos(2 * pi * freqPJ2 * tt)));
case 3 % with (PJ + RJ)
vol = amplSig * sin(2 * pi * modFreq(ii) * (tt + amplPJ1 * cos(2 * pi * freqPJ1 * tt) + intpolRJ));
otherwise % WITHOUT jitter
vol = amplSig * sin(2 * pi * modFreq(ii) * tt);
end
% % To check the RJ:
figure;
plot(tt, vol, '.-');
grid on;
data = [tt', vol'];
pathDir = strcat('E:\GA\Matlab\JitterAnalysis\wavesSR\jitsCalculation\SR4000\sSG8\', foldName); % sSG6, sSG7, sSG8
cd(pathDir);
filename = strcat('sSG', num2str(10 * ii), '_8GSa', num2str(nn), '.csv');
csvwrite(filename, data, 0, 0);
end
end
pathDir = 'E:\GA\Matlab\JitterAnalysis\wavesSR\jitsCalculation\SR4000\';
cd(pathDir);
% % End of the file: genSG10to50.m --------
3.1.2 提取抖动：rstSG10to50.m
% % Extracting the jitter which the simulated data by genSG10to50.m
% % File name : rstSG10to50.m --> includes for RJ, PJ and RJ + PJ;
% % Prepared by Gehan Anthonys -on September, 2017;
% %
% % DATA SAVED in .xlsx are NOT in t−domain, they are in f−domain abs(jitter) values, i.e. mean jitter of FFT spectra
clear all;
folTOF = 8;
dataDir = strcat('E:\GA\Matlab\JitterAnalysis\wavesSR\jitsCalculation\SR4000\sSG', num2str(folTOF), '\');
cd(dataDir);
splitDataDir = strsplit(dataDir, '\');
simData = 1; % set 1 for simulation (generated) data and 0 for SigGen data
nSet = 20;
sigFrq = (10:10:50) * 1e+6;
ifSmth = 0; %% 1 for smoothing, 0 for NOT --> specially for simulated data only
ifWndw = 0; %% 1 for windowing, 0 for NOT, before fft of jitter
if ifWndw == 1
strWndw = ', with windowing';
crrctFctrEny = sqrt(8 / 3); % for Hann, energy
crrctFctrAmp = 2; % for Hann, amplitude
else
strWndw = ', without windowing';
crrctFctrEny = 1;
crrctFctrAmp = 1; % for Hann, amplitude
end
for i = 1:5
if simData == 1 % simulation (generated) data
sheet = strcat('sSG', num2str(i * 10));
if ifSmth == 1
if ifWndw == 1
dataSG = abs(xlsread('rstWndwSimSG10to50.xlsx', sheet, 'B:B'));
% in order to generalise take the absolute
else
dataSG = abs(xlsread('rstNoWndwSimSG10to50.xlsx', sheet, 'B:B'));
end
else
if ifWndw == 1
dataSG = abs(xlsread('rstWndwSimSG10to50wo.xlsx', sheet, 'B:B'));
else
dataSG = abs(xlsread('rstNoWndwSimSG10to50wo.xlsx', sheet, 'B:B'));
end
end
else
sheet = strcat('SG', num2str(i * 10)); % Sig Gen data
if ifWndw == 1
dataSG = abs(xlsread('rstWndwSG10to50.xlsx', sheet, 'B:B')); % Sig Gen data: fft
else
dataSG = abs(xlsread('rstNoWndwSG10to50.xlsx', sheet, 'B:B')); % Sig Gen data: fft
end
end
valFfSG{i} = dataSG;
lenFfSG(i) = length(dataSG);
meanFfSG(i) = mean(dataSG);
stdFfSG(i) = std(dataSG);
rmsFfSG(i) = rms(dataSG);
trshldLvlDataSG = meanFfSG(i);
jFrq = sigFrq(i) * (0:lenFfSG(i) - 1) / lenFfSG(i);
jFrqHlf = sigFrq(i) * (0:lenFfSG(i) / 2 - 1) / lenFfSG(i);
maxFfSG(i) = max(dataSG);
% % To find the area under the curve (total, RJ, PJ) by using fft
% To CHECK: whether have the PJ or not, in the signals
figure;
plot(jFrqHlf / 1e+6, dataSG(1:end / 2) * 1e+12, 'o-', jFrqHlf / 1e+6, repmat(trshldLvlDataSG, 1, lenFfSG(i) / 2) * 1e+12, '.-');
grid on;
sheet = 'SET.A'; %'SET.A -for RJ, 'SET.B -for PJ, 'SET.C -for RJ+PJ
if strcmp(upper(sheet), 'SET.A') == 1
thrshld = 3;
else
thrshld = 1;
end
indMaxFfSG(i) = find(dataSG == maxFfSG(i), 1, 'first'); % peak index,
% (i.e. @ PJ)
if indMaxFfSG(i) < length(jFrqHlf) && fix(maxFfSG(i) ./ meanFfSG(i)) > thrshld
% if more than 3, ASSUMED exist PJ. 'fix: Round toward zero'
existPJ = 1;
freqMaxSG(i) = jFrqHlf(indMaxFfSG(i));
indPrevMax(i) = indMaxFfSG(i) - 1;
indNextMax(i) = indMaxFfSG(i) + 1;
indsLvlPrevMax = find(dataSG(1:indMaxFfSG(i)) > trshldLvlDataSG);
indsLvlNextMax = find(dataSG(indMaxFfSG(i):end / 2) > trshldLvlDataSG);
% to find prev index wrt max
if isempty(find(diff(indsLvlPrevMax) == 1, 1))
indPrevMax(i) = indsLvlPrevMax(end) - 1;
else
for aa = length(indsLvlPrevMax):-1:1
if dataSG(indsLvlPrevMax(aa) - 1) > trshldLvlDataSG
indPrevMax(i) = indPrevMax(i) - 1;
else
break;
end
end
end
% to find next index wrt max
if isempty(find(diff(indsLvlNextMax) == 1, 1, 'last'))
indNextMax(i) = indMaxFfSG(i) + 1;
else
for zz = 1:length(indsLvlNextMax)
if dataSG(indMaxFfSG(i) + indsLvlNextMax(zz)) > trshldLvlDataSG
if dataSG(indMaxFfSG(i) + indsLvlNextMax(zz)) < dataSG(indMaxFfSG(i) + indsLvlNextMax(zz) - 1)
indNextMax(i) = indNextMax(i) + 1;
else
indNextMax(i) = indNextMax(i) - 1;
break;
end
else
break;
end
end
end
else
existPJ = 0;
freqMaxSG(i) = 0;
indPrevMax(i) = 1;
indNextMax(i) = floor(lenFfSG(i) / 2);
dataLvlPJ_RJ = 0;
end
% % To find the PSD of partial RJ,
% % we need the intermediate points inside the PJ bar -
% % ---by using linear interpolation
intpolPartRJ{i} = interp1([jFrqHlf(indPrevMax(i)) jFrqHlf(indNextMax(i))], ...
[dataSG(indPrevMax(i)) dataSG(indNextMax(i))], ...
jFrqHlf(indPrevMax(i):indNextMax(i)));
% % Calculating areas & power of the signal
% % Half: total area under the curve (ie. PJ + RJ)
areaFfSG(i) = 1 * trapz(dataSG(1:end / 2));
% % Total power spectrum
psFfSG(i) = sum(dataSG(1:end / 2) .* dataSG(1:end / 2)); % by half
psFfSG_A(i) = sum(dataSG(1:end / 2)); % by amplitude, by half
psFfSG_byFull(i) = 1 * sum(dataSG .* dataSG); % by full,
% BOTH ARE SAME when twice of half as, sqrt(2 * ...)
if existPJ == 1
figure;
plot(jFrqHlf / 1e+6, dataSG(1:end / 2) * 1e+12, 'o-', 'color', 'blue');
grid on;
if indNextMax(i) <= length(jFrqHlf)
hold on;
plot(jFrqHlf(indPrevMax(i):indNextMax(i)) / 1e6, intpolPartRJ{i} * 1e12, '+ -');
xlabel('Frequecy (MHz)');
ylabel('Amplitude (ps)');
legend(strcat(num2str(sigFrq(i) / 1e6), ' MHz'), 'interpolated level');
end
hold on;
fill(jFrqHlf(indPrevMax(i):indNextMax(i)) / 1e6, dataSG(indPrevMax(i):indNextMax(i)) * 1e12, 'b');
% OR, direct calculation from plot (half side)
lenPJpart = length(dataSG(indPrevMax(i):indNextMax(i)));
psPJFfSG(i) = sum(dataSG(indPrevMax(i):indNextMax(i)) .* dataSG(indPrevMax(i):indNextMax(i)) - ...
intpolPartRJ{i}' .* intpolPartRJ{i}');
psPJFfSG_A(i) = sum(dataSG(indPrevMax(i):indNextMax(i)) - intpolPartRJ{i}'); % by amplitude
else
figure;
plot(jFrqHlf / 1e+6, dataSG(1:end / 2) * 1e+12, 'o-', 'color', 'blue');
grid on;
xlabel('Frequecy (MHz)');
ylabel('Amplitude (ps)');
legend(strcat(num2str(sigFrq(i) / 1e6), ' MHz'));
lenPJpart = 0;
psPJFfSG(i) = 0;
psPJFfSG_A(i) = 0;
end
psdRJFfSG(i) = psFfSG(i) - psPJFfSG(i);
psdRJFfSG_A(i) = psFfSG_A(i) - psPJFfSG_A(i);
end
meanFfSG;
stdFfSG;
% % Total jitter, from spectrum
rmsFfSgnl = rmsFfSG .* sqrt(lenFfSG);
% from spectrum, to check the accuracy of total jitter,
ampFfSgnl = sqrt(2) * rmsFfSgnl;
% This should be EQUAL to 'calAmpFfSG' given by 'sqrt(signal's POWER)'
% % with correction factor
rmsFfSgnl_cF = rmsFfSG .* sqrt(crrctFctrEny) .* sqrt(lenFfSG);
% corrected on Nov 17, 2018 --> factor should in sqrt()
ampFfSgnl_cF = sqrt(2) * rmsFfSgnl_cF;
%%
% % Calculating jitter parameters:
% NOTE: FROM THE SPECTRUM IT GIVES "RMS value of each freq of time signal",
% therefore, to get the amplitude we need to MULTIPLE by sqrt(2) for PJ and
% for Total J. But for RJ it can be get as it is, since RJ measured in RMS
% %----on Jan 31, 2018.
% Total jitter
if strcmp(upper(sheet), 'SET.A') == 1
%'SET.A -for RJ, 'SET.B -for PJ, 'SET.C -for RJ+PJ
calRmsFfSG = sqrt(2 * psFfSG); % since two sided
calRmsFfSG_cF = sqrt(2 * psFfSG) * sqrt(crrctFctrEny);
% corrected on Nov 17, 2018 --> factor should in sqrt()
else % 'SET.B' or 'SET.C'
calAmpFfSG = sqrt(2) * sqrt(2 * psFfSG); % since two sided
calAmpFfSG_cF = sqrt(2) * sqrt(2 * psFfSG) * sqrt(crrctFctrEny);
% corrected on Nov 17, 2018 --> factor should in sqrt()
end
% % for PJ --> [ampl (or rms), freq], by (in f and PSD - domains)
freqPJSG = freqMaxSG;
% by PSD is given amplitude, but spectrum is given the rms amounts
calAmplPJSG = sqrt(2) * sqrt(2 * psPJFfSG); % since two sided
calAmplPJSG_cF = sqrt(2) * sqrt(2 * psPJFfSG) * sqrt(crrctFctrEny);
% corrected on Nov 17, 2018 --> factor should in sqrt()
% % for RJ --> [rms], by (in f and PSD - domains)
if strcmp(upper(sheet), 'SET.B') == 1 % if PJ only --> SET.B
calRmsRJSG = 0;
calRmsRJSG_cF = 0;
elseif strcmp(upper(sheet), 'SET.A') == 1 % if RJ only --> SET.B
calRmsRJSG = sqrt(2) * sqrt(2 * psdRJFfSG); % since two sided
calRmsRJSG_cF = sqrt(2) * sqrt(2 * psdRJFfSG) * crrctFctrEny;
else % if PJ + RJ --> SET.C
calRmsRJSG = sqrt(2 * psdRJFfSG); % since two sided
calRmsRJSG_cF = sqrt(2 * psdRJFfSG) * crrctFctrEny;
% corrected on Nov 17, 2018 --> factor should in sqrt()
end
%%
% % for plotting
for i = 1:5
padFfSG{i} = num2cell(padarray(valFfSG{i}(1:end / 2), (max(lenFfSG) - lenFfSG(i)) / 2, 0, 'post'));
end
dataSG10 = cell2mat(padFfSG{1});
dataSG20 = cell2mat(padFfSG{2});
dataSG30 = cell2mat(padFfSG{3});
dataSG40 = cell2mat(padFfSG{4});
dataSG50 = cell2mat(padFfSG{5});
jFs = 30e+6;
jffN = max(lenFfSG);
jFreq = jFs * (0:jffN - 1) / jffN;
jFreqHalf = jFs * (0:jffN / 2 - 1) / jffN;
if simData == 1 % simulation (generated) data
deviceSG = 'Simulation data: ';
measureDev = 'using computer';
else % for Sig Gen data
deviceSG = 'Sig/Gen Agilent 8648B: ';
measureDev = 'by using oscilloscope: Hp infiniium 2.25 GHz, 8 GSa/s';
end
textStr = num2str(calRmsRJSG);
% % FFT spectrum for all together, THIS is better
figure;
plot(jFreqHalf / 1e+6, dataSG10 * 1e+12, 'o-', 'color', 'blue');
hold on;
%text(1, 13, num2str(calRmsRJSG(1)), 'FontSize', 16);
plot(jFreqHalf / 1e+6, dataSG20 * 1e+12, '+ -', 'color', 'magenta');
hold on;
%text(6, 12, num2str(calRmsRJSG(2)), 'FontSize', 16);
plot(jFreqHalf / 1e+6, dataSG30 * 1e+12, 's-', 'color', 'green');
hold on;
%text(11, 11, num2str(calRmsRJSG(3)), 'FontSize', 16);
plot(jFreqHalf / 1e+6, dataSG40 * 1e+12, 'p-', 'color', 'red');
hold on;
%text(16, 10, num2str(calRmsRJSG(4)), 'FontSize', 16);
plot(jFreqHalf / 1e+6, dataSG50 * 1e+12, 'x-', 'color', 'cyan');
grid on;
%text(21, 9, num2str(calRmsRJSG(5)), 'FontSize', 16);%ylim([0 5] * 1)
title({[deviceSG, 'FFT spectrum (1st half) of the total jitter (mean) of data sets: ', num2str(nSet), ' @ ', cell2mat(splitDataDir(end - 1))], ...
[measureDev, ' for 10 MHz to 50 MHz', strWndw]}, 'FontSize', 12);
xlabel('Frequency (MHz)', 'FontSize', 24, 'FontWeight', 'bold');
ylabel('Amplitude (ps)', 'FontSize', 24, 'FontWeight', 'bold');
legend({'10 MHz', '20 MHz', '30 MHz', '40 MHz', '50 MHz'}, 'FontSize', 20);
set(gca, 'FontSize', 20);
%%
% % save the data
if simData == 1 % simulation (generated) data
if ifWndw == 1
xlsName = 'rstWndwSimAEF.xlsx';
% with windowing before fft of jitter
else
xlsName = 'rstNoWndwSimAEF.xlsx';
end
[Num, Txt, Raw] = xlsread(xlsName, sheet);
nRows = size(Raw, 1);
dataAtOut = strcat('A', num2str(nRows + 1));
dataCal = [sigFrq; freqPJSG; calAmplPJSG; calRmsRJSG];
xlswrite(xlsName, dataCal', sheet, dataAtOut);
end
dataDir = 'E:\GA\Matlab\JitterAnalysis\wavesSR\jitsCalculation\SR4000\';
cd(dataDir);
% % End of the file: rstSG10to50.m --------
3.2 操作步骤
3.2.1 生成模拟数据
- 设置参数，如采样率、调制频率等。
-
根据
boolPjRjBoth的值选择生成的数据类型（RJ 单独、PJ 单独或两者都有）。 - 生成随机抖动（RJ）和周期性抖动（PJ）的数据。
- 循环生成多个数据集，并将数据保存为 CSV 文件。
3.2.2 提取抖动
- 清空工作区，设置数据目录和相关参数。
-
根据
simData的值选择读取模拟数据还是 SigGen 数据。 -
根据
ifSmth和ifWndw的值选择是否进行平滑和加窗处理。 - 读取数据，计算相关统计量，如均值、标准差、均方根等。
- 判断是否存在周期性抖动（PJ），并计算相关参数。
- 计算抖动的总面积、功率谱等。
- 绘制 FFT 频谱图，进行可视化分析。
-
根据
simData和ifWndw的值选择保存数据的文件，并将计算结果保存到 Excel 文件中。
通过以上数学推导、统计定理的证明以及 Matlab 脚本的实现，可以更好地理解和分析相关的数学和统计问题，同时利用 Matlab 工具进行数据的生成、处理和可视化。
4. 模拟数据生成与抖动提取流程总结
4.1 数据生成与抖动提取流程图
graph LR;
A[开始] --> B[设置参数]
B --> C{选择数据类型}
C -- RJ单独 --> D[生成RJ数据]
C -- PJ单独 --> E[生成PJ数据]
C -- RJ+PJ --> F[生成RJ和PJ数据]
D --> G[生成数据集]
E --> G
F --> G
G --> H[保存为CSV文件]
H --> I[结束数据生成]
I --> J[清空工作区并设置参数]
J --> K{选择数据来源}
K -- 模拟数据 --> L[读取模拟数据文件]
K -- SigGen数据 --> M[读取SigGen数据文件]
L --> N{是否平滑和加窗}
M --> N
N -- 是 --> O[平滑和加窗处理]
N -- 否 --> P[直接处理数据]
O --> Q[计算统计量]
P --> Q
Q --> R{是否存在PJ}
R -- 是 --> S[计算PJ相关参数]
R -- 否 --> T[跳过PJ计算]
S --> U[计算抖动参数]
T --> U
U --> V[绘制FFT频谱图]
V --> W{保存数据}
W -- 是 --> X[保存到Excel文件]
W -- 否 --> Y[结束]
X --> Y
4.2 关键参数说明
| 参数 | 说明 |
|---|---|
sampRate
| 采样率，单位为 GSa/s |
Ts
| 采样周期，单位为 ps |
modFreq
| 调制频率，单位为 MHz |
boolPjRjBoth
| 选择生成的数据类型，1 表示 RJ 单独，2 表示 PJ 单独，3 表示 RJ + PJ |
rmsRJ
| 随机抖动（RJ）的均方根值，单位为 ps |
amplPJ1
、
amplPJ2
| 周期性抖动（PJ）的振幅，单位为 ps |
freqPJ1
、
freqPJ2
| 周期性抖动（PJ）的频率，单位为 Hz |
simData
| 选择数据来源，1 表示模拟数据，0 表示 SigGen 数据 |
ifSmth
| 是否进行平滑处理，1 表示是，0 表示否 |
ifWndw
| 是否进行加窗处理，1 表示是，0 表示否 |
5. 抖动分析的实际应用与意义
5.1 实际应用场景
抖动分析在许多领域都有重要的应用，例如：
-
通信系统
：在高速数据传输中，抖动会影响信号的准确性和可靠性，通过分析抖动可以优化通信系统的性能，提高数据传输的质量。
-
电子设备
：在时钟信号、电源信号等方面，抖动可能导致设备的不稳定运行，通过对抖动的监测和分析，可以及时发现潜在的问题，保证设备的正常工作。
-
测试测量
：在使用示波器等仪器进行信号测量时，抖动分析可以帮助准确评估信号的质量，为测试结果提供更可靠的依据。
5.2 抖动分析的意义
- 性能评估 ：通过对抖动的分析，可以量化信号的稳定性和准确性，评估系统或设备的性能。
- 故障诊断 ：抖动的异常变化可能是系统故障或潜在问题的表现，通过监测抖动可以及时发现故障并进行诊断。
- 优化设计 ：了解抖动的来源和特性，可以为系统或设备的设计提供参考，优化设计方案，降低抖动的影响。
6. 总结与展望
6.1 总结
本文介绍了数学推导中的误差函数和单调函数的相关知识，以及统计定理中的无意识统计学家定律、累积分布函数、高斯分布等内容，并对这些定理进行了证明。同时，详细阐述了模拟数据的生成和抖动提取的 Matlab 脚本，包括具体的代码实现和操作步骤。通过流程图和表格对数据生成与抖动提取的流程进行了总结，明确了关键参数的含义。最后，探讨了抖动分析在实际应用中的场景和意义。
6.2 展望
随着科技的不断发展，抖动分析在更多领域的应用将越来越广泛。未来可以进一步研究更复杂的抖动模型，提高抖动分析的精度和效率。同时，可以结合机器学习等技术，实现对抖动的自动监测和预测，为系统的稳定性和可靠性提供更有力的保障。此外，还可以探索抖动分析在新兴领域（如量子通信、人工智能硬件等）的应用，为这些领域的发展提供支持。
通过对抖动分析的深入研究和应用，可以更好地应对各种复杂的信号处理问题，推动相关领域的技术进步。
扫一扫
6
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
