1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
use super::super::{Domain, DspVec, MetaData, RealNumberSpace, ToSliceMut, Vector};
use crate::multicore_support::*;
use crate::numbers::*;
use crate::simd_extensions::*;

/// Operations on real types.
///
/// # Failures
///
/// If one of the methods is called on complex data then `self.len()` will be set to `0`.
/// To avoid this it's recommended to use the `to_real_time_vec`, `to_real_freq_vec`
/// `to_complex_time_vec` and `to_complex_freq_vec` constructor methods since
/// the resulting types will already check at compile time (using the type system)
/// that the data is real.
pub trait RealOps {
    /// Gets the absolute value of all vector elements.
    /// # Example
    ///
    /// ```
    /// use basic_dsp_vector::*;
    /// let mut vector = vec!(1.0, -2.0).to_real_time_vec();
    /// vector.abs();
    /// assert_eq!([1.0, 2.0], vector[..]);
    /// ```
    fn abs(&mut self);
}

/// Operations on real types.
///
/// # Failures
///
/// If one of the methods is called on complex data then `self.len()` will be set to `0`.
/// To avoid this it's recommended to use the `to_real_time_vec`, `to_real_freq_vec`
/// `to_complex_time_vec` and `to_complex_freq_vec` constructor methods since
/// the resulting types will already check at compile time (using the type system)
/// that the data is real.
pub trait ModuloOps<T>
where
    T: RealNumber,
{
    /// Each value in the vector is dividable by the divisor and the remainder
    /// is stored in the resulting
    /// vector. This the same a modulo operation or to phase wrapping.
    ///
    /// # Example
    ///
    /// ```
    /// use basic_dsp_vector::*;
    /// let mut vector = vec!(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0).to_real_time_vec();
    /// vector.wrap(4.0);
    /// assert_eq!([1.0, 2.0, 3.0, 0.0, 1.0, 2.0, 3.0, 0.0], vector[..]);
    /// ```
    fn wrap(&mut self, divisor: T);

    /// This function corrects the jumps in the given vector which occur due
    /// to wrap or modulo operations.
    /// This will undo a wrap operation only if the deltas are smaller than half the divisor.
    ///
    /// # Example
    ///
    /// ```
    /// use basic_dsp_vector::*;
    /// let mut vector = vec!(1.0, 2.0, 3.0, 0.0, 1.0, 2.0, 3.0, 0.).to_real_time_vec();
    /// vector.unwrap(4.0);
    /// assert_eq!([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0], vector[..]);
    /// ```
    fn unwrap(&mut self, divisor: T);
}

