Skip to content

Commit

Permalink
Added fftfreq() generator
Browse files Browse the repository at this point in the history
  • Loading branch information
cliffburdick committed Jun 6, 2023
1 parent e016476 commit 31298b6
Show file tree
Hide file tree
Showing 7 changed files with 200 additions and 71 deletions.
1 change: 1 addition & 0 deletions docs_input/api/tensorgenerators.rst
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,7 @@ since their only purpose is to compute a single value at a particular location b
.. doxygenfunction:: matx::linspace(const index_t (&s)[RANK], T first, T last)
.. doxygenfunction:: matx::logspace(ShapeType &&s, T first, T last)
.. doxygenfunction:: matx::logspace(const index_t (&s)[RANK], T first, T last)
.. doxygenfunction:: matx::fftfreq(index_t n, float d = 1.0)
.. doxygenfunction:: matx::meshgrid(Ts&&... ts)
.. doxygenfunction:: matx::chirp(SpaceOp t, FreqType f0, typename SpaceOp::scalar_type t1, FreqType f1, ChirpMethod method )
.. doxygenfunction:: matx::chirp(index_t num, TimeType last, FreqType f0, TimeType t1, FreqType f1, ChirpMethod method)
Expand Down
136 changes: 69 additions & 67 deletions include/matx/core/tensor_utils.h
Original file line number Diff line number Diff line change
Expand Up @@ -671,6 +671,75 @@ namespace detail {

static constexpr bool PRINT_ON_DEVICE = false; ///< print() uses printf on device


/**
* @brief Print a tensor's values to stdout
*
* This is the interal `Print()` takes integral values for each index, and prints that as many values
* in each dimension as the arguments specify. For example:
*
* `a.Print(2, 3, 2);`
*
* Will print 2 values of the first, 3 values of the second, and 2 values of the third dimension
* of a 3D tensor. The number of parameters must match the rank of the tensor. A special value of
* 0 can be used if the entire tensor should be printed:
*
* `a.Print(0, 0, 0);` // Prints the whole tensor
*
* For more fine-grained printing, see the over `Print()` overloads.
*
* @tparam Args Integral argument types
* @param op input Operator
* @param dims Number of values to print for each dimension
*/
template <typename Op, typename... Args,
std::enable_if_t<((std::is_integral_v<Args>)&&...) &&
(Op::Rank() == 0 || sizeof...(Args) > 0),
bool> = true>
void PrintData(const Op &op, Args... dims) {
MATX_NVTX_START("", matx::MATX_NVTX_LOG_API)

#ifdef __CUDACC__
cudaDeviceSynchronize();
if constexpr (is_tensor_view_v<Op>) {
auto kind = GetPointerKind(op.Data());

// Try to get pointer from cuda
if (kind == MATX_INVALID_MEMORY) {
CUmemorytype mtype;
void *data[] = {&mtype};
CUpointer_attribute attrs[] = {CU_POINTER_ATTRIBUTE_MEMORY_TYPE};
auto ret = cuPointerGetAttributes(1,
&attrs[0],
data,
reinterpret_cast<CUdeviceptr>(op.Data()));
MATX_ASSERT_STR_EXP(ret, CUDA_SUCCESS, matxCudaError, "Failed to get memory type");
MATX_ASSERT_STR(mtype == CU_MEMORYTYPE_HOST || mtype == 0, matxNotSupported, "Invalid memory type for printing");

detail::InternalPrint(op, dims...);
}
else if (kind == MATX_INVALID_MEMORY || HostPrintable(kind)) {
detail::InternalPrint(op, dims...);
}
else if (DevicePrintable(kind)) {
if constexpr (PRINT_ON_DEVICE) {
PrintKernel<<<1, 1>>>(op, dims...);
}
else {
auto tmpv = make_tensor<typename Op::scalar_type>(op.Shape());
(tmpv = op).run();
PrintData(tmpv, dims...);
}
}
}
else {
InternalPrint(op, dims...);
}
#else
InternalPrint(op, dims...);
#endif
}

/**
* @brief Print a tensor's values to stdout
*
Expand Down Expand Up @@ -777,73 +846,6 @@ void print(const Op &op, Args... dims)

}

/**
* @brief Print a tensor's values to stdout
*
* This is the interal `Print()` takes integral values for each index, and prints that as many values
* in each dimension as the arguments specify. For example:
*
* `a.Print(2, 3, 2);`
*
* Will print 2 values of the first, 3 values of the second, and 2 values of the third dimension
* of a 3D tensor. The number of parameters must match the rank of the tensor. A special value of
* 0 can be used if the entire tensor should be printed:
*
* `a.Print(0, 0, 0);` // Prints the whole tensor
*
* For more fine-grained printing, see the over `Print()` overloads.
*
* @tparam Args Integral argument types
* @param op input Operator
* @param dims Number of values to print for each dimension
*/
template <typename Op, typename... Args,
std::enable_if_t<((std::is_integral_v<Args>)&&...) &&
(Op::Rank() == 0 || sizeof...(Args) > 0),
bool> = true>
void PrintData(const Op &op, Args... dims) {
MATX_NVTX_START("", matx::MATX_NVTX_LOG_API)

#ifdef __CUDACC__
cudaDeviceSynchronize();
if constexpr (is_tensor_view_v<Op>) {
auto kind = GetPointerKind(op.Data());

// Try to get pointer from cuda
if (kind == MATX_INVALID_MEMORY) {
CUmemorytype mtype;
void *data[] = {&mtype};
CUpointer_attribute attrs[] = {CU_POINTER_ATTRIBUTE_MEMORY_TYPE};
auto ret = cuPointerGetAttributes(1,
&attrs[0],
data,
reinterpret_cast<CUdeviceptr>(op.Data()));
MATX_ASSERT_STR_EXP(ret, CUDA_SUCCESS, matxCudaError, "Failed to get memory type");
MATX_ASSERT_STR(mtype == CU_MEMORYTYPE_HOST || mtype == 0, matxNotSupported, "Invalid memory type for printing");

detail::InternalPrint(op, dims...);
}
else if (kind == MATX_INVALID_MEMORY || HostPrintable(kind)) {
detail::InternalPrint(op, dims...);
}
else if (DevicePrintable(kind)) {
if constexpr (PRINT_ON_DEVICE) {
PrintKernel<<<1, 1>>>(op, dims...);
}
else {
auto tmpv = make_tensor<typename Op::scalar_type>(op.Shape());
(tmpv = op).run();
PrintData(tmpv, dims...);
}
}
}
else {
InternalPrint(op, dims...);
}
#else
InternalPrint(op, dims...);
#endif
}

/**
* @brief Print a tensor's all values to stdout
Expand Down
82 changes: 82 additions & 0 deletions include/matx/generators/fftfreq.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,82 @@
////////////////////////////////////////////////////////////////////////////////
// BSD 3-Clause License
//
// Copyright (c) 2021, NVIDIA Corporation
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// 1. Redistributions of source code must retain the above copyright notice, this
// list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
//
// 3. Neither the name of the copyright holder nor the names of its
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
// DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
// FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
// DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
// CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
// OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
/////////////////////////////////////////////////////////////////////////////////

#pragma once

#include "matx/generators/generator1d.h"

namespace matx
{
namespace detail {
template <class T> class FFTFreqOp {
private:
index_t n_;
float d_;

public:
using scalar_type = T;
using matxop = bool;

__MATX_INLINE__ std::string str() const { return "fftfreq"; }

inline FFTFreqOp(index_t n, float d = 1.0)
{
n_ = n;
d_ = d;
}

__MATX_DEVICE__ __MATX_HOST__ __MATX_INLINE__ T operator()(index_t idx) const {
index_t offset = idx >= (n_+1)/2 ? -n_ : 0;
return static_cast<T>(idx + offset) / (d_*(T)n_);
}
};
}


/**
* @brief Return FFT sample frequencies
*
* Returns the bin centers in cycles/unit of the sampling frequency known by the user.
*
* @tparam T Type of output
* @param n Number of elements
* @param d Sample spacing (defaults to 1.0)
* @return Operator with sampling frequencies
*/
template <typename T = float>
inline auto fftfreq(index_t n, float d = 1.0)
{
detail::FFTFreqOp<T> l(n, d);
std::array<index_t, 1> s{n};
return detail::matxGenerator1D_t<detail::FFTFreqOp<T>, 0, decltype(s)>(std::move(s), l);
}
} // end namespace matx
1 change: 1 addition & 0 deletions include/matx/generators/generators.h
Original file line number Diff line number Diff line change
Expand Up @@ -47,3 +47,4 @@
#include "matx/generators/random.h"
#include "matx/generators/range.h"
#include "matx/generators/zeros.h"
#include "matx/generators/fftfreq.h"
6 changes: 3 additions & 3 deletions include/matx/generators/linspace.h
Original file line number Diff line number Diff line change
Expand Up @@ -103,10 +103,10 @@ namespace matx
* spaced apart over the set of values. Distance is determined
* by the shape and selected dimension.
*
* @tparam T Operator type
* @tparam Dim Dimension to operate over
* @tparam ShapeType Shape type
* @param s Shape object
* @tparam RANK Rank of shape
* @tparam T Operator type
* @param s Array of sizes
* @param first First value
* @param last Last value
* @return Operator with linearly-spaced values
Expand Down
25 changes: 25 additions & 0 deletions test/00_operators/GeneratorTests.cu
Original file line number Diff line number Diff line change
Expand Up @@ -218,6 +218,31 @@ TEST(OperatorTests, MeshGrid)
MATX_EXIT_HANDLER();
}

