From b5d856c1574d86d4d0af19c91e7aae30b2330c77 Mon Sep 17 00:00:00 2001 From: Thanasis Date: Tue, 2 Jun 2026 08:46:38 +0200 Subject: [PATCH] Update README.md --- README.md | 626 +++++++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 625 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index b678cab..4c587e1 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,626 @@ -# druckluft_methodik +# Workflow: Calibration-Based Distance Estimate from Welch Spectral Ratios +This workflow estimates the distance of an unknown excitation from sensor 1 by comparing its relative Welch peak-band levels to calibration measurements taken at known source distances. + +The calibration table is built from a **known calibration signal** using **cross-spectral transfer-function estimates**, rather than from raw sensor Welch auto-spectral levels. + +The method should be understood as a **calibration-based spectral-ratio estimator**. It is valid only within the calibrated measurement geometry and for comparable excitation conditions. + +--- + +## 1. Use sensor 1 as the reference sensor + +Use sensor 1 as the reference for all relative level comparisons. + +This workflow assumes that sensor 1 contains a clear, stable, high-SNR peak at the excitation frequency used for matching. If sensor 1 does not show the peak clearly, then it is a poor reference sensor for this workflow, because both ratios $r_{21}$ and $r_{31}$ depend on the reliability of $A_1$. + +The goal is not to use the absolute Welch level at one sensor. Instead, the goal is to compare how strong the same spectral peak appears at sensors 2 and 3 relative to sensor 1. + +This reduces sensitivity to unknown excitation strength, but only if the reference peak in sensor 1 is reliable. + +--- + +## 2. Compute spectra with fixed settings + +Use the same spectral settings for calibration and unknown measurements: + +* sampling frequency, +* segment length, +* window type, +* overlap, +* FFT length, +* detrending method, +* scaling convention, +* averaging method. + +For the unknown measurement, compute Welch spectra for all three sensors: + +$$ +P_1(f),\quad P_2(f),\quad P_3(f). +$$ + +For the calibration measurements, compute cross spectra between each sensor and the known calibration signal, as described in Section 6. + +Let the discrete frequency bins be + +$$ +f_m,\qquad m=0,1,2,\ldots,M-1, +$$ + +with bin spacing + +$$ +\Delta f = \frac{f_s}{N_{\mathrm{FFT}}}. +$$ + +--- + +## 3. Use the known calibration frequency as the matching frequency + +Use the known calibration signal to define the matching frequency. If the known calibration signal is a tone or narrowband excitation at frequency $f_0$, choose the FFT bin closest to that known frequency: + +$$ +m_0 = \arg\min_m |f_m-f_0|. +$$ + +Use this same bin for every calibration distance and for the unknown measurement. + +Do not substitute an automatically selected sensor-response peak as a fallback. That would make the calibration feature depend on the sensor/structure response instead of the known calibration input, and it is a different method. + +The selected frequency should be: + +* known from the calibration signal, +* clearly present in the calibration reference signal $u(t)$, +* coherent with the sensor responses during calibration, +* visible in the relevant sensors, +* isolated from neighboring peaks, +* repeatable between calibration and unknown measurements. + +--- + +## 4. Use a narrow peak-band level around the selected frequency + +Even if the Welch/FFT peaks look very sharp and clear, do not use only the single bin $m_0$ unless there is a specific reason to do so. A single-bin value can change due to small frequency shifts, finite FFT resolution, leakage, windowing, or slight differences between repeated measurements. + +For this workflow, use a **three-bin peak band** centered on the selected bin: + +$$ +\mathcal{B}(m_0)=\{m_0-1,\;m_0,\;m_0+1\}. +$$ + +This is narrow enough to represent the same sharp peak, while being less sensitive to one-bin jitter. + +For the unknown measurement, compute the peak-band power as + +$$ +A_i = \sum_{m\in\mathcal{B}(m_0)} P_i(f_m)\Delta f, +\qquad i=1,2,3. +$$ + +If the Welch output is already a power spectrum rather than a power spectral density, use + +$$ +A_i = \sum_{m\in\mathcal{B}(m_0)} P_i(f_m), +\qquad i=1,2,3. +$$ + +Because the same frequency bins are used for all sensors, the factor $\Delta f$ cancels in the ratios. It is included only for dimensional correctness when $P_i(f_m)$ is a power spectral density. + +Use the same peak-band rule for all calibration and unknown measurements. + +Implementation note: when Welch spectra are computed from `np.fft.rfft` and frequencies from `np.fft.rfftfreq`, the three-bin integration is just the sum of the already averaged Welch power/PSD values over indices `m0-1:m0+2`. Use the same peak index $m_0$ for all sensors. Sum power or PSD values, not raw complex FFT coefficients. If the Welch result is a PSD, multiplying by $\Delta f$ is dimensionally correct, but this factor cancels in the ratios. + +--- + +## 5. Convert unknown peak-band powers to dB ratios + +Convert the unknown peak-band powers to dB: + +$$ +L_{i,\mathrm{dB}} = 10\log_{10} A_i, +\qquad i=1,2,3. +$$ + +Then compute the two relative level differences: + +$$ +r_{21}^{\mathrm{unknown}} = L_{2,\mathrm{dB}} - L_{1,\mathrm{dB}}, +$$ + +$$ +r_{31}^{\mathrm{unknown}} = L_{3,\mathrm{dB}} - L_{1,\mathrm{dB}}. +$$ + +Equivalently, + +$$ +r_{21}^{\mathrm{unknown}}=10\log_{10}\frac{A_2}{A_1}, +$$ + +$$ +r_{31}^{\mathrm{unknown}}=10\log_{10}\frac{A_3}{A_1}. +$$ + +These ratios compare sensors 2 and 3 to the same reference sensor and therefore reduce sensitivity to unknown excitation amplitude. + +--- + +## 6. Build the calibration table using a known calibration signal and cross spectra + +Do not build the calibration table from raw Welch output levels. Use a measured reference signal for the known calibration input, for example the shaker drive, actuator command, hammer force, injected signal, or another signal directly related to the known calibration excitation. + +Let the calibration reference be + +$$ +u(t). +$$ + +For each sensor signal $x_i(t)$ at known distance $d_k$, compute the cross spectrum between the sensor and the calibration reference: + +$$ +S_{x_i u}(f;d_k) += +\left\langle X_i(f;d_k)U^*(f)\right\rangle, +$$ + +and the auto spectrum of the calibration reference: + +$$ +S_{uu}(f) += +\left\langle U(f)U^*(f)\right\rangle. +$$ + +Here $X_i(f;d_k)$ and $U(f)$ are the windowed FFTs of matching segments of $x_i(t)$ and $u(t)$, $(\cdot)^*$ denotes complex conjugation, and $\langle \cdot \rangle$ denotes Welch averaging over segments. + +Estimate the transfer function from the calibration input to sensor $i$ as + +$$ +\hat H_i(f;d_k) += +\frac{S_{x_i u}(f;d_k)}{S_{uu}(f)}, +\qquad i=1,2,3. +$$ + +Then use the same three-bin peak band around the known calibration frequency: + +$$ +\mathcal B(m_0)=\{m_0-1,m_0,m_0+1\}. +$$ + +Compute the peak-band transfer magnitude for each sensor: + +$$ +G_i(d_k) += +\sum_{m\in\mathcal B(m_0)} +\left|\hat H_i(f_m;d_k)\right|^2, +\qquad i=1,2,3. +$$ + +The calibration ratios are then + +$$ +r_{21}^{\mathrm{cal}}(d_k) += +10\log_{10} +\frac{G_2(d_k)}{G_1(d_k)}, +$$ + +$$ +r_{31}^{\mathrm{cal}}(d_k) += +10\log_{10} +\frac{G_3(d_k)}{G_1(d_k)}. +$$ + +Equivalently, this is the same as comparing RMS transfer magnitudes in dB: + +$$ +r_{21}^{\mathrm{cal}}(d_k) += +20\log_{10} +\frac{\sqrt{G_2(d_k)}}{\sqrt{G_1(d_k)}}, +$$ + +$$ +r_{31}^{\mathrm{cal}}(d_k) += +20\log_{10} +\frac{\sqrt{G_3(d_k)}}{\sqrt{G_1(d_k)}}. +$$ + +Repeat this for each known calibration distance $d_k$ and build the calibration table: + +| Known distance from sensor 1 | $r_{21}^{\mathrm{cal}}$ | $r_{31}^{\mathrm{cal}}$ | +| ---------------------------: | ----------------------: | ----------------------: | +| $d_1$ | ... dB | ... dB | +| $d_2$ | ... dB | ... dB | +| $d_3$ | ... dB | ... dB | +| $\vdots$ | $\vdots$ | $\vdots$ | +| $d_N$ | ... dB | ... dB | + +This calibration represents the response to the known calibration input, not the mixed response of calibration plus unknown excitation. The approach assumes that the unknown excitation is not coherent with the calibration reference $u(t)$. Components uncorrelated with $u(t)$ are reduced by Welch averaging in the cross-spectral estimate. + +--- + +And the corresponding calculation with `np.fft.rfft` is essentially: + +```python +# For each Welch segment: +X_i = np.fft.rfft(window * x_i_segment, n=nfft) +U = np.fft.rfft(window * u_segment, n=nfft) + +S_xiu += X_i * np.conj(U) +S_uu += U * np.conj(U) + +# After all segments: +S_xiu /= n_segments +S_uu /= n_segments + +H_i = S_xiu / S_uu +``` + +Use the same segmentation, window, overlap, and `nfft` for all sensors and for the reference signal $u(t)$. The exact PSD scaling constants cancel in the transfer-function ratio, as long as they are applied consistently. + +--- + +## 7. Compute the unknown spectral-ratio vector + +For the unknown excitation, compute Welch spectra using the same settings and the same three-bin peak-band rule. + +Obtain + +$$ +r_{21}^{\mathrm{unknown}}, +\qquad +r_{31}^{\mathrm{unknown}}. +$$ + +Define the unknown ratio vector as + +$$ +\mathbf r^{\mathrm{unknown}} += +\begin{bmatrix} +r_{21}^{\mathrm{unknown}} \\ +r_{31}^{\mathrm{unknown}} +\end{bmatrix}. +$$ + +For each calibration distance, define + +$$ +\mathbf r^{\mathrm{cal}}(d_k) += +\begin{bmatrix} +r_{21}^{\mathrm{cal}}(d_k) \\ +r_{31}^{\mathrm{cal}}(d_k) +\end{bmatrix}. +$$ + +--- + +## 8. Match the unknown measurement to the calibration table + +Use the squared error in dB-ratio space: + +$$ +J(d_k)= +\left(r_{21}^{\mathrm{unknown}}-r_{21}^{\mathrm{cal}}(d_k)\right)^2+ +\left(r_{31}^{\mathrm{unknown}}-r_{31}^{\mathrm{cal}}(d_k)\right)^2. +$$ + +Select the distance with the smallest error: + +$$ +\hat d = \arg\min_{d_k} J(d_k). +$$ + +Report + +$$ +\boxed{\hat d}. +$$ + +This is nearest-neighbor matching in the two-dimensional feature space + +$$ +(r_{21},r_{31}). +$$ + +--- + +## 9. Uncertainty-weighted matching when repeated calibration data are available + +If repeated calibration measurements are available, use an uncertainty-weighted cost function instead of the unweighted squared error. + +For each calibration distance $d_k$, define the mean calibration vector + +$$ +\overline{\mathbf r}^{\mathrm{cal}}(d_k) += +\begin{bmatrix} +\bar r_{21}^{\mathrm{cal}}(d_k) \\ +\bar r_{31}^{\mathrm{cal}}(d_k) +\end{bmatrix}. +$$ + +Define the difference vector + +$$ +\Delta\mathbf r(d_k) += +\mathbf r^{\mathrm{unknown}} +- +\overline{\mathbf r}^{\mathrm{cal}}(d_k). +$$ + +Then use the covariance-weighted cost + +$$ +J(d_k) += +\Delta\mathbf r(d_k)^{\mathrm T} +\Sigma(d_k)^{-1} +\Delta\mathbf r(d_k), +$$ + +where $\Sigma(d_k)$ is the covariance matrix of the repeated calibration ratio measurements at distance $d_k$. + +If only separate standard deviations are estimated and the two ratios are treated as independent, use + +$$ +J(d_k)= +\frac{\left(r_{21}^{\mathrm{unknown}}-\bar r_{21}^{\mathrm{cal}}(d_k)\right)^2}{\sigma_{21}^2(d_k)} ++ +\frac{\left(r_{31}^{\mathrm{unknown}}-\bar r_{31}^{\mathrm{cal}}(d_k)\right)^2}{\sigma_{31}^2(d_k)}. +$$ + +This version is more defensible when one ratio is noisier or less repeatable than the other. A ratio with larger calibration scatter receives less weight in the distance estimate. + +--- + +## 10. Interpretation of the result + +The estimate + +$$ +\hat d +$$ + +should be interpreted as the calibrated distance whose relative spectral levels best match the unknown measurement. + +If the excitation is constrained to a known one-dimensional path, line, rail, beam, or test coordinate, then $\hat d$ can be interpreted as an estimated distance from sensor 1 along that path. + +If the excitation can occur at arbitrary two-dimensional or three-dimensional positions, then the same distance from sensor 1 may produce different ratio pairs depending on direction. Conversely, different positions may produce similar ratio pairs. In that case, the method should be described as spectral-ratio localization or calibration matching, not simply as distance estimation. + +--- + +## 11. Practical validation checks + +Before relying on the method, check the following: + +1. **Repeatability:** Repeated measurements at the same distance should produce similar $(r_{21},r_{31})$ values. + +2. **Separability:** Different distances should produce ratio vectors that are separated by more than the measurement scatter. + +3. **Uniqueness:** No two calibration distances should produce nearly identical ratio pairs. + +4. **Peak stability:** The selected frequency should remain at approximately the same frequency bin and should not merge with other peaks. + +5. **Reference quality:** Sensor 1 must show a clear, stable response at the chosen frequency. It should not be near a spectral node at the chosen peak. If $A_1$ or $G_1$ is very small, noisy, or not clearly peak-like, the ratios become unstable. + +6. **Known-signal coherence:** During calibration, the sensor responses should be coherent with the known calibration reference $u(t)$ in the chosen band. + +7. **Consistent excitation:** The unknown excitation should be comparable to the calibration excitation in direction, coupling, and spectral content. + +--- + +## 12. Compact algorithm + +1. Use the known calibration frequency $f_0$ to choose the matching bin: + + $$ + m_0=\arg\min_m |f_m-f_0|. + $$ + +2. Use the three-bin peak band + + $$ + \mathcal B(m_0)=\{m_0-1,m_0,m_0+1\}. + $$ + +3. For each calibration distance $d_k$, compute cross spectra with the known calibration reference $u(t)$: + + $$ + S_{x_i u}(f;d_k)=\left\langle X_i(f;d_k)U^*(f)\right\rangle, + \qquad + S_{uu}(f)=\left\langle U(f)U^*(f)\right\rangle. + $$ + +4. Estimate the calibration transfer functions: + + $$ + \hat H_i(f;d_k)=\frac{S_{x_i u}(f;d_k)}{S_{uu}(f)}, + \qquad i=1,2,3. + $$ + +5. Compute calibration peak-band transfer magnitudes: + + $$ + G_i(d_k)=\sum_{m\in\mathcal B(m_0)}\left|\hat H_i(f_m;d_k)\right|^2, + \qquad i=1,2,3. + $$ + +6. Compute calibration ratios: + + $$ + r_{21}^{\mathrm{cal}}(d_k)=10\log_{10}\frac{G_2(d_k)}{G_1(d_k)}, + \qquad + r_{31}^{\mathrm{cal}}(d_k)=10\log_{10}\frac{G_3(d_k)}{G_1(d_k)}. + $$ + +7. For the unknown excitation, compute Welch spectra $P_1(f_m),P_2(f_m),P_3(f_m)$ using the same settings. + +8. Compute unknown peak-band powers: + + $$ + A_i=\sum_{m\in\mathcal B(m_0)}P_i(f_m)\Delta f, + \qquad i=1,2,3. + $$ + +9. Compute unknown relative dB ratios: + + $$ + r_{21}^{\mathrm{unknown}}=10\log_{10}\frac{A_2}{A_1}, + \qquad + r_{31}^{\mathrm{unknown}}=10\log_{10}\frac{A_3}{A_1}. + $$ + +10. Compare to the calibration table using + + $$ + J(d_k)= + \left(r_{21}^{\mathrm{unknown}}-r_{21}^{\mathrm{cal}}(d_k)\right)^2+ + \left(r_{31}^{\mathrm{unknown}}-r_{31}^{\mathrm{cal}}(d_k)\right)^2. + $$ + +11. Estimate the distance as + + $$ + \boxed{ + \hat d=\arg\min_{d_k}J(d_k) + }. + $$ + +If repeated calibration data are available, use the uncertainty-weighted cost from Section 9. + +--- + +# Appendix A: Theoretical justification + +For a narrow frequency region around a dominant excitation peak, the measured response at sensor $i$ can be approximated as + +$$ +X_i(f) = H_i(f,d)S(f) + N_i(f), +$$ + +where: + +* $S(f)$ is the unknown excitation spectrum, +* $H_i(f,d)$ is the transfer path from the excitation point to sensor $i$, +* $d$ is the calibrated distance coordinate from sensor 1, +* $N_i(f)$ is measurement noise. + +If the peak has high signal-to-noise ratio and the system behaves approximately linearly, then near the peak frequency + +$$ +P_i(f) \approx |H_i(f,d)|^2 P_S(f). +$$ + +Taking ratios relative to sensor 1 approximately cancels the unknown source strength: + +$$ +\frac{P_i(f)}{P_1(f)} +\approx +\frac{|H_i(f,d)|^2}{|H_1(f,d)|^2}. +$$ + +In dB form, + +$$ +10\log_{10}\frac{P_i(f)}{P_1(f)} +\approx +20\log_{10}\frac{|H_i(f,d)|}{|H_1(f,d)|}. +$$ + +For calibration, the known reference signal $u(t)$ provides the input used to estimate the transfer path: + +$$ +\hat H_i(f;d_k) += +\frac{S_{x_i u}(f;d_k)}{S_{uu}(f)}. +$$ + +Therefore, relative calibration transfer magnitudes can be compared to the unknown relative spectral levels, provided the calibration and unknown measurements are comparable. + +The method is theoretically justifiable under the following assumptions: + +1. The structure and sensors behave approximately linearly and time-invariantly during calibration and testing. + +2. The calibration reference $u(t)$ is known, measured, and coherent with the sensor responses during calibration. + +3. The unknown excitation is comparable to the calibration excitations. + +4. The same peak frequency region is used for calibration and testing. + +5. The selected peak has high signal-to-noise ratio. + +6. Sensor gains and mountings remain fixed or are calibrated out. + +7. Sensor 1 is a valid reference sensor: it has a clear, stable, high-SNR peak at the selected frequency in both calibration and unknown measurements. + +8. The source position is constrained so that the calibrated distance coordinate is meaningful. + +9. The mapping from distance to $(r_{21},r_{31})$ is repeatable and sufficiently unique. + +Under these assumptions, the method is not just a heuristic. It is a calibrated inverse problem using relative spectral magnitudes. + +--- + +# Appendix B: Alternative peak-band definitions + +The main workflow uses a three-bin peak band because the peaks are assumed to be sharp, clean, and isolated. Other definitions are possible, but the chosen rule should be fixed before comparing calibration and unknown data. + +## B.1 Single-bin peak level + +Use + +$$ +\mathcal B(m_0)=\{m_0\}. +$$ + +This is the simplest approach, but it is sensitive to small peak shifts and bin alignment. It is usually not recommended unless the peak frequency is extremely stable and the FFT bin alignment is identical across all measurements. + +## B.2 Three-bin peak band + +Use + +$$ +\mathcal{B}(m_0)=\{m_0-1,m_0,m_0+1\}. +$$ + +This is the recommended default for sharp, isolated, visually clean peaks. It protects against one-bin jitter while still representing the peak itself. + +## B.3 Five-bin peak band + +Use + +$$ +\mathcal{B}(m_0)=\{m_0-2,m_0-1,m_0,m_0+1,m_0+2\}. +$$ + +This can be useful if the peak location shifts slightly more between measurements. It should only be used if the neighboring bins still belong to the same peak and not to noise or another mode. + +## B.4 Half-power peak support + +Define the band as the contiguous bins around $m_0$ for which the sensor-1 Welch level remains within 3 dB of the peak: + +$$ +10\log_{10}P_1(f_m) +\geq +10\log_{10}P_1(f_{m_0}) - 3. +$$ + +This corresponds to the half-power region of the peak. For very sharp Welch peaks, this may reduce to a single bin. In that case, use at least a three-bin band to avoid excessive sensitivity to bin alignment: + +$$ +\mathcal{B}(m_0) += +\mathcal{B}_{3\text{-bin}}(m_0) +\cup +\mathcal{B}_{\text{half-power}}(m_0). +$$ + +For the present workflow, where the peaks appear very sharp and clear, the recommended rule remains the three-bin peak band.