Kernels.icc

#include <cmath>
#include <complex>
#include <type_traits>
 
#include <isce3/except/Error.h>
#include <isce3/math/Sinc.h>
 
namespace isce3 { namespace core {
 
/*
 * Bartlett
 */
 
// call
template<typename T>
T BartlettKernel<T>::operator()(double t) const
{
    double t2 = fabs(t / this->_halfwidth);
    if (t2 > 1.0) {
        return T(0.0);
    }
    return T(1.0 - t2);
}
 
/*
 * Knab, sampling window from 1983 paper
 */
 
template<typename T>
T _sampling_window(T t, T halfwidth, T bandwidth)
{
    using isce3::except::RuntimeError;
    if (!((0.0 < bandwidth) && (bandwidth < 1.0))) {
        throw RuntimeError(ISCE_SRCINFO(), "Require 0 < bandwidth < 1");
    }
    const T c = M_PI * halfwidth * (1.0 - bandwidth);
    const T tf = t / halfwidth;
    std::complex<T> y = sqrt((std::complex<T>) (1.0 - tf * tf));
    T window = std::real<T>(cosh(c * y) / cosh(c));
    if (!std::isfinite(window)) {
        throw RuntimeError(ISCE_SRCINFO(), "Invalid window parameters.");
    }
    return window;
}
 
// call
template<typename T>
T KnabKernel<T>::operator()(double t) const
{
    auto st = isce3::math::sinc<T>(t);
    return _sampling_window(t, this->_halfwidth, this->_bandwidth) * st;
}
 
/*
 * NFFT
 */
 
// constructor
template<typename T>
NFFTKernel<T>::NFFTKernel(int m, int n, int fft_size)
    : Kernel<T>(2 * m + 1), _m(m), _n(n), _fft_size(fft_size)
{
    if ((m < 1) || (n < 1) || (fft_size < 1)) {
        throw isce3::except::LengthError(ISCE_SRCINFO(),
                                        "NFFT parameters must be positive.");
    }
    _b = M_PI * (2.0 - 1.0 * n / fft_size);
    _scale = 1.0 / (M_PI * isce3::math::bessel_i0(_m * _b));
}
 
// call
template<typename T>
T NFFTKernel<T>::operator()(double t) const
{
    T x2 = t * t - _m * _m;
    // x=0
    if (std::abs(x2) < std::numeric_limits<double>::epsilon()) {
        return _b * _scale;
    }
    T out = 1.0;
    if (x2 < 0.0) {
        T x = std::sqrt(std::abs(x2));
        out = std::sinh(_b * x) / x;
    } else {
        T x = std::sqrt(x2);
        out = std::sin(_b * x) / x;
    }
    return _scale * out;
}
 
/*
 * Azimuth autocorrelation
 */
template<typename T>
T AzimuthKernel<T>::operator()(double t) const
{
    const double x = std::abs(t * 2 / this->_halfwidth);
    if (x > 2.0) {
        return T(0.0);
    }
    // Use Horner's method to save a few multiplies.
    if (x < 1.0) {
        return x * x * (0.75 * x - 1.5) + 1.0;
    }
    // else 1 <= abs(x) <= 2
    return x * (x * (-0.25 * x + 1.5) - 3.0) + 2.0;
}
 
/*
 * Tabulated kernel.
 */
 
// Constructor
template<typename T>
template<typename TI>
TabulatedKernel<T>::TabulatedKernel(const Kernel<TI>& kernel, int n)
    : Kernel<T>(kernel.width())
{
    // Need at least two points for linear interpolation.
    if (n < 2) {
        throw isce3::except::LengthError(ISCE_SRCINFO(),
                                        "Require table size >= 2.");
    }
    // Need i+1 < n so linear interp doesn't run off of table.
    _imax = n - 2;
    // Allocate table.
    _table.resize(n);
    // Assume Kernel is even and fill table with f(x) for 0 <= x <= halfwidth.
    const double dx = this->_halfwidth / (n - 1.0);
    _1_dx = 1.0 / dx;
    for (int i = 0; i < n; ++i) {
        double x = i * dx;
        _table[i] = static_cast<T>(kernel(x));
    }
}
 
// call
template<typename T>
T TabulatedKernel<T>::operator()(double x) const
{
    // Return zero outside table.
    auto ax = std::abs(x);
    if (ax > this->_halfwidth) {
        return T(0);
    }
    // Normalize to table sample index.
    auto axn = ax * _1_dx;
    // Determine left side of interval.
    int i = std::floor(axn);
    // Make sure floating point multiply doesn't cause off-by-one.
    i = std::min(i, _imax);
    // Linear interpolation.
    return _table[i] + (axn - i) * (_table[i + 1] - _table[i]);
}
 
template<typename T>
template<typename Tin>
ChebyKernel<T>::ChebyKernel(const Kernel<Tin>& kernel, int n)
    : Kernel<T>(kernel.width())
{
    if (n < 1) {
        throw isce3::except::LengthError(ISCE_SRCINFO(),
                                        "Need at least one coefficient.");
    }
    // Fit a kernel with DCT of fn at Chebyshev zeros.
    // Assume even function and fit on interval [0,width/2] to avoid a bunch
    // of zero coefficients.
    std::vector<T> q(n), fx(n);
    _scale = 4.0 / kernel.width();
    for (int i = 0; i < n; ++i) {
        q[i] = M_PI * (2.0 * i + 1.0) / (2.0 * n);
        // shift & scale [-1,1] to [0,width/2].
        T x = (std::cos(q[i]) + 1.0) / _scale;
        fx[i] = static_cast<T>(kernel(x));
    }
    // FFTW provides DCT with plan_r2r(REDFT10) but this isn't exposed in
    // isce3::core::signal.  Typically we're only fitting a few coefficients
    // anyway so just implement O(n^2) algorithm.
    _coeffs.resize(n);
    for (int i = 0; i < n; ++i) {
        _coeffs[i] = 0.0;
        for (int j = 0; j < n; ++j) {
            T w = std::cos(i * q[j]);
            _coeffs[i] += w * fx[j];
        }
        _coeffs[i] *= 2.0 / n;
    }
    _coeffs[0] *= 0.5;
}
 
template<typename T>
T ChebyKernel<T>::operator()(double x) const
{
    // Careful to avoid weird stuff outside [-1,1] definition.
    const auto ax = std::abs(x);
    if (ax > this->_halfwidth) {
        return T(0);
    }
    // Map [0,L/2] to [-1,1]
    const T q = (ax * _scale) - T(1);
    const T twoq = T(2) * q;
    const int n = _coeffs.size();
    // Clenshaw algorithm for two term recurrence.
    T bk = 0, bk1 = 0, bk2 = 0;
    for (int i = n - 1; i > 0; --i) {
        bk = _coeffs[i] + twoq * bk1 - bk2;
        bk2 = bk1;
        bk1 = bk;
    }
    return _coeffs[0] + q * bk1 - bk2;
}
 
}} // namespace isce3::core