Sequence1.hpp

Go to the documentation of this file.

// Copyright (c) 2011 Michael D. Adams
// All rights reserved.

// __START_OF_LICENSE__
// 
// Copyright (c) 2015 Michael D. Adams
// All rights reserved.
// 
// This file is part of the Signal Processing Library (SPL).
// 
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 3,
// or (at your option) any later version.
// 
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
// 
// You should have received a copy of the GNU General Public
// License along with this program; see the file LICENSE.  If not,
// see <http://www.gnu.org/licenses/>.
// 
// __END_OF_LICENSE__

#ifndef SPL_Sequence1_hpp
#define SPL_Sequence1_hpp

//

#define SPL_SEQUENCE1_DEBUG

#define SPL_SEQUENCE1_USE_NEW_CONV

#if defined(SPL_SEQUENCE1_DEBUG)
#define SPL_SEQUENCE1_INLINE
#else
#define SPL_SEQUENCE1_INLINE inline
#endif

// Header Files.

#include <SPL/config.hpp>
#include <iostream>
#include <vector>
#include <SPL/Array1.hpp>
#include <SPL/math.hpp>
#include <SPL/Sequence.hpp>

//

namespace SPL {

// One-Dimensional Sequence Template Class.

template <class T>
class Sequence1
{
public:

        // Types

        typedef T ElemType;

        typedef typename Array1<T>::ConstIterator ConstIterator;

        typedef typename Array1<T>::Iterator Iterator;

        // Constructors and destructor

        Sequence1();

        Sequence1(int startInd, int size);

        Sequence1(int startInd, int size, const T& value);

        template <class InputIterator>
        Sequence1(int startInd, int size, InputIterator data);

        Sequence1(const Sequence1& f);

        template <class OtherT>
        Sequence1(const Sequence1<OtherT>& f);

        Sequence1(const Array1<T>& data);

        Sequence1(int startInd, const Array1<T>& data);

        ~Sequence1();

        // Assignment operators

        Sequence1& operator=(const Sequence1& f);

        template <class OtherT>
        Sequence1& operator=(const Sequence1<OtherT>& f);

        // Compound assignment operators

        Sequence1& operator+=(const Sequence1& f);

        Sequence1& operator-=(const Sequence1& f);

        Sequence1& operator*=(const Sequence1& f);

        Sequence1& operator/=(const Sequence1& f);

        Sequence1& operator+=(const T& value);

        Sequence1& operator-=(const T& value);

        Sequence1& operator*=(const T& value);

        Sequence1& operator/=(const T& value);

        // Basic queries

        int getStartInd() const;

        int getEndInd() const;

        int getSize() const;

        bool isShared() const;

        // Accessors

        const T& operator()(int i) const;

        T& operator()(int i);

        // Iterators

        ConstIterator begin() const;

        Iterator begin();

        ConstIterator end() const;

        Iterator end();

        // Resizing

        // Get various statistics

        T min() const;

        T max() const;

        T sum() const;

        // Input/ouput

        // Miscellaneous

        Array1<T> getArray() const;

        void swapArray(Array1<T>& data);

        void fill(const T& value);

        Sequence1& translate(int delta);

private:

        // Do not look at the private members.
        // Doing so will make your eyes bleed.

        template <class> friend class Sequence1;

        template <class U>
        friend std::ostream& operator<<(std::ostream& out, const Sequence1<U>&);
        template <class U>
        friend Sequence1<U> translate(const Sequence1<U>&, int i);

        typedef Array1<T> Array;

        int startInd_;

