diff --git a/crates/processing/src/craft.rs b/crates/processing/src/craft.rs index b7d441b..fb0f326 100644 --- a/crates/processing/src/craft.rs +++ b/crates/processing/src/craft.rs @@ -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, @@ -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, @@ -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(); @@ -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, @@ -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(); @@ -568,7 +576,11 @@ fn fit_region( .copied() .collect(); let times: Vec = (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 { diff --git a/crates/processing/src/craft_tests.rs b/crates/processing/src/craft_tests.rs index ce84017..03eff81 100644 --- a/crates/processing/src/craft_tests.rs +++ b/crates/processing/src/craft_tests.rs @@ -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); diff --git a/docs/src/content/docs/guides/craft.md b/docs/src/content/docs/guides/craft.md index 483a332..795369b 100644 --- a/docs/src/content/docs/guides/craft.md +++ b/docs/src/content/docs/guides/craft.md @@ -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 diff --git a/docs/src/content/docs/zh-cn/guides/craft.md b/docs/src/content/docs/zh-cn/guides/craft.md index bf2b2d3..3e4be97 100644 --- a/docs/src/content/docs/zh-cn/guides/craft.md +++ b/docs/src/content/docs/zh-cn/guides/craft.md @@ -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。它只决定宽信号组如何分段计算,不会增加