Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
22 changes: 17 additions & 5 deletions crates/processing/src/craft.rs
Original file line number Diff line number Diff line change
Expand Up @@ -130,7 +130,7 @@ impl CraftParams {
Self {
profile: CraftProfile::Conventional,
regions: Vec::new(),
max_components_per_fit_window: 7,
max_components_per_fit_window: 15,
min_amplitude_to_noise: 3.3,
linewidth_hz: (0.05, 10.0),
filter_taps: 499,
Expand Down Expand Up @@ -329,7 +329,7 @@ pub fn process_craft_cancellable(
if cancelled() {
return Err(CraftError::Cancelled);
}
let result = fit_region(input, skip, sw, region, params, cancelled)?;
let result = fit_region(input, skip, data.group_delay, sw, region, params, cancelled)?;
if let Some((kind, message)) = result.warning {
warnings.push(CraftWarning {
kind,
Expand Down Expand Up @@ -455,7 +455,10 @@ pub fn process_craft_cancellable(
.points
.iter()
.enumerate()
.map(|(index, &sample)| sample - model_at(&components, index as f64 / sw))
.map(|(index, &sample)| {
let time = (index as f64 - data.group_delay) / sw;
sample - model_at(&components, time)
})
.collect();
let residual_rss: f64 = residual_fid[skip..].iter().map(Complex64::norm_sqr).sum();
let input_rss: f64 = input.iter().map(Complex64::norm_sqr).sum();
Expand Down Expand Up @@ -518,6 +521,7 @@ pub fn process_craft_cancellable(
fn fit_region(
input: &[Complex64],
skipped_points: usize,
group_delay_points: f64,
sw: f64,
region: HzRegion,
params: &CraftParams,
Expand All @@ -535,7 +539,11 @@ fn fit_region(
.iter()
.enumerate()
.map(|(index, &value)| {
let time = (skipped_points + index) as f64 / sw;
// Bruker points after the digital-filter transient are already
// samples of the FID starting at `ceil(delay) - delay`. Keeping the
// raw point number here would reintroduce the removed delay as an
// amplitude extrapolation and a first-order phase ramp.
let time = (skipped_points as f64 + index as f64 - group_delay_points) / sw;
value * Complex64::from_polar(1.0, -TAU * center_hz * time)
})
.collect();
Expand Down Expand Up @@ -568,7 +576,11 @@ fn fit_region(
.copied()
.collect();
let times: Vec<f64> = (0..samples.len())
.map(|index| (skipped_points + phase_start + index * decimation) as f64 / sw)
.map(|index| {
(skipped_points as f64 + phase_start as f64 + (index * decimation) as f64
- group_delay_points)
/ sw
})
.collect();
if samples.len() < 16 {
return Ok(RegionResult {
Expand Down
70 changes: 70 additions & 0 deletions crates/processing/src/craft_tests.rs
Original file line number Diff line number Diff line change
Expand Up @@ -98,6 +98,76 @@ fn ssfp_skip_extrapolates_amplitude_to_time_zero() {
);
}

#[test]
fn group_delay_uses_the_physical_fid_time_origin() {
let sw = 2_000.0;
let delay = 67.985_885_620_117_2;
let expected_amplitude = 5.0;
let expected_phase = 0.4;
let mut input = data(&[], 4096, sw);
input.group_delay = delay;
input.points = (0..input.points.len())
.map(|index| {
let time = (index as f64 - delay) / sw;
Complex64::from_polar(
expected_amplitude * (-PI * 2.0 * time).exp(),
expected_phase + TAU * 100.0 * time,
)
})
.collect();
let params = CraftParams {
regions: vec![CraftRegion::new(CraftRegionId(1), 0.18, 0.22)],
filter_taps: 127,
..CraftParams::default()
};

let result = process_craft_cancellable(
&input,
&CraftInvocation::acquisition(&input, params),
&|| false,
)
.unwrap();

assert_eq!(result.components.len(), 1, "{:?}", result.components);
let component = &result.components[0];
assert!((component.amplitude_t0 - expected_amplitude).abs() < 0.02);
let phase_error = (component.phase_rad - expected_phase).sin().abs();
assert!(phase_error < 0.01, "phase error was {phase_error}");
}

#[test]
fn default_model_limit_resolves_non_lorentzian_multiplet_quantitation() {
let mut components = vec![
(-108.0, 0.75, 0.3, 0.5),
(-100.0, 1.50, 0.3, 0.6),
(-92.0, 0.75, 0.3, 0.7),
];
for (frequency, amplitude) in [(80.0, 0.25), (88.0, 0.75), (96.0, 0.75), (104.0, 0.25)] {
components.push((frequency, amplitude * 0.6, 0.3, 1.0));
components.push((frequency + 0.8, amplitude * 0.4, 0.3, 1.8));
}
let input = data(&components, 4096, 2_000.0);
let params = CraftParams {
regions: vec![
CraftRegion::new(CraftRegionId(0), -0.23, -0.17),
CraftRegion::new(CraftRegionId(1), 0.14, 0.23),
],
filter_taps: 127,
..CraftParams::default()
};

let result = process_craft_cancellable(
&input,
&CraftInvocation::acquisition(&input, params),
&|| false,
)
.unwrap();

assert!(result.diagnostics.fit_windows[1].selected_model_order > 7);
let ratio = result.region_ratio.unwrap().value;
assert!((ratio - 1.5).abs() < 0.03, "ratio was {ratio}");
}

#[test]
fn rejects_frequency_domain_input() {
let mut input = data(&[], 128, 1000.0);
Expand Down
2 changes: 1 addition & 1 deletion docs/src/content/docs/guides/craft.md
Original file line number Diff line number Diff line change
Expand Up @@ -60,7 +60,7 @@ profile uses these defaults:

- **Minimum A/N**: 3.3. Lower values retain weaker candidates but increase the
chance of fitting noise; values below 3.3 are flagged for review.
- **Max components / fit window**: 7 (allowed range 1–64). Reaching the limit
- **Max components / fit window**: 15 (allowed range 1–64). Reaching the limit
is reported as a diagnostic warning.
- **Linewidth range (Hz)**: 0.05–10 Hz.
- **Fit window width (Hz)**: 500 Hz. This controls how wide groups are divided
Expand Down
2 changes: 1 addition & 1 deletion docs/src/content/docs/zh-cn/guides/craft.md
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,7 @@ CRAFT 直接拟合一维 NMR 采集中的原始复数 FID,并报告共振分

- **Minimum A/N**:3.3。较低值会保留更弱的候选分量,但也更容易拟合噪声;低于
3.3 时会标记为需要复核。
- **Max components / fit window**:7(允许范围 1–64)。达到上限会在诊断中给出
- **Max components / fit window**:15(允许范围 1–64)。达到上限会在诊断中给出
警告。
- **Linewidth range (Hz)**:0.05–10 Hz。
- **Fit window width (Hz)**:500 Hz。它只决定宽信号组如何分段计算,不会增加
Expand Down
Loading