        Array data_;

};

// Constructors and destructor

template <class T>
SPL_SEQUENCE1_INLINE Sequence1<T>::Sequence1() : startInd_(0), data_()
{
}

template <class T>
SPL_SEQUENCE1_INLINE Sequence1<T>::Sequence1(int startInd, int size) :
  startInd_(startInd), data_(size)
{
        assert(size >= 0);
}

template <class T>
SPL_SEQUENCE1_INLINE Sequence1<T>::Sequence1(int startInd, int size,
  const T& data) : startInd_(startInd), data_(size, data)
{
        assert(size >= 0);
}

template <class T>
template <class InputIterator>
SPL_SEQUENCE1_INLINE Sequence1<T>::Sequence1(int startInd, int size,
  InputIterator data) : startInd_(startInd), data_(size, data)
{
        assert(size >= 0);
}

template <class T>
SPL_SEQUENCE1_INLINE Sequence1<T>::Sequence1(const Array1<T>& f) :
  startInd_(0), data_(f)
{
}

template <class T>
SPL_SEQUENCE1_INLINE Sequence1<T>::Sequence1(int startInd,
  const Array1<T>& data) : startInd_(startInd), data_(data)
{
}

template <class T>
SPL_SEQUENCE1_INLINE Sequence1<T>::Sequence1(const Sequence1<T>& f) :
  startInd_(f.startInd_), data_(f.data_)
{
}

template <class T>
SPL_SEQUENCE1_INLINE Sequence1<T>::~Sequence1()
{
}

template <class T>
template <class OtherT>
SPL_SEQUENCE1_INLINE Sequence1<T>::Sequence1(const Sequence1<OtherT>& f) :
  startInd_(f.startInd_), data_(f.data_)
{
}

// Assignment operators

template <class T>
SPL_SEQUENCE1_INLINE Sequence1<T>& Sequence1<T>::operator=(
  const Sequence1<T>& f)
{
        startInd_ = f.startInd_;
        data_ = f.data_;
        return *this;
}

template <class T>
template <class OtherT>
SPL_SEQUENCE1_INLINE Sequence1<T>& Sequence1<T>::operator=(
  const Sequence1<OtherT>& f)
{
        startInd_ = f.startInd_;
        data_ = f.data_;
        return *this;
}

// Compound assignment operators

template <class T>
SPL_SEQUENCE1_INLINE Sequence1<T>& Sequence1<T>::operator+=(
  const Sequence1<T>& f)
{
        // Ensure that both sequences have the same domain.
        assert((getSize() == 0 && f.getSize() == 0) ||
          (getStartInd() == f.getStartInd() && getEndInd() == f.getEndInd()));

        data_ += f.data_;
        return *this;
}

template <class T>
SPL_SEQUENCE1_INLINE Sequence1<T>& Sequence1<T>::operator-=(
  const Sequence1<T>& f)
{
        // Ensure that both sequences have the same domain.
        assert((getSize() == 0 && f.getSize() == 0) ||
          (getStartInd() == f.getStartInd() && getEndInd() == f.getEndInd()));

        data_ -= f.data_;
        return *this;
}

template <class T>
SPL_SEQUENCE1_INLINE Sequence1<T>& Sequence1<T>::operator*=(
  const Sequence1<T>& f)
{
        // Ensure that both sequences have the same domain.
        assert((getSize() == 0 && f.getSize() == 0) ||
          (getStartInd() == f.getStartInd() && getEndInd() == f.getEndInd()));

        data_ *= f.data_;
        return *this;
}

template <class T>
SPL_SEQUENCE1_INLINE Sequence1<T>& Sequence1<T>::operator/=(
  const Sequence1<T>& f)
{
        // Ensure that both sequences have the same domain.
        assert((getSize() == 0 && f.getSize() == 0) ||
          (getStartInd() == f.getStartInd() && getEndInd() == f.getEndInd()));

        data_ /= f.data_;
        return *this;
}

template <class T>
SPL_SEQUENCE1_INLINE Sequence1<T>& Sequence1<T>::operator+=(const T& value)
{
        data_ += value;
        return *this;
}

template <class T>
SPL_SEQUENCE1_INLINE Sequence1<T>& Sequence1<T>::operator-=(const T& value)
{
        data_ -= value;
        return *this;
}

template <class T>
SPL_SEQUENCE1_INLINE Sequence1<T>& Sequence1<T>::operator*=(const T& value)
{
        data_ *= value;
        return *this;
}

template <class T>
SPL_SEQUENCE1_INLINE Sequence1<T>& Sequence1<T>::operator/=(const T& value)
{
        data_ /= value;
        return *this;
}

// Basic queries

template <class T>
SPL_SEQUENCE1_INLINE int Sequence1<T>::getStartInd() const
{
        return startInd_;
}

template <class T>
SPL_SEQUENCE1_INLINE int Sequence1<T>::getEndInd() const
{
        return startInd_ + getSize();
}

template <class T>
SPL_SEQUENCE1_INLINE int Sequence1<T>::getSize() const
{
        return data_.getSize();
}

template <class T>
SPL_SEQUENCE1_INLINE bool Sequence1<T>::isShared() const
{
        return data_.isShared();
}

// Accessors

template <class T>
SPL_SEQUENCE1_INLINE const T& Sequence1<T>::operator()(int i) const
{
        assert(i >= startInd_ && i < startInd_ + getSize());
        return data_(i - startInd_);
}

template <class T>
SPL_SEQUENCE1_INLINE T& Sequence1<T>::operator()(int i)
{
        assert(i >= startInd_ && i < startInd_ + getSize());
        return data_(i - startInd_);
}

// Iterators

template <class T>
SPL_SEQUENCE1_INLINE typename Sequence1<T>::ConstIterator Sequence1<T>::begin()
  const
{
        return data_.begin();
}

template <class T>
SPL_SEQUENCE1_INLINE typename Sequence1<T>::ConstIterator Sequence1<T>::end()
  const
{
        return data_.end();
}

template <class T>
SPL_SEQUENCE1_INLINE typename Sequence1<T>::Iterator Sequence1<T>::begin()
{
        return data_.begin();
}

template <class T>
SPL_SEQUENCE1_INLINE typename Sequence1<T>::Iterator Sequence1<T>::end()
{
        return data_.end();
}

// Resizing

// Get various statistics

template <class T>
SPL_SEQUENCE1_INLINE T Sequence1<T>::min() const
{
        assert(getSize() > 0);
        return data_.min();
}

template <class T>
SPL_SEQUENCE1_INLINE T Sequence1<T>::max() const
{
        assert(getSize() > 0);
        return data_.max();
}

template <class T>
SPL_SEQUENCE1_INLINE T Sequence1<T>::sum() const
{
        return data_.sum();
}


// Input/output

template <class T>
std::ostream& operator<<(std::ostream& out, const Sequence1<T>& f)
{
        out << f.getStartInd() << " ";
        out << f.data_;
        return out;
}

template <class T>
std::istream& operator>>(std::istream& in, Sequence1<T>& f)
{
        int startInd;
        int size;
        if (!(in >> startInd)) {
                return in;
        }
        if (!(in >> size)) {
                return in;
        }
        if (size < 0) {
                in.setstate(std::ios::failbit);
                return in;
        }
        typename std::vector<T> data;
        data.reserve(size);
        for (int i = 0; i < size; ++i) {
                T value;
                if (!(in >> value)) {
                        return in;
                }
                data.push_back(value);
        }
        f = Sequence1<T>(startInd, data.size(), data.begin());
        return in;
}

// Miscellaneous

template <class T>
SPL_SEQUENCE1_INLINE void Sequence1<T>::swapArray(Array1<T>& data)
{
        data_.swap(data);
}

template<class T>
SPL_SEQUENCE1_INLINE void Sequence1<T>::fill(const T& value)
{
        data_.fill(value);
}

template<class T>
SPL_SEQUENCE1_INLINE Array1<T> Sequence1<T>::getArray() const
{
        return data_;
}

template<class T>
SPL_SEQUENCE1_INLINE Sequence1<T>& Sequence1<T>::translate(int delta)
{
        startInd_ += delta;
        return *this;
}

//

template <class T>
SPL_SEQUENCE1_INLINE Sequence1<T> operator+(const Sequence1<T>& f,
  const Sequence1<T>& g)
{
        return Sequence1<T>(f) += g;
}

template <class T>
SPL_SEQUENCE1_INLINE Sequence1<T> operator-(const Sequence1<T>& f,
  const Sequence1<T>& g)
{
        return Sequence1<T>(f) -= g;
}

template <class T>
SPL_SEQUENCE1_INLINE Sequence1<T> operator*(const Sequence1<T>& f,
  const Sequence1<T>& g)
{
        return Sequence1<T>(f) *= g;
}

template <class T>
SPL_SEQUENCE1_INLINE Sequence1<T> operator/(const Sequence1<T>& f,
  const Sequence1<T>& g)
{
        return Sequence1<T>(f) /= g;
}

template <class T>
Sequence1<T> add(const Sequence1<T>& f, const Sequence1<T>& g)
{
        int startInd;
        int endInd;
        if (!f.getSize()) {
                startInd = g.getStartInd();
                endInd = g.getEndInd();
        } else if (!g.getSize()) {
                startInd = f.getStartInd();
                endInd = f.getEndInd();
        } else {
                startInd = std::min(f.getStartInd(), g.getStartInd());
                endInd = std::max(f.getStartInd() + f.getSize(), g.getStartInd() +
                  g.getSize());
        }
        Sequence1<T> result(startInd, endInd - startInd, T(0));
        int i0 = std::max(f.getStartInd(), startInd);
        int i1 = std::min(f.getStartInd() + f.getSize(), endInd);
        for (int i = i0; i < i1; ++i) {
                result(i) += f(i);
        }
        i0 = std::max(g.getStartInd(), startInd);
        i1 = std::min(g.getStartInd() + g.getSize(), endInd);
        for (int i = i0; i < i1; ++i) {
                result(i) += g(i);
        }
        return result;
}

template <class T>
SPL_SEQUENCE1_INLINE Sequence1<T> operator+(const T& a, const Sequence1<T>& f)
{
        return Sequence1<T>(f) += a;
}

template <class T>
SPL_SEQUENCE1_INLINE Sequence1<T> operator+(const Sequence1<T>& f, const T& a)
{
        return Sequence1<T>(f) += a;
}

template <class T>
SPL_SEQUENCE1_INLINE Sequence1<T> operator-(const Sequence1<T>& f, const T& a)
{
        return Sequence1<T>(f) -= a;
}

template <class T>
SPL_SEQUENCE1_INLINE Sequence1<T> operator*(const T& a, const Sequence1<T>& f)
{
        return Sequence1<T>(f) *= a;
}

template <class T>
SPL_SEQUENCE1_INLINE Sequence1<T> operator*(const Sequence1<T>& f, const T& a)
{
        return Sequence1<T>(f) *= a;
}

template <class T>
SPL_SEQUENCE1_INLINE Sequence1<T> operator/(const Sequence1<T>& f, const T& a)
{
        return Sequence1<T>(f) /= a;
}

template <class T>
SPL_SEQUENCE1_INLINE bool operator==(const Sequence1<T>& f,
  const Sequence1<T>& g)
{
        bool result;
        if (!f.getSize() && !g.getSize()) {
                result = true;
        } else if (f.getStartInd() != g.getStartInd() ||
          f.getEndInd() != g.getEndInd()) {
                result = false;
        } else {
                result = std::equal(f.begin(), f.end(), g.begin());
        }
        return result;
}

template <class T>
SPL_SEQUENCE1_INLINE bool operator!=(const Sequence1<T>& f,
  const Sequence1<T>& g)
{
        return !(f == g);
}

template <class T>
SPL_SEQUENCE1_INLINE bool approxEqual(const Sequence1<T>& f,
  const Sequence1<T>& g, T threshold = 1e-9)
{
        bool result;
        if (!f.getSize() && !g.getSize()) {
                // Empty sequence.
                result = true;
        } else if (f.getStartInd() != g.getStartInd() ||
          f.getEndInd() != g.getEndInd()) {
                // Domains do not match.
                result = false;
        } else {
                result = true;
                typename Sequence1<T>::ConstIterator src = f.begin();
                typename Sequence1<T>::ConstIterator dst = g.begin();
                while (src != f.end()) {
                        if (absVal((*src) - (*dst)) > threshold) {
                                result = false;
                                break;
                        }
                        ++src;
                        ++dst;
                }
        }
        return result;
}

//

template <class T>
Sequence1<T> subsequence(const Sequence1<T>& f, int startInd, int size)
{
        assert(startInd >= f.getStartInd() && startInd + size <= f.getEndInd());
        Sequence1<T> result(startInd, size);
        for (int i = 0; i < size; ++i) {
                result(startInd + i) = f(startInd + i);
        }
        return result;
}

template <class T>
SPL_SEQUENCE1_INLINE Sequence1<T> translate(const Sequence1<T>& f, int delta)
{
        return Sequence1<T>(f).translate(delta);
}

// Convolution


template<class T>
SPL_SEQUENCE1_INLINE T getExtSeqValue(const Sequence1<T>& f, int i, int mode)
{
        assert(mode == ConvolveMode::full ||
          mode == ConvolveMode::sameDomainConstExt ||
          mode == ConvolveMode::sameDomainZeroExt ||
          mode == ConvolveMode::sameDomainPerExt ||
          mode == ConvolveMode::sameDomainSymExt0);
        T value;
        switch (mode) {
        case ConvolveMode::full:
        case ConvolveMode::sameDomainZeroExt:
                if (i >= f.getStartInd() && i < f.getEndInd()) {
                        value = f(i);
                } else {
                        value = T(0);
                }
                break;
        case ConvolveMode::sameDomainConstExt:
                if (i < f.getStartInd()) {
                        i = f.getStartInd();
                } else if (i >= f.getEndInd()) {
                        i = f.getEndInd() - 1;
                }
                assert(i >= f.getStartInd() && i < f.getEndInd());
                value = f(i);
                break;
        case ConvolveMode::sameDomainPerExt:
                if (i < f.getStartInd() || i >= f.getEndInd()) {
                        i = f.getStartInd() + mod(i - f.getStartInd(), f.getSize());
                }
                assert(i >= f.getStartInd() && i < f.getEndInd());
                value = f(i);
                break;
        case ConvolveMode::sameDomainSymExt0:
                if (i < f.getStartInd() || i >= f.getEndInd()) {
                        int tmp = mod(i - f.getStartInd(), 2 * f.getSize() - 2);
                        i = f.getStartInd() + std::min(tmp, 2 * f.getSize() - 2 - tmp);
                }
                assert(i >= f.getStartInd() && i < f.getEndInd());
                value = f(i);
                break;
        default:
                assert(0);
                abort();
                break;
        }
        return value;
}

template <class T>
Sequence1<T> convolveFull(const Sequence1<T>& f, const Sequence1<T>& g)
{
        Sequence1<T> seq(f);
        Sequence1<T> filt(g);
        Sequence1<T> result(f.getStartInd() + g.getStartInd(),
          f.getSize() + g.getSize() - 1);
        convolveHelper(seq.begin(), seq.end(), filt.begin(), filt.end(),
          result.begin(), T());
        return result;
}

template <class T>
Sequence1<T> convolveSameDomainZeroExt(const Sequence1<T>& f,
  const Sequence1<T>& g)
{
        Sequence1<T> result(f.getStartInd(), f.getSize());
        int skip;
        int shift;
        if (g.getStartInd() >= 0) {
                shift = g.getStartInd();
                skip = 0;
        } else {
                shift = 0;
                skip = -g.getStartInd();
        }
        convolveHelper2(f.begin(), f.end(), g.begin(), g.end(),
          result.begin() + shift, skip, result.getSize() - shift, T());
        return result;
}

template <class T>
Sequence1<T> convolve(const Sequence1<T>& f, const Sequence1<T>& g,
  int mode = ConvolveMode::full)
{
#if defined(SPL_SEQUENCE1_USE_NEW_CONV)
        switch (mode) {
        case ConvolveMode::sameDomainZeroExt:
                return convolveSameDomainZeroExt(f, g);
                break;
        case ConvolveMode::full:
                return convolveFull(f, g);
                break;
        default:
                assert(0);
                abort();
                break;
        }
#endif

        Sequence1<T> result;

        switch (mode) {
        case ConvolveMode::full:
                result = Sequence1<T>(f.getStartInd() + g.getStartInd(),
                  f.getSize() + g.getSize() - 1);
                break;
        case ConvolveMode::sameDomainZeroExt:
        case ConvolveMode::sameDomainConstExt:
        case ConvolveMode::sameDomainPerExt:
        case ConvolveMode::sameDomainSymExt0:
                result = Sequence1<T>(f.getStartInd(), f.getSize());
                break;
        default:
                assert(0);
                abort();
                break;
        }

        typename Sequence1<T>::Iterator resultIter = result.begin();
        for (int n = result.getStartInd(); n < result.getEndInd(); ++n) {
                T sum = T(0);
                for (int k = g.getStartInd(); k < g.getEndInd(); ++k) {
                        sum += getExtSeqValue(f, n - k, mode) * g(k);
                }
                *resultIter = sum;
                ++resultIter;
        }
        assert(resultIter == result.end());
        return result;
}

// Downsampling and upsampling

template <class T>
Sequence1<T> downsample(const Sequence1<T>& f, int factor)
{
        assert(factor >= 1);
        int startInd = ceilDiv(f.getStartInd(), factor) * factor;
        int n = ceilDiv(f.getEndInd() - startInd, factor);
        Sequence1<T> g(roundTowardZeroDiv(startInd, factor), n);
        typename Sequence1<T>::ConstIterator srcIter = f.begin() +
          (startInd - f.getStartInd());
        typename Sequence1<T>::Iterator dstIter = g.begin();
        while (--n >= 0) {
                *dstIter = *srcIter;
                ++dstIter;
                srcIter += factor;
        }
        return g;
}

template <class T>
Sequence1<T> upsample(const Sequence1<T>& f, int factor, int pad = 0)
{
        assert(factor >= 1);
        assert(pad >= 0 && pad < factor);
        int n = f.getSize() + std::max((f.getSize() - 1), 0) * (factor - 1) + pad;
        int m = f.getSize();
        Sequence1<T> g(factor * f.getStartInd(), n);
        typename Sequence1<T>::ConstIterator srcIter = f.begin();
        typename Sequence1<T>::Iterator dstIter = g.begin();
        while (m > 1) {
                *dstIter = *srcIter;
                ++dstIter;
                ++srcIter;
                for (int i = factor - 1; i > 0; --i) {
                        *dstIter = T(0);
                        ++dstIter;
                }
                --m;
        }
        if (m > 0) {
                *dstIter = *srcIter;
                ++dstIter;
        }
        for (int i = pad; i > 0; --i) {
                *dstIter = T(0);
                ++dstIter;
        }
        return g;
}

// Polyphase Decompositions.


inline int getCosetOffset(int type, int numPhases, int index)
{
        int result;
        switch (type) {
        case 1:
                result = -index;
                break;
        case 2:
                result = index - (numPhases - 1);
                break;
        case 3:
                result = index;
                break;
        case 4:
                result = (numPhases - 1) - index;
                break;
        default:
                assert(0);
                abort();
                break;
        }
        return result;
}

template <class T>
Array1<Sequence1<T> > polyphaseSplit(const Sequence1<T>& seq, int type,
  int numPhases)
{
        assert(numPhases >= 2);
        Array1<Sequence1<T> > result(numPhases);
        for (int i = 0; i < numPhases; ++i) {
                int offset = getCosetOffset(type, numPhases, i);
                result(i) = downsample(translate(seq, offset), numPhases);
        }
        return result;
}

template <class T>
Sequence1<T> polyphaseJoin(const Array1<Sequence1<T> >& comps, int type)
{
        const int numPhases = comps.getSize();
        assert(numPhases >= 2);
        Sequence1<T> seq;
        for (int i = 0; i < numPhases; ++i) {
                int offset = getCosetOffset(type, numPhases, i);
                seq = add(seq, translate(upsample(comps(i), numPhases), -offset));
        }
        return seq;
}

// Types.

typedef Sequence1<double> RealSequence1;

typedef Sequence1<int> IntSequence1;


}

#endif