1use crate::signal::{
2 FreqSignal,
3 data::{ComplexData, RealData},
4 math::{gain_to_db, wrap_phase_2pi},
5};
6
7const NEG_INF_DB: f64 = -100.0;
8
9impl ComplexData {
10 pub fn to_magnitude(&self) -> RealData {
11 RealData::from_axis(
12 self.x_axis().clone(),
13 self.y_data().mapv(|value| value.norm()),
14 )
15 .expect("x/y data originated from ComplexData")
16 }
17
18 pub fn to_magnitude_db(&self) -> RealData {
19 RealData::from_axis(
20 self.x_axis().clone(),
21 self.y_data()
22 .mapv(|value| gain_to_db(value.norm()).max(NEG_INF_DB)),
23 )
24 .expect("x/y data originated from ComplexData")
25 }
26
27 pub fn to_phase(&self) -> RealData {
28 RealData::from_axis(
29 self.x_axis().clone(),
30 self.y_data().mapv(|value| value.arg()),
31 )
32 .expect("x/y data originated from ComplexData")
33 }
34
35 pub fn to_phase_unwrapped(&self) -> RealData {
36 let mut phase = self.y_data().mapv(|value| value.arg());
37
38 for mut channel in phase.outer_iter_mut() {
39 unwrap_phase_in_place(channel.as_slice_mut().expect("channel view is contiguous"));
40 }
41
42 RealData::from_axis(self.x_axis().clone(), phase)
43 .expect("x/y data originated from ComplexData")
44 }
45
46 pub fn to_phase_2pi(&self) -> RealData {
47 RealData::from_axis(
48 self.x_axis().clone(),
49 self.y_data().mapv(|value| wrap_phase_2pi(value.arg())),
50 )
51 .expect("x/y data originated from ComplexData")
52 }
53
54 pub fn to_phase_degrees(&self) -> RealData {
55 RealData::from_axis(
56 self.x_axis().clone(),
57 self.y_data().mapv(|value| value.arg().to_degrees()),
58 )
59 .expect("x/y data originated from ComplexData")
60 }
61
62 pub fn to_phase_degrees_unwrapped(&self) -> RealData {
63 let mut phase = self.y_data().mapv(|value| value.arg());
64
65 for mut channel in phase.outer_iter_mut() {
66 unwrap_phase_in_place(channel.as_slice_mut().expect("channel view is contiguous"));
67 }
68
69 RealData::from_axis(self.x_axis().clone(), phase.mapv(f64::to_degrees))
70 .expect("x/y data originated from ComplexData")
71 }
72}
73
74impl FreqSignal {
75 pub fn to_magnitude(&self) -> RealData {
76 self.data().to_magnitude()
77 }
78
79 pub fn to_magnitude_db(&self) -> RealData {
80 self.data().to_magnitude_db()
81 }
82
83 pub fn to_phase(&self) -> RealData {
84 self.data().to_phase()
85 }
86
87 pub fn to_phase_unwrapped(&self) -> RealData {
88 self.data().to_phase_unwrapped()
89 }
90
91 pub fn to_phase_2pi(&self) -> RealData {
92 self.data().to_phase_2pi()
93 }
94
95 pub fn to_phase_degrees(&self) -> RealData {
96 self.data().to_phase_degrees()
97 }
98
99 pub fn to_phase_degrees_unwrapped(&self) -> RealData {
100 self.data().to_phase_degrees_unwrapped()
101 }
102}
103
104fn unwrap_phase_in_place(phase: &mut [f64]) {
105 if phase.len() < 2 {
106 return;
107 }
108
109 let two_pi = 2.0 * std::f64::consts::PI;
110 let mut offset = 0.0;
111 let mut previous = phase[0];
112
113 for current in phase.iter_mut().skip(1) {
114 let current_value = *current;
115 let delta = current_value - previous;
116 if delta > std::f64::consts::PI {
117 offset -= two_pi;
118 } else if delta < -std::f64::consts::PI {
119 offset += two_pi;
120 }
121 *current = current_value + offset;
122 previous = current_value;
123 }
124}
125
126#[cfg(test)]
127mod tests {
128 use std::f64::consts::{FRAC_PI_2, FRAC_PI_4, PI};
129
130 use approx::assert_abs_diff_eq;
131 use ndarray::{arr1, arr2};
132 use num::complex::Complex64;
133
134 use crate::signal::data::ComplexData;
135
136 #[test]
137 fn complex_data_conversions() {
138 let data = ComplexData::new(
139 arr1(&[10.0, 20.0, 30.0]),
140 arr2(&[[
141 Complex64::new(3.0, 4.0),
142 Complex64::new(0.0, 1.0),
143 Complex64::new(-1.0, 0.0),
144 ]]),
145 )
146 .unwrap();
147
148 let magnitude = data.to_magnitude();
149 assert_abs_diff_eq!(magnitude.channel(0)[0], 5.0, epsilon = 1e-12);
150 assert_abs_diff_eq!(magnitude.channel(0)[1], 1.0, epsilon = 1e-12);
151
152 let magnitude_db = data.to_magnitude_db();
153 assert_abs_diff_eq!(
154 magnitude_db.channel(0)[0],
155 20.0 * 5.0_f64.log10(),
156 epsilon = 1e-12
157 );
158 assert_abs_diff_eq!(magnitude_db.channel(0)[1], 0.0, epsilon = 1e-12);
159
160 let phase = data.to_phase();
161 assert_abs_diff_eq!(phase.channel(0)[0], (4.0_f64).atan2(3.0), epsilon = 1e-12);
162 assert_abs_diff_eq!(phase.channel(0)[1], FRAC_PI_2, epsilon = 1e-12);
163 assert_abs_diff_eq!(phase.channel(0)[2], PI, epsilon = 1e-12);
164 }
165
166 #[test]
167 fn complex_data_phase_quadrants() {
168 let data = ComplexData::new(
169 arr1(&[0.0, 1.0, 2.0, 3.0]),
170 arr2(&[[
171 Complex64::new(1.0, 0.0),
172 Complex64::new(0.0, 1.0),
173 Complex64::new(1.0, 1.0),
174 Complex64::new(1.0, -1.0),
175 ]]),
176 )
177 .unwrap();
178
179 let phase = data.to_phase();
180 assert_abs_diff_eq!(phase.channel(0)[0], 0.0, epsilon = 1e-12);
181 assert_abs_diff_eq!(phase.channel(0)[1], FRAC_PI_2, epsilon = 1e-12);
182 assert_abs_diff_eq!(phase.channel(0)[2], FRAC_PI_4, epsilon = 1e-12);
183 assert_abs_diff_eq!(phase.channel(0)[3], -FRAC_PI_4, epsilon = 1e-12);
184 }
185
186 #[test]
187 fn complex_data_phase_variants() {
188 let data = ComplexData::new(
189 arr1(&[0.0, 1.0, 2.0]),
190 arr2(&[[
191 Complex64::from_polar(1.0, 3.0 * PI / 4.0),
192 Complex64::from_polar(1.0, -3.0 * PI / 4.0),
193 Complex64::from_polar(1.0, -PI / 2.0),
194 ]]),
195 )
196 .unwrap();
197
198 let unwrapped = data.to_phase_unwrapped();
199 let wrapped_2pi = data.to_phase_2pi();
200 let degrees = data.to_phase_degrees();
201 let degrees_unwrapped = data.to_phase_degrees_unwrapped();
202
203 assert_abs_diff_eq!(unwrapped.channel(0)[0], 3.0 * PI / 4.0, epsilon = 1e-12);
204 assert_abs_diff_eq!(unwrapped.channel(0)[1], 5.0 * PI / 4.0, epsilon = 1e-12);
205 assert_abs_diff_eq!(unwrapped.channel(0)[2], 3.0 * PI / 2.0, epsilon = 1e-12);
206
207 assert_abs_diff_eq!(wrapped_2pi.channel(0)[1], 5.0 * PI / 4.0, epsilon = 1e-12);
208 assert_abs_diff_eq!(wrapped_2pi.channel(0)[2], 3.0 * PI / 2.0, epsilon = 1e-12);
209
210 assert_abs_diff_eq!(degrees.channel(0)[0], 135.0, epsilon = 1e-12);
211 assert_abs_diff_eq!(degrees.channel(0)[1], -135.0, epsilon = 1e-12);
212 assert_abs_diff_eq!(degrees_unwrapped.channel(0)[1], 225.0, epsilon = 1e-12);
213 assert_abs_diff_eq!(degrees_unwrapped.channel(0)[2], 270.0, epsilon = 1e-12);
214 }
215}