TYPED_TEST(BasicGeneratorTestsFloatNonComplex, FFTFreq)
{
MATX_ENTER_HANDLER();
auto pb = std::make_unique<detail::MatXPybind>();
pb->template InitAndRunTVGenerator<TypeParam>(
"01_signal", "fftfreq", "run", {100});


auto t1 = make_tensor<TypeParam>({100});
auto t2 = make_tensor<TypeParam>({101});

(t1 = fftfreq(t1.Size(0))).run();
cudaStreamSynchronize(0);
MATX_TEST_ASSERT_COMPARE(pb, t1, "F1", 0.1);

(t2 = fftfreq(t2.Size(0))).run();
cudaStreamSynchronize(0);
MATX_TEST_ASSERT_COMPARE(pb, t2, "F2", 0.1);

(t1 = fftfreq(t1.Size(0), 0.5)).run();
cudaStreamSynchronize(0);
MATX_TEST_ASSERT_COMPARE(pb, t1, "F3", 0.1);

MATX_EXIT_HANDLER();
}


TYPED_TEST(BasicGeneratorTestsAll, Zeros)
Expand Down
20 changes: 19 additions & 1 deletion test/test_vectors/generators/01_signal.py
Original file line number Diff line number Diff line change
Expand Up @@ -40,4 +40,22 @@ def run(self):

return {
'Y': Y
}
}

class fftfreq:
def __init__(self, dtype: str, size: List[int]):
self.size = size
self.dtype = dtype

def run(self):
N = self.size[0]

F1 = sf.fftfreq(N)
F2 = sf.fftfreq(N+1)
F3 = sf.fftfreq(N, 0.5)

return {
'F1':F1,
'F2':F2,
'F3':F3,
}

0 comments on commit 31298b6

Please sign in to comment.