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

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

/// 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::*;
    /// # fn main() {
    /// let mut vector = vec!(1.0, 2.0, 3.0, 4.0).to_complex_time_vec();
    /// vector.multiply_complex_exponential(2.0, 3.0);
    /// let actual = &vector[..];
    /// let expected = &[-1.2722325, -1.838865, 4.6866837, -1.7421241];
    /// assert_eq!(actual.len(), expected.len());
    /// for i in 0..actual.len() {
    ///        assert!(f64::abs(actual[i] - expected[i]) < 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::*;
    /// # fn main() {
    /// let mut vector = vec!(1.0, 2.0, 3.0, 4.0).to_complex_time_vec();
    /// vector.conj();
    /// assert_eq!([1.0, -2.0, 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_.valid_len = 0;
        }
    }
}

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 mut array = self.data.to_slice_mut();
        Chunk::execute_with_range(Complexity::Small,
                                  &self.multicore_settings,
                                  &mut array[0..data_length],
                                  T::Reg::len(),
                                  (a, b),
                                  move |array, range, args| {
            let two = T::one() + T::one();
            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).unwrap() / two));
            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 = T::Reg::from_complex(Complex::<T>::new(T::one(), -T::one()));
        self.simd_complex_operation(|x, y| x * y, |x, _| x.conj(), factor, Complexity::Small)
    }
}