/// Recommended to be only used with the CPU feature flags `sse` or `avx`.
///
/// This trait provides alternative implementations for some standard functions which
/// are less accurate but perform faster. Those approximations are written for SSE2 or AVX2 processors.
/// If a processor supports neither of those instruction sets, then the standard functions will be used.
///
/// Information on the error of the approximation and their performance are rough numbers.
/// A detailed table can be obtained by running the `approx_accuracy` example.
///
/// # Failures
///
/// If one of the methods is called on complex data then `self.len()` will be set to `0`.
/// To avoid this it's recommended to use the `to_real_time_vec`, `to_real_freq_vec`
/// `to_complex_time_vec` and `to_complex_freq_vec` constructor methods since
/// the resulting types will already check at compile time (using the type system)
/// that the data is real.
pub trait ApproximatedOps<T>
where
    T: RealNumber,
{
    /// Computes the principal value approximation of natural logarithm of every element in the vector.
    ///
    /// Error should be below `1%` as long as the values in the vector are larger than `1`.
    /// Single core performance should be about `5x` as fast.
    ///
    /// # Example
    ///
    /// ```no_run
    /// # use std::f64;
    /// use basic_dsp_vector::*;
    /// let mut vector = vec!(2.718281828459045, 7.389056, 20.085537).to_real_time_vec();
    /// vector.ln_approx();
    /// let actual = &vector[0..];
    /// let expected = &[1.0, 2.0, 3.0];
    /// assert_eq!(actual.len(), expected.len());
    /// for i in 0..actual.len() {
    ///        assert!(f64::abs(actual[i] - expected[i]) < 1e-2);
    /// }
    /// ```
    fn ln_approx(&mut self);

    /// Calculates the natural exponential approximation for every vector element.
    ///
    /// Error should be less than `1%`` as long as the values in the vector are small
    /// (e.g. in the range between -10 and 10).
    /// Single core performance should be about `50%` faster.
    ///
    /// # Example
    ///
    /// ```no_run
    /// # use std::f64;
    /// use basic_dsp_vector::*;
    /// let mut vector = vec!(1.0, 2.0, 3.0).to_real_time_vec();
    /// vector.exp_approx();
    /// let actual = &vector[0..];
    /// let expected = &[2.718281828459045, 7.389056, 20.085537];
    /// assert_eq!(actual.len(), expected.len());
    /// for i in 0..actual.len() {
    ///        assert!(f64::abs(actual[i] - expected[i]) < 1e-4);
    /// }
    /// ```
    fn exp_approx(&mut self);

    /// Calculates the sine approximation of each element in radians.
    ///
    /// Error should be below `1E-6`.
    /// Single core performance should be about `2x` as fast.
    ///
    /// # Example
    ///
    /// ```
    /// use std::f32;
    /// use basic_dsp_vector::*;
    /// let mut vector = vec!(f32::consts::PI/2.0, -f32::consts::PI/2.0).to_real_time_vec();
    /// vector.sin_approx();
    /// assert_eq!([1.0, -1.0], vector[..]);
    /// ```
    fn sin_approx(&mut self);

    /// Calculates the cosine approximation of each element in radians
    ///
    /// Error should be below `1E-6`.
    /// Single core performance should be about `2x` as fast.
    ///
    /// # Example
    ///
    /// ```
    /// use std::f32;
    /// use basic_dsp_vector::*;
    /// let mut vector = vec!(2.0 * f32::consts::PI, f32::consts::PI).to_real_time_vec();
    /// vector.cos_approx();
    /// assert_eq!([1.0, -1.0], vector[..]);
    /// ```
    fn cos_approx(&mut self);

    /// Calculates the approximated logarithm to the given base for every vector element.
    ///
    /// Error should be below `1%` as long as the values in the vector are larger than `1`.
    /// Single core performance should be about `5x` as fast.
    ///
    /// # Example
    ///
    /// ```no_run
    /// # use std::f64;
    /// use basic_dsp_vector::*;
    /// let mut vector = vec!(10.0, 100.0, 1000.0).to_real_time_vec();
    /// vector.log_approx(10.0);
    /// let actual = &vector[0..];
    /// let expected = &[1.0, 2.0, 3.0];
    /// assert_eq!(actual.len(), expected.len());
    /// for i in 0..actual.len() {
    ///        assert!(f64::abs(actual[i] - expected[i]) < 1e-4);
    /// }
    /// ```
    fn log_approx(&mut self, base: T);

    /// Calculates the approximated exponential to the given base for every vector element.
    ///
    /// Error should be less than `1%`` as long as the values in the vector are small
    /// (e.g. in the range between -10 and 10).
    /// Single core performance should be about `5x` as fast.
    ///
    /// # Example
    ///
    /// ```no_run
    /// use basic_dsp_vector::*;
    /// use std::f32;
    /// let vector: Vec<f32> = vec!(1.0, 2.0, 3.0);
    /// let mut vector = vector.to_real_time_vec();
    /// vector.expf_approx(10.0);
    /// assert!((vector[0] - 10.0).abs() < 1e-3);
    /// assert!((vector[1] - 100.0).abs() < 1e-3);
    /// assert!((vector[2] - 1000.0).abs() < 1e-3);
    /// ```
    fn expf_approx(&mut self, base: T);

    /// Raises every vector element to approximately a floating point power.
    ///
    /// Error should be less than `1%`` as long as the values in the vector are really small
    /// (e.g. in the range between -0.1 and 0.1).
    /// Single core performance should be about `5x` as fast.
    ///
    /// # Example
    ///
    /// ```no_run
    /// use basic_dsp_vector::*;
    /// use std::f32;
    /// let vector: Vec<f32> = vec!(1.0, 2.0, 3.0);
    /// let mut vector = vector.to_real_time_vec();
    /// vector.powf_approx(3.0);
    /// assert!((vector[0] - 1.0).abs() < 1e-3);
    /// assert!((vector[1] - 8.0).abs() < 1e-3);
    /// assert!((vector[2] - 27.0).abs() < 1e-3);
    /// ```
    fn powf_approx(&mut self, exponent: T);
}

