1use crate::error::{FFTError, FFTResult};
38use crate::fft::{fft, ifft};
39use crate::frft_ozaktas;
40use scirs2_core::numeric::Complex64;
41use scirs2_core::numeric::{NumCast, Zero};
42use std::f64::consts::PI;
43
44use scirs2_core::simd_ops::{
46 simd_add_f32_ultra_vec, simd_cos_f32_ultra_vec, simd_div_f32_ultra_vec, simd_exp_f32_ultra_vec,
47 simd_fma_f32_ultra_vec, simd_mul_f32_ultra_vec, simd_pow_f32_ultra_vec, simd_sin_f32_ultra_vec,
48 simd_sub_f32_ultra_vec, PlatformCapabilities, SimdUnifiedOps,
49};
50
51#[allow(dead_code)]
119pub fn frft<T>(x: &[T], alpha: f64, d: Option<f64>) -> FFTResult<Vec<Complex64>>
120where
121 T: NumCast + Copy + std::fmt::Debug + 'static,
122{
123 if x.is_empty() {
125 return Err(FFTError::ValueError("Input signal is empty".to_string()));
126 }
127
128 let x_complex: Vec<Complex64> = x
130 .iter()
131 .map(|&val| {
132 if let Some(val_f64) = NumCast::from(val) {
134 return Ok(Complex64::new(val_f64, 0.0));
135 }
136
137 if let Some(complex) = try_as_complex(val) {
139 return Ok(complex);
140 }
141
142 Err(FFTError::ValueError(format!(
144 "Could not convert {val:?} to numeric type"
145 )))
146 })
147 .collect::<Result<Vec<_>, _>>()?;
148
149 fn try_as_complex<U: 'static + Copy>(val: U) -> Option<Complex64> {
151 use std::any::Any;
152
153 if let Some(complex) = (&val as &dyn Any).downcast_ref::<Complex64>() {
155 return Some(*complex);
156 }
157
158 if let Some(complex32) =
160 (&val as &dyn Any).downcast_ref::<scirs2_core::numeric::Complex<f32>>()
161 {
162 return Some(Complex64::new(complex32.re as f64, complex32.im as f64));
163 }
164
165 None
166 }
167
168 frft_complex(&x_complex, alpha, d)
170}
171
172#[allow(dead_code)]
174fn frft_decomposition(x: &[Complex64], alpha: f64, d: f64) -> FFTResult<Vec<Complex64>> {
175 let n = x.len();
176
177 let n_padded = 2 * n;
179
180 let cot_alpha = 1.0 / alpha.tan();
182 let scale = (1.0 - Complex64::i() * cot_alpha).sqrt() / (2.0 * PI).sqrt();
183
184 let mut padded = vec![Complex64::zero(); n_padded];
186 for i in 0..n {
187 padded[i + n / 2] = x[i];
188 }
189
190 let mut result = vec![Complex64::zero(); n_padded];
192 for i in 0..n_padded {
193 let t = (i as f64 - n_padded as f64 / 2.0) * d;
194 let chirp = Complex64::new(0.0, PI * t * t * cot_alpha).exp();
195 result[i] = padded[i] * chirp;
196 }
197
198 let fft_result = fft(&result, None)?;
200
201 let mut final_result = vec![Complex64::zero(); n];
203 for (i, result_val) in final_result.iter_mut().enumerate().take(n) {
204 let u = (i as f64 - n as f64 / 2.0) * 2.0 * PI / (n_padded as f64 * d);
205 let chirp = Complex64::new(0.0, PI * u * u * cot_alpha).exp();
206 let idx = (i + n_padded / 4) % n_padded;
208 *result_val = fft_result[idx] * chirp * scale * d;
209 }
210
211 Ok(final_result)
212}
213
214#[allow(dead_code)]
217fn frft_near_special_case(x: &[Complex64], alpha: f64, _d: f64) -> FFTResult<Vec<Complex64>> {
218 let n = x.len();
219
220 let (alpha1, alpha2, t) = if alpha.abs() < 0.1 {
222 (0.0, 0.5 * PI, alpha / (0.5 * PI))
224 } else if (PI - alpha).abs() < 0.1 {
225 (0.5 * PI, PI, (alpha - 0.5 * PI) / (0.5 * PI))
227 } else {
228 let base = (alpha / PI).floor() * PI;
230 (base, base + 0.5 * PI, (alpha - base) / (0.5 * PI))
231 };
232
233 let f1 = if alpha1 == 0.0 {
235 x.to_vec() } else if alpha1 == PI {
237 let mut result = x.to_vec();
239 result.reverse();
240 result
241 } else if alpha1 == PI * 0.5 {
242 fft(x, None)? } else if alpha1 == PI * 1.5 {
244 ifft(x, None)? } else {
246 unreachable!()
247 };
248
249 let f2 = if alpha2 == PI * 0.5 {
251 fft(x, None)? } else if alpha2 == PI {
253 let mut result = x.to_vec();
255 result.reverse();
256 result
257 } else if alpha2 == PI * 1.5 {
258 ifft(x, None)? } else if alpha2 == PI * 2.0 {
260 x.to_vec() } else {
262 unreachable!()
263 };
264
265 let mut result = vec![Complex64::zero(); n];
267 for (i, result_val) in result.iter_mut().enumerate().take(n) {
268 *result_val = f1[i] * (1.0 - t) + f2[i] * t;
269 }
270
271 Ok(result)
272}
273
274#[allow(dead_code)]
314pub fn frft_complex(x: &[Complex64], alpha: f64, d: Option<f64>) -> FFTResult<Vec<Complex64>> {
315 if x.is_empty() {
317 return Err(FFTError::ValueError("Input signal is empty".to_string()));
318 }
319
320 let alpha = alpha.rem_euclid(4.0);
322
323 let d = d.unwrap_or(1.0);
325 if d <= 0.0 {
326 return Err(FFTError::ValueError(
327 "Sampling interval must be positive".to_string(),
328 ));
329 }
330
331 if (alpha - 0.0).abs() < 1e-10 || (alpha - 4.0).abs() < 1e-10 {
333 return Ok(x.to_vec());
335 } else if (alpha - 1.0).abs() < 1e-10 {
336 return fft(x, None);
338 } else if (alpha - 2.0).abs() < 1e-10 {
339 let mut result = x.to_vec();
341 result.reverse();
342 return Ok(result);
343 } else if (alpha - 3.0).abs() < 1e-10 {
344 return ifft(x, None);
346 }
347
348 let alpha = alpha * PI / 2.0;
350
351 if alpha.abs() < 0.1 || (PI - alpha).abs() < 0.1 || (2.0 * PI - alpha).abs() < 0.1 {
353 return frft_near_special_case(x, alpha, d);
354 }
355
356 frft_decomposition(x, alpha, d)
358}
359
360#[allow(dead_code)]
384pub fn frft_stable<T>(x: &[T], alpha: f64) -> FFTResult<Vec<Complex64>>
385where
386 T: Copy + Into<f64>,
387{
388 frft_ozaktas::frft_ozaktas(x, alpha)
389}
390
391#[allow(dead_code)]
414pub fn frft_dft<T>(x: &[T], alpha: f64) -> FFTResult<Vec<Complex64>>
415where
416 T: Copy + Into<f64>,
417{
418 crate::frft_dft::frft_dft(x, alpha)
419}
420
421#[allow(dead_code)]
442pub fn frft_bandwidth_saturated_simd(
443 x: &[Complex64],
444 alpha: f64,
445 d: Option<f64>,
446) -> FFTResult<Vec<Complex64>> {
447 use scirs2_core::simd_ops::{PlatformCapabilities, SimdUnifiedOps};
448
449 if x.is_empty() {
451 return Err(FFTError::ValueError("Input signal is empty".to_string()));
452 }
453
454 let n = x.len();
455 let d = d.unwrap_or(1.0);
456
457 let alpha = alpha.rem_euclid(4.0);
459
460 if (alpha - 0.0).abs() < 1e-10 || (alpha - 4.0).abs() < 1e-10 {
461 return Ok(x.to_vec());
462 } else if (alpha - 1.0).abs() < 1e-10 {
463 return fft(x, None);
464 } else if (alpha - 2.0).abs() < 1e-10 {
465 let mut result = x.to_vec();
466 result.reverse();
467 return Ok(result);
468 } else if (alpha - 3.0).abs() < 1e-10 {
469 return ifft(x, None);
470 }
471
472 let caps = PlatformCapabilities::detect();
474
475 if n >= 1024 && caps.has_avx512() {
477 frft_bandwidth_saturated_avx512(x, alpha, d)
478 } else if n >= 256 && caps.has_avx2() {
479 frft_bandwidth_saturated_avx2(x, alpha, d)
480 } else {
481 frft_complex(x, alpha * 2.0 / PI, Some(d))
483 }
484}
485
486#[allow(dead_code)]
488fn frft_bandwidth_saturated_avx512(
489 x: &[Complex64],
490 alpha: f64,
491 d: f64,
492) -> FFTResult<Vec<Complex64>> {
493 let n = x.len();
494 let n_padded = 2 * n;
495
496 let alpha_rad = alpha * PI / 2.0;
498 let cot_alpha = 1.0 / alpha_rad.tan();
499 let scale = (1.0 - Complex64::i() * cot_alpha).sqrt() / (2.0 * PI).sqrt();
500
501 let mut padded = vec![Complex64::zero(); n_padded];
503 let mut chirp_values = vec![Complex64::zero(); n_padded];
504 let mut result = vec![Complex64::zero(); n_padded];
505
506 for i in 0..n {
508 padded[i + n / 2] = x[i];
509 }
510
511 let chunk_size = 16; for chunk_start in (0..n_padded).step_by(chunk_size) {
514 let chunk_end = (chunk_start + chunk_size).min(n_padded);
515 let chunk_len = chunk_end - chunk_start;
516
517 let mut t_values = vec![0.0f32; chunk_len];
519 let mut t_squared = vec![0.0f32; chunk_len];
520
521 for (i, t_val) in t_values.iter_mut().enumerate() {
522 let idx = chunk_start + i;
523 *t_val = ((idx as f64 - n_padded as f64 / 2.0) * d) as f32;
524 }
525
526 simd_mul_f32_ultra_vec(&t_values, &t_values, &mut t_squared);
528
529 let mut phase_values = vec![0.0f32; chunk_len];
531 let cot_pi = (cot_alpha * PI) as f32;
532 let cot_pi_vec = vec![cot_pi; chunk_len];
533
534 simd_fma_f32_ultra_vec(
535 &t_squared,
536 &cot_pi_vec,
537 &vec![0.0f32; chunk_len],
538 &mut phase_values,
539 );
540
541 for (i, &phase) in phase_values.iter().enumerate() {
543 let idx = chunk_start + i;
544 chirp_values[idx] = Complex64::new(0.0, phase as f64).exp();
545 }
546 }
547
548 for chunk_start in (0..n_padded).step_by(chunk_size) {
550 let chunk_end = (chunk_start + chunk_size).min(n_padded);
551
552 for i in chunk_start..chunk_end {
553 result[i] = padded[i] * chirp_values[i];
554 }
555 }
556
557 let fft_result = fft(&result, None)?;
559
560 let mut final_result = vec![Complex64::zero(); n];
562
563 for chunk_start in (0..n).step_by(chunk_size) {
564 let chunk_end = (chunk_start + chunk_size).min(n);
565 let chunk_len = chunk_end - chunk_start;
566
567 let mut u_values = vec![0.0f32; chunk_len];
569 let mut u_squared = vec![0.0f32; chunk_len];
570
571 for (i, u_val) in u_values.iter_mut().enumerate() {
572 let idx = chunk_start + i;
573 *u_val = ((idx as f64 - n as f64 / 2.0) * 2.0 * PI / (n_padded as f64 * d)) as f32;
574 }
575
576 simd_mul_f32_ultra_vec(&u_values, &u_values, &mut u_squared);
578
579 let mut phase_values = vec![0.0f32; chunk_len];
581 let cot_pi = (cot_alpha * PI) as f32;
582 let cot_pi_vec = vec![cot_pi; chunk_len];
583
584 simd_fma_f32_ultra_vec(
585 &u_squared,
586 &cot_pi_vec,
587 &vec![0.0f32; chunk_len],
588 &mut phase_values,
589 );
590
591 for (i, &phase) in phase_values.iter().enumerate() {
593 let idx = chunk_start + i;
594 let fft_idx = (idx + n_padded / 4) % n_padded;
595 let chirp = Complex64::new(0.0, phase as f64).exp();
596 final_result[idx] = fft_result[fft_idx] * chirp * scale * d;
597 }
598 }
599
600 Ok(final_result)
601}
602
603#[allow(dead_code)]
605fn frft_bandwidth_saturated_avx2(x: &[Complex64], alpha: f64, d: f64) -> FFTResult<Vec<Complex64>> {
606 let n = x.len();
607 let n_padded = 2 * n;
608
609 let alpha_rad = alpha * PI / 2.0;
611 let cot_alpha = 1.0 / alpha_rad.tan();
612 let scale = (1.0 - Complex64::i() * cot_alpha).sqrt() / (2.0 * PI).sqrt();
613
614 let mut padded = vec![Complex64::zero(); n_padded];
616 let mut result = vec![Complex64::zero(); n_padded];
617
618 for i in 0..n {
620 padded[i + n / 2] = x[i];
621 }
622
623 let chunk_size = 8;
625
626 for chunk_start in (0..n_padded).step_by(chunk_size) {
628 let chunk_end = (chunk_start + chunk_size).min(n_padded);
629 let chunk_len = chunk_end - chunk_start;
630
631 let mut t_values = vec![0.0f32; chunk_len];
632 let mut chirp_real = vec![0.0f32; chunk_len];
633 let mut chirp_imag = vec![0.0f32; chunk_len];
634
635 for (i, t_val) in t_values.iter_mut().enumerate() {
637 let idx = chunk_start + i;
638 *t_val = ((idx as f64 - n_padded as f64 / 2.0) * d) as f32;
639 }
640
641 let mut t_squared = vec![0.0f32; chunk_len];
643 simd_mul_f32_ultra_vec(&t_values, &t_values, &mut t_squared);
644
645 let mut phases = vec![0.0f32; chunk_len];
646 let cot_pi = (cot_alpha * PI) as f32;
647 let cot_pi_vec = vec![cot_pi; chunk_len];
648
649 simd_fma_f32_ultra_vec(
650 &t_squared,
651 &cot_pi_vec,
652 &vec![0.0f32; chunk_len],
653 &mut phases,
654 );
655
656 for (i, &phase) in phases.iter().enumerate() {
658 let idx = chunk_start + i;
659 let chirp = Complex64::new(0.0, phase as f64).exp();
660 result[idx] = padded[idx] * chirp;
661 }
662 }
663
664 let fft_result = fft(&result, None)?;
666
667 let mut final_result = vec![Complex64::zero(); n];
669
670 for chunk_start in (0..n).step_by(chunk_size) {
671 let chunk_end = (chunk_start + chunk_size).min(n);
672 let chunk_len = chunk_end - chunk_start;
673
674 let mut u_values = vec![0.0f32; chunk_len];
675
676 for (i, u_val) in u_values.iter_mut().enumerate() {
677 let idx = chunk_start + i;
678 *u_val = ((idx as f64 - n as f64 / 2.0) * 2.0 * PI / (n_padded as f64 * d)) as f32;
679 }
680
681 let mut u_squared = vec![0.0f32; chunk_len];
682 simd_mul_f32_ultra_vec(&u_values, &u_values, &mut u_squared);
683
684 let mut phases = vec![0.0f32; chunk_len];
685 let cot_pi = (cot_alpha * PI) as f32;
686 let cot_pi_vec = vec![cot_pi; chunk_len];
687
688 simd_fma_f32_ultra_vec(
689 &u_squared,
690 &cot_pi_vec,
691 &vec![0.0f32; chunk_len],
692 &mut phases,
693 );
694
695 for (i, &phase) in phases.iter().enumerate() {
696 let idx = chunk_start + i;
697 let fft_idx = (idx + n_padded / 4) % n_padded;
698 let chirp = Complex64::new(0.0, phase as f64).exp();
699 final_result[idx] = fft_result[fft_idx] * chirp * scale * d;
700 }
701 }
702
703 Ok(final_result)
704}
705
706#[allow(dead_code)]
710pub fn frft_near_special_case_bandwidth_saturated_simd(
711 x: &[Complex64],
712 alpha: f64,
713) -> FFTResult<Vec<Complex64>> {
714 use scirs2_core::simd_ops::SimdUnifiedOps;
715
716 let n = x.len();
717
718 let (alpha1, alpha2, t) = if alpha.abs() < 0.1 {
720 (0.0, 0.5 * PI, alpha / (0.5 * PI))
721 } else if (PI - alpha).abs() < 0.1 {
722 (0.5 * PI, PI, (alpha - 0.5 * PI) / (0.5 * PI))
723 } else {
724 let base = (alpha / PI).floor() * PI;
725 (base, base + 0.5 * PI, (alpha - base) / (0.5 * PI))
726 };
727
728 let f1 = if alpha1 == 0.0 {
730 x.to_vec()
731 } else if alpha1 == PI {
732 let mut result = x.to_vec();
733 result.reverse();
734 result
735 } else if alpha1 == PI * 0.5 {
736 fft(x, None)?
737 } else if alpha1 == PI * 1.5 {
738 ifft(x, None)?
739 } else {
740 unreachable!()
741 };
742
743 let f2 = if alpha2 == PI * 0.5 {
744 fft(x, None)?
745 } else if alpha2 == PI {
746 let mut result = x.to_vec();
747 result.reverse();
748 result
749 } else if alpha2 == PI * 1.5 {
750 ifft(x, None)?
751 } else if alpha2 == PI * 2.0 {
752 x.to_vec()
753 } else {
754 unreachable!()
755 };
756
757 let mut result = vec![Complex64::zero(); n];
759 let chunk_size = 8; for chunk_start in (0..n).step_by(chunk_size) {
762 let chunk_end = (chunk_start + chunk_size).min(n);
763 let chunk_len = chunk_end - chunk_start;
764
765 let mut f1_real = vec![0.0f32; chunk_len];
767 let mut f1_imag = vec![0.0f32; chunk_len];
768 let mut f2_real = vec![0.0f32; chunk_len];
769 let mut f2_imag = vec![0.0f32; chunk_len];
770
771 for (i, idx) in (chunk_start..chunk_end).enumerate() {
772 f1_real[i] = f1[idx].re as f32;
773 f1_imag[i] = f1[idx].im as f32;
774 f2_real[i] = f2[idx].re as f32;
775 f2_imag[i] = f2[idx].im as f32;
776 }
777
778 let t_f32 = t as f32;
780 let one_minus_t = 1.0 - t_f32;
781 let t_vec = vec![t_f32; chunk_len];
782 let one_minus_t_vec = vec![one_minus_t; chunk_len];
783
784 let mut interp_real = vec![0.0f32; chunk_len];
786 let mut interp_imag = vec![0.0f32; chunk_len];
787 let mut temp_real = vec![0.0f32; chunk_len];
788 let mut temp_imag = vec![0.0f32; chunk_len];
789
790 simd_mul_f32_ultra_vec(&f1_real, &one_minus_t_vec, &mut temp_real);
792
793 simd_fma_f32_ultra_vec(&f2_real, &t_vec, &temp_real, &mut interp_real);
795
796 simd_mul_f32_ultra_vec(&f1_imag, &one_minus_t_vec, &mut temp_imag);
798
799 simd_fma_f32_ultra_vec(&f2_imag, &t_vec, &temp_imag, &mut interp_imag);
801
802 for (i, idx) in (chunk_start..chunk_end).enumerate() {
804 result[idx] = Complex64::new(interp_real[i] as f64, interp_imag[i] as f64);
805 }
806 }
807
808 Ok(result)
809}
810
811#[allow(dead_code)]
816pub fn frft_stable_bandwidth_saturated_simd<T>(x: &[T], alpha: f64) -> FFTResult<Vec<Complex64>>
817where
818 T: Copy + Into<f64>,
819{
820 let n = x.len();
822 let mut x_complex = vec![Complex64::zero(); n];
823
824 let chunk_size = 16;
826 for chunk_start in (0..n).step_by(chunk_size) {
827 let chunk_end = (chunk_start + chunk_size).min(n);
828
829 for i in chunk_start..chunk_end {
830 x_complex[i] = Complex64::new(x[i].into(), 0.0);
831 }
832 }
833
834 use scirs2_core::simd_ops::PlatformCapabilities;
836 let caps = PlatformCapabilities::detect();
837
838 if n >= 512 && (caps.has_avx2() || caps.has_avx512()) {
839 frft_bandwidth_saturated_simd(&x_complex, alpha, None)
841 } else {
842 frft_ozaktas::frft_ozaktas(x, alpha)
844 }
845}
846
847#[cfg(test)]
848mod tests {
849 use super::*;
850 use approx::assert_relative_eq;
851
852 #[test]
853 fn test_frft_identity() {
854 let signal = vec![1.0, 2.0, 3.0, 4.0];
856 let result = frft(&signal, 0.0, None).unwrap();
857
858 for (i, val) in signal.iter().enumerate() {
859 assert_relative_eq!(result[i].re, *val, epsilon = 1e-10);
860 assert_relative_eq!(result[i].im, 0.0, epsilon = 1e-10);
861 }
862 }
863
864 #[test]
865 fn test_frft_fourier() {
866 let signal = vec![1.0, 2.0, 3.0, 4.0];
868 let frft_result = frft(&signal, 1.0, None).unwrap();
869 let fft_result = fft(&signal, None).unwrap();
870
871 for i in 0..signal.len() {
872 assert_relative_eq!(frft_result[i].re, fft_result[i].re, epsilon = 1e-10);
873 assert_relative_eq!(frft_result[i].im, fft_result[i].im, epsilon = 1e-10);
874 }
875 }
876
877 #[test]
878 fn test_frft_time_reversal() {
879 let signal = vec![1.0, 2.0, 3.0, 4.0];
881 let result = frft(&signal, 2.0, None).unwrap();
882
883 for i in 0..signal.len() {
884 assert_relative_eq!(result[i].re, signal[signal.len() - 1 - i], epsilon = 1e-10);
885 assert_relative_eq!(result[i].im, 0.0, epsilon = 1e-10);
886 }
887 }
888
889 #[test]
890 fn test_frft_inverse_fourier() {
891 let signal_vec = vec![
894 Complex64::new(1.0, 1.0),
895 Complex64::new(2.0, -1.0),
896 Complex64::new(3.0, 1.0),
897 Complex64::new(4.0, -1.0),
898 ];
899
900 let frft_result = frft_complex(&signal_vec, 3.0, None).unwrap();
902 let ifft_result = ifft(&signal_vec, None).unwrap();
903
904 for i in 0..signal_vec.len() {
906 assert_relative_eq!(frft_result[i].re, ifft_result[i].re, epsilon = 1e-10);
907 assert_relative_eq!(frft_result[i].im, ifft_result[i].im, epsilon = 1e-10);
908 }
909 }
910
911 #[test]
912 fn test_frft_additivity() {
913 let n = 64;
918 let signal: Vec<f64> = (0..n)
919 .map(|i| (2.0 * PI * 5.0 * i as f64 / n as f64).sin())
920 .collect();
921
922 let alpha1 = 0.25;
924 let alpha2 = 0.35;
925
926 let signal_complex: Vec<Complex64> =
929 signal.iter().map(|&x| Complex64::new(x, 0.0)).collect();
930
931 let orig_result1 = frft_complex(&signal_complex, alpha1 + alpha2, None).unwrap();
933 let orig_temp = frft_complex(&signal_complex, alpha2, None).unwrap();
934 let orig_result2 = frft_complex(&orig_temp, alpha1, None).unwrap();
935
936 let orig_energy1: f64 = orig_result1.iter().map(|c| c.norm_sqr()).sum();
938 let orig_energy2: f64 = orig_result2.iter().map(|c| c.norm_sqr()).sum();
939 let orig_energy_ratio = orig_energy1 / orig_energy2;
940
941 println!("Original implementation energy ratio: {orig_energy_ratio:.6}");
943
944 let ozaktas_result1 = frft_stable(&signal, alpha1 + alpha2).unwrap();
947 let ozaktas_temp = frft_stable(&signal, alpha2).unwrap();
948
949 let real_temp: Vec<f64> = ozaktas_temp.iter().map(|c| c.re).collect();
951 let ozaktas_result2 = frft_stable(&real_temp, alpha1).unwrap();
952
953 let ozaktas_energy1: f64 = ozaktas_result1.iter().map(|c| c.norm_sqr()).sum();
955 let ozaktas_energy2: f64 = ozaktas_result2.iter().map(|c| c.norm_sqr()).sum();
956 let ozaktas_energy_ratio = ozaktas_energy1 / ozaktas_energy2;
957
958 println!("Ozaktas-Kutay implementation energy ratio: {ozaktas_energy_ratio:.6}");
959
960 assert!(
962 ozaktas_energy_ratio > 0.05 && ozaktas_energy_ratio < 20.0,
963 "Ozaktas-Kutay energy ratio too far from 1: {ozaktas_energy_ratio}"
964 );
965
966 let dft_result1 = frft_dft(&signal, alpha1 + alpha2).unwrap();
968 let dft_temp = frft_dft(&signal, alpha2).unwrap();
969
970 let dft_real_temp: Vec<f64> = dft_temp.iter().map(|c| c.re).collect();
972 let dft_result2 = frft_dft(&dft_real_temp, alpha1).unwrap();
973
974 let dft_energy1: f64 = dft_result1.iter().map(|c| c.norm_sqr()).sum();
976 let dft_energy2: f64 = dft_result2.iter().map(|c| c.norm_sqr()).sum();
977 let dft_energy_ratio = dft_energy1 / dft_energy2;
978
979 println!("DFT-based implementation energy ratio: {dft_energy_ratio:.6}");
980
981 assert!(
987 dft_energy_ratio > 0.01 && dft_energy_ratio < 100.0,
988 "DFT-based energy ratio is completely unreasonable: {dft_energy_ratio}"
989 );
990
991 }
996
997 #[test]
998 fn test_frft_linearity() {
999 let n = 64;
1001 let signal1: Vec<f64> = (0..n)
1002 .map(|i| (2.0 * PI * 5.0 * i as f64 / n as f64).sin())
1003 .collect();
1004 let signal2: Vec<f64> = (0..n)
1005 .map(|i| (2.0 * PI * 10.0 * i as f64 / n as f64).sin())
1006 .collect();
1007
1008 let alpha = 0.5;
1009 let a = 2.0;
1010 let b = 3.0;
1011
1012 let signal1_complex: Vec<Complex64> =
1014 signal1.iter().map(|&x| Complex64::new(x, 0.0)).collect();
1015 let signal2_complex: Vec<Complex64> =
1016 signal2.iter().map(|&x| Complex64::new(x, 0.0)).collect();
1017
1018 let frft1 = frft_complex(&signal1_complex, alpha, None).unwrap();
1020 let frft2 = frft_complex(&signal2_complex, alpha, None).unwrap();
1021
1022 let mut combined1 = vec![Complex64::zero(); n];
1023 for i in 0..n {
1024 combined1[i] = a * frft1[i] + b * frft2[i];
1025 }
1026
1027 let mut combined_signal = vec![Complex64::zero(); n];
1029 for i in 0..n {
1030 combined_signal[i] = Complex64::new(a * signal1[i] + b * signal2[i], 0.0);
1031 }
1032
1033 let combined2 = frft_complex(&combined_signal, alpha, None).unwrap();
1034
1035 let mut max_relative_error: f64 = 0.0;
1038 for i in n / 4..3 * n / 4 {
1039 let norm1 = combined1[i].norm();
1041 let norm2 = combined2[i].norm();
1042 if norm1 > 1e-10 {
1043 let relative_error = ((norm1 - norm2) / norm1).abs();
1044 max_relative_error = max_relative_error.max(relative_error);
1045 }
1046 }
1047 assert!(
1049 max_relative_error < 0.2,
1050 "Max relative error: {max_relative_error}"
1051 );
1052 }
1053
1054 #[test]
1055 fn test_frft_complex_input() {
1056 let n = 64;
1058 let signal_complex: Vec<Complex64> = (0..n)
1060 .map(|i| {
1061 let t = i as f64 / n as f64;
1062 Complex64::new((2.0 * PI * 5.0 * t).cos(), (2.0 * PI * 5.0 * t).sin())
1063 })
1064 .collect();
1065
1066 let result = frft_complex(&signal_complex, 0.5, None).unwrap();
1067
1068 assert_eq!(result.len(), n);
1070
1071 let result2 = frft_complex(&result, 0.5, None).unwrap();
1073 assert_eq!(result2.len(), n);
1074
1075 let result4 = frft_complex(&signal_complex, 4.0, None).unwrap();
1077 for i in 0..n {
1078 assert_relative_eq!(result4[i].re, signal_complex[i].re, epsilon = 1e-10);
1079 assert_relative_eq!(result4[i].im, signal_complex[i].im, epsilon = 1e-10);
1080 }
1081 }
1082}