bowei
ea4541ee90
- simulink_rtl_dev / rtl-ip-replace: 从工程级 .claude/skills 迁到用户级并去工程耦合 (模型名/Notes 文档依赖/工程实测基线移出; addpath 与脚本路径改用户级) - rtl-ip-replace/scripts/gen_cmp_tb.py: Vivado 路径不再写死, 改读 VIVADO_PATH / XILINX_VIVADO (兼容 install root / bin / bin\vivado.bat 三种写法), 回落时提示 去设 skills_env.md 的 VIVADO_PATH - simulink-rtl-redev: 新增。原生模型改动后的仿真模型同步规范, 含两个实测过的工具 check_sync.py (三条时序链新鲜度闸门) 与 measure_k.m (时域/频域/包络三估计器互证) - sim-rtl-flow/SKILL.md: 加入 sync_sim:* 步骤, 引用上述三个 skill
200 lines
8.7 KiB
Matlab
200 lines
8.7 KiB
Matlab
function [K, R] = measure_k(ref, dut, varargin)
|
|
%MEASURE_K Measure the integer delay K that <dut> lags <ref>, by TWO INDEPENDENT
|
|
% methods that must agree (Task16-3 section 5: "延迟必须实测,不能按设计意图推算.
|
|
% 用独立微台、至少两种互不依赖的方法测出一致值").
|
|
%
|
|
% K = measure_k(ref, dut)
|
|
% [K, R] = measure_k(ref, dut, 'Name', value, ...)
|
|
%
|
|
% ref, dut : real or COMPLEX vectors, same sample rate. dut is the signal that
|
|
% comes out K samples later than ref.
|
|
%
|
|
% Options
|
|
% 'MaxLag' (default 4096) search window, in samples
|
|
% 'OnsetTol' (default 1e-3) onset threshold, RELATIVE to each signal's own
|
|
% RMS -- so the two signals may differ in gain
|
|
% 'FracTol' (default 0.05) how far the interpolated correlation peak may sit
|
|
% from an integer before K is declared non-integer
|
|
% 'Verbose' (default true)
|
|
%
|
|
% Returns
|
|
% K : the agreed integer delay, or NaN if the methods disagree
|
|
% R : struct with per-method results and the verdict --
|
|
% .k_corr argmax of the normalised cross-correlation
|
|
% .k_onset difference of first-significant-sample indices
|
|
% .k_frac parabolic-interpolated peak (fractional)
|
|
% .frac_resid |k_frac - round(k_frac)| -- large => not a pure integer shift
|
|
% .rho |correlation| at the peak, 0..1
|
|
% .agree k_corr == k_onset
|
|
% .pass agree && frac_resid <= FracTol
|
|
%
|
|
% WHY TWO METHODS: a correlation argmax alone will happily return a number for
|
|
% two unrelated signals. Onset detection is a structurally different estimator
|
|
% (amplitude threshold, not inner product), so agreement between them is real
|
|
% evidence. The fractional peak is the third check: it catches the case where
|
|
% the change did NOT introduce a clean integer pipeline delay (e.g. it altered a
|
|
% fractional interpolator), which an integer-only estimator would silently round.
|
|
%
|
|
% Base MATLAB only -- no Signal Processing Toolbox (no xcorr/finddelay).
|
|
|
|
p = inputParser;
|
|
p.addParameter('MaxLag', 4096, @(x) isnumeric(x) && isscalar(x) && x > 0);
|
|
p.addParameter('OnsetTol', 1e-3, @(x) isnumeric(x) && isscalar(x) && x > 0);
|
|
p.addParameter('FracTol', 0.05, @(x) isnumeric(x) && isscalar(x) && x > 0);
|
|
p.addParameter('Verbose', true, @(x) islogical(x) || isnumeric(x));
|
|
p.parse(varargin{:});
|
|
o = p.Results;
|
|
|
|
ref = ref(:); dut = dut(:);
|
|
n = min(numel(ref), numel(dut));
|
|
if n < 16
|
|
error('measure_k: need at least 16 samples, got %d', n);
|
|
end
|
|
ref = ref(1:n); dut = dut(1:n);
|
|
|
|
maxlag = min(o.MaxLag, n - 1);
|
|
|
|
% ---- method 1: normalised cross-correlation, argmax over 0..maxlag -----------
|
|
% FFT-based so it stays base MATLAB and is fast for long vectors.
|
|
nfft = 2^nextpow2(2*n);
|
|
Rxy = ifft(fft(dut, nfft) .* conj(fft(ref, nfft))); % dut delayed by k => peak at k
|
|
c = abs(Rxy(1:maxlag+1));
|
|
den = norm(ref) * norm(dut);
|
|
if den == 0
|
|
error('measure_k: one of the signals is all zeros');
|
|
end
|
|
c = c / den;
|
|
[rho, idx] = max(c);
|
|
k_corr = idx - 1;
|
|
|
|
% ---- method 3 (on the same curve): parabolic interpolation of the peak -------
|
|
% Done here because it needs c; it is reported separately and is NOT one of the
|
|
% "two independent methods" -- it only tests integer-ness.
|
|
if idx > 1 && idx < numel(c)
|
|
y0 = c(idx-1); y1 = c(idx); y2 = c(idx+1);
|
|
d = (y0 - y2) / (2*(y0 - 2*y1 + y2));
|
|
if ~isfinite(d), d = 0; end
|
|
k_frac = k_corr + d;
|
|
else
|
|
k_frac = k_corr;
|
|
end
|
|
frac_resid = abs(k_frac - round(k_frac));
|
|
|
|
% ---- method 2: cross-spectral phase slope (FREQUENCY domain, always available)
|
|
% A pure delay of K samples makes angle(Y.*conj(X)) a straight line of slope -K
|
|
% against normalised frequency. Fitting that slope uses no argmax at all, so it
|
|
% is a genuinely different estimator from method 1 -- and it is fractional, so it
|
|
% doubles as an integer-ness check that cannot be fooled by rounding.
|
|
k_phase = local_phase_slope(ref, dut, nfft);
|
|
|
|
% ---- method 3: onset of the first significant sample (only if one exists) ----
|
|
% Threshold is a fraction of each signal's own PEAK envelope, not its RMS: a
|
|
% noise/quantisation floor in the DUT during the reference's silent lead-in must
|
|
% not trigger it. Returns NaN for continuously excited signals -- that is not a
|
|
% failure, it just means this method does not apply.
|
|
i_ref = local_onset(ref, o.OnsetTol);
|
|
i_dut = local_onset(dut, o.OnsetTol);
|
|
if isnan(i_ref) || isnan(i_dut)
|
|
k_onset = NaN;
|
|
else
|
|
k_onset = i_dut - i_ref;
|
|
end
|
|
|
|
% ---- adjudication -----------------------------------------------------------
|
|
% Require methods 1 and 2 (always available, different domains) to agree, and
|
|
% method 3 to agree as well WHEN it applies.
|
|
agree_phase = isfinite(k_phase) && abs(k_phase - k_corr) <= o.FracTol;
|
|
agree_onset = isnan(k_onset) || (k_onset == k_corr); % NaN = not applicable
|
|
n_methods = 2 + double(~isnan(k_onset));
|
|
agree = agree_phase && agree_onset;
|
|
pass = agree && (frac_resid <= o.FracTol);
|
|
if pass, K = k_corr; else, K = NaN; end
|
|
|
|
R = struct('k_corr', k_corr, 'k_phase', k_phase, 'k_onset', k_onset, ...
|
|
'k_frac', k_frac, 'frac_resid', frac_resid, 'rho', rho, ...
|
|
'agree_phase', agree_phase, 'agree_onset', agree_onset, ...
|
|
'n_methods', n_methods, 'agree', agree, 'pass', pass, ...
|
|
'n', n, 'maxlag', maxlag);
|
|
|
|
if o.Verbose
|
|
fprintf('measure_k: n=%d maxlag=%d\n', n, maxlag);
|
|
fprintf(' method 1 xcorr argmax (time domain) K = %d |rho| = %.6f\n', k_corr, rho);
|
|
if isfinite(k_phase)
|
|
fprintf(' method 2 phase slope (freq domain) K = %.3f\n', k_phase);
|
|
else
|
|
fprintf(' method 2 phase slope (freq domain) K = n/a (not enough usable bins)\n');
|
|
end
|
|
if isnan(k_onset)
|
|
fprintf(' method 3 onset diff (amplitude) K = n/a (no clear onset -- continuous signal)\n');
|
|
else
|
|
fprintf(' method 3 onset diff (amplitude) K = %d\n', k_onset);
|
|
end
|
|
fprintf(' check interpolated peak %.3f (residual %.4f, tol %.4f)\n', ...
|
|
k_frac, frac_resid, o.FracTol);
|
|
if pass
|
|
fprintf(' VERDICT: PASS -- %d independent methods agree, K = %d\n', n_methods, K);
|
|
else
|
|
fprintf(' VERDICT: FAIL -- K NOT established\n');
|
|
if ~agree_phase
|
|
fprintf(' methods 1 and 2 disagree (%d vs %.3f). Different domains giving\n', ...
|
|
k_corr, k_phase);
|
|
fprintf(' different answers means these are not a pure delayed copy of each\n');
|
|
fprintf(' other -- check you compared the right taps.\n');
|
|
end
|
|
if ~agree_onset
|
|
fprintf(' method 3 disagrees (%d vs %d). A noise/quantisation floor in the\n', ...
|
|
k_corr, k_onset);
|
|
fprintf(' DUT during the reference lead-in can do this: raise OnsetTol.\n');
|
|
end
|
|
if frac_resid > o.FracTol
|
|
fprintf(' peak sits %.3f off an integer: the change did not introduce a\n', frac_resid);
|
|
fprintf(' clean integer pipeline delay (fractional/resampling path?).\n');
|
|
end
|
|
fprintf(' Do NOT fall back to the design-intent value (Task16-3 section 5).\n');
|
|
end
|
|
end
|
|
end
|
|
|
|
|
|
function kp = local_phase_slope(ref, dut, nfft)
|
|
%LOCAL_PHASE_SLOPE Delay from the slope of the cross-spectrum phase.
|
|
% Only bins where the reference has real energy are used, and the phase is
|
|
% unwrapped over those bins before a weighted least-squares line fit.
|
|
X = fft(ref, nfft); Y = fft(dut, nfft);
|
|
half = floor(nfft/2);
|
|
X = X(2:half); Y = Y(2:half); % drop DC and Nyquist
|
|
w = abs(X);
|
|
if isempty(w) || max(w) == 0, kp = NaN; return, end
|
|
keep = w > 0.10 * max(w); % usable bins only
|
|
if nnz(keep) < 16, kp = NaN; return, end
|
|
k = (1:numel(X)).';
|
|
ph = unwrap(angle(Y(keep) .* conj(X(keep))));
|
|
kk = k(keep);
|
|
ww = w(keep);
|
|
% weighted least squares: ph ~ a + b*(2*pi*kk/nfft), delay = -b
|
|
x = 2*pi*kk/nfft;
|
|
sw = sum(ww);
|
|
mx = sum(ww.*x)/sw; my = sum(ww.*ph)/sw;
|
|
b = sum(ww.*(x-mx).*(ph-my)) / sum(ww.*(x-mx).^2);
|
|
kp = -b;
|
|
end
|
|
|
|
|
|
function i0 = local_onset(x, tol)
|
|
%LOCAL_ONSET First index where the SMOOTHED envelope exceeds tol * its own PEAK.
|
|
% Threshold is relative to the peak, not the RMS, so a noise or quantisation
|
|
% floor well below the signal does not trigger it. Returns NaN when the signal
|
|
% is already at level from the first sample (continuously excited) -- that means
|
|
% "this method does not apply", not "failure".
|
|
a = abs(x);
|
|
w = 8; % short moving average
|
|
e = filter(ones(w,1)/w, 1, a);
|
|
pk = max(e);
|
|
if pk <= 0, i0 = NaN; return, end
|
|
thr = tol * pk;
|
|
i0 = find(e > thr, 1, 'first');
|
|
if isempty(i0), i0 = NaN; return, end
|
|
% already hot at the very start => no usable onset
|
|
if i0 <= w, i0 = NaN; end
|
|
end
|