macro_rules! assert_real {
    ($self_: ident) => {
        if $self_.is_complex() {
            $self_.mark_vector_as_invalid();
            return;
        }
    };
}

impl<S, T, N, D> RealOps for DspVec<S, T, N, D>
where
    S: ToSliceMut<T>,
    T: RealNumber,
    N: RealNumberSpace,
    D: Domain,
{
    fn abs(&mut self) {
        assert_real!(self);
        sel_reg!(self.simd_real_operation::<T>(
            |x, _arg| (x * x).sqrt(),
            |x, _arg| x.abs(),
            (),
            Complexity::Small
        ));
    }
}

impl<S, T, N, D> ModuloOps<T> for DspVec<S, T, N, D>
where
    S: ToSliceMut<T>,
    T: RealNumber,
    N: RealNumberSpace,
    D: Domain,
{
    fn wrap(&mut self, divisor: T) {
        assert_real!(self);
        self.pure_real_operation(|x, y| x % y, divisor, Complexity::Small);
    }

    fn unwrap(&mut self, divisor: T) {
        assert_real!(self);
        let data_length = self.len();
        let data = self.data.to_slice_mut();
        let mut i = 0;
        let mut j = 1;
        let half = divisor / T::from(2.0).unwrap();
        while j < data_length {
            let mut diff = data[j] - data[i];
            if diff > half {
                diff = diff % divisor;
                diff = diff - divisor;
                data[j] = data[i] + diff;
            } else if diff < -half {
                diff = diff % divisor;
                diff = diff + divisor;
                data[j] = data[i] + diff;
            }

            i += 1;
            j += 1;
        }
    }
}

/// Complexity of all approximation functions. Small, because with the SIMD optimizations
/// threading becomes less beneficial.
const APPROX_COMPLEXITY: Complexity = Complexity::Small;

impl<S, T, N, D> ApproximatedOps<T> for DspVec<S, T, N, D>
where
    S: ToSliceMut<T>,
    T: RealNumber,
    N: RealNumberSpace,
    D: Domain,
{
    fn ln_approx(&mut self) {
        assert_real!(self);
        sel_reg!(self.simd_real_operation::<T>(
            |x, _arg| x.ln_approx(),
            |x, _arg| x.ln(),
            (),
            APPROX_COMPLEXITY
        ));
    }

    fn exp_approx(&mut self) {
        assert_real!(self);
        sel_reg!(self.simd_real_operation::<T>(
            |x, _arg| x.exp_approx(),
            |x, _arg| x.exp(),
            (),
            APPROX_COMPLEXITY
        ));
    }

    fn sin_approx(&mut self) {
        assert_real!(self);
        sel_reg!(self.simd_real_operation::<T>(
            |x, _arg| x.sin_approx(),
            |x, _arg| x.sin(),
            (),
            APPROX_COMPLEXITY
        ));
    }

    fn cos_approx(&mut self) {
        assert_real!(self);
        sel_reg!(self.simd_real_operation::<T>(
            |x, _arg| x.cos_approx(),
            |x, _arg| x.cos(),
            (),
            APPROX_COMPLEXITY
        ));
    }

    fn log_approx(&mut self, base: T) {
        assert_real!(self);
        sel_reg!(self.simd_real_operationf::<T>(
            |x, b| x.ln_approx() / b,
            |x, b| x.ln() / b.extract(0),
            base.ln(),
            APPROX_COMPLEXITY
        ));
    }

    fn expf_approx(&mut self, base: T) {
        assert_real!(self);
        // Transform base with:
        // x^y = e^(ln(x)*y)
        sel_reg!(self.simd_real_operationf::<T>(
            |y, x| (x * y).exp_approx(),
            |y, x| (x.extract(0) * y).exp(),
            base.ln(),
            APPROX_COMPLEXITY
        ));
    }

    fn powf_approx(&mut self, exponent: T) {
        assert_real!(self);
        // Transform base with the same equation as for `expf_approx`
        sel_reg!(self.simd_real_operationf::<T>(
            |x, y| (x.ln_approx() * y).exp_approx(),
            |x, y| (x.ln() * y.extract(0)).exp(),
            exponent,
            APPROX_COMPLEXITY
        ));
    }
}