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

use super::super::{ComplexNumberSpace, Domain, DspVec, MetaData, ToSliceMut, Vector};
use crate::array_to_complex_mut;
use crate::multicore_support::*;
use crate::numbers::*;
use crate::simd_extensions::*;

/// Operations on complex types.
///
/// # Failures
///
/// If one of the methods is called on real 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 complex.
pub trait ComplexOps<T>
where
    T: RealNumber,
{
    /// Multiplies each vector element with `exp(j*(a*idx*self.delta() + b))`
    /// where `a` and `b` are arguments and `idx` is the index of the data points
    /// in the vector ranging from `0 to self.points() - 1`. `j` is the imaginary number and
    /// `exp` the exponential function.
    ///
    /// This method can be used to perform a frequency shift in time domain.
    ///
    /// # Example
    ///
    /// ```
    /// # use std::f64;
    /// # extern crate num_complex;
    /// # extern crate basic_dsp_vector;
    /// use basic_dsp_vector::*;
    /// # use num_complex::Complex;
    /// # fn main() {
    /// let mut vector = vec!(Complex::new(1.0, 2.0), Complex::new(3.0, 4.0)).to_complex_time_vec();
    /// vector.multiply_complex_exponential(2.0, 3.0);
    /// let actual = &vector[..];
    /// let expected = &[Complex::new(-1.2722325, -1.838865), Complex::new(4.6866837, -1.7421241)];
    /// assert_eq!(actual.len(), expected.len());
    /// for i in 0..actual.len() {
    ///        assert!((actual[i] - expected[i]).norm() < 1e-4);
    /// }
    /// # }
    /// ```
    fn multiply_complex_exponential(&mut self, a: T, b: T);

    /// Calculates the complex conjugate of the vector.
    /// # Example
    ///
    /// ```
    /// # extern crate num_complex;
    /// # extern crate basic_dsp_vector;
    /// use basic_dsp_vector::*;
    /// # use num_complex::Complex;
    /// # fn main() {
    /// let mut vector = vec!(Complex::new(1.0, 2.0), Complex::new(3.0, 4.0)).to_complex_time_vec();
    /// vector.conj();
    /// assert_eq!([Complex::new(1.0, -2.0), Complex::new(3.0, -4.0)], vector[..]);
    /// # }
    /// ```
    fn conj(&mut self);
}

macro_rules! assert_complex {
    ($self_: ident) => {
        if !$self_.is_complex() {
            $self_.number_space.to_real();
            $self_.mark_vector_as_invalid();
        }
    };
}

impl<S, T, N, D> ComplexOps<T> for DspVec<S, T, N, D>
where
    S: ToSliceMut<T>,
    T: RealNumber,
    N: ComplexNumberSpace,
    D: Domain,
{
    fn multiply_complex_exponential(&mut self, a: T, b: T) {
        assert_complex!(self);
        let a = a * self.delta();
        let b = b * self.delta();
        let data_length = self.len();
        let array = self.data.to_slice_mut();
        Chunk::execute_with_range(
            Complexity::Small,
            &self.multicore_settings,
            &mut array[0..data_length],
            2,
            (a, b),
            move |array, range, args| {
                let (a, b) = args;
                let mut exponential = Complex::<T>::from_polar(T::one(), b)
                    * Complex::<T>::from_polar(T::one(), a * T::from(range.start / 2).unwrap());
                let increment = Complex::<T>::from_polar(T::one(), a);
                let array = array_to_complex_mut(array);
                for complex in array {
                    *complex = (*complex) * exponential;
                    exponential = exponential * increment;
                }
            },
        );
    }

    fn conj(&mut self) {
        assert_complex!(self);
        let factor = Complex::<T>::new(T::one(), -T::one());
        sel_reg!(self.simd_complex_operationf::<T>(
            |x, y| x * y,
            |x, _| x.conj(),
            factor,
            Complexity::Small
        ))
    }
}