diff --git a/CMakeLists.txt b/CMakeLists.txt
index 9d9ed8d33dbbee15ea00770ae4d71248528cecbd..1cabb5c6f0908840c152c7bc70e01f3d93033ae3 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -4,7 +4,7 @@
 # Top level CMakeLists.txt file for SchaapCommon
 cmake_minimum_required(VERSION 3.13) # CMP0079 needs CMake 3.13
 
-set(SCHAAPCOMMON_ALL_MODULES facets fft fitters h5parm ducc0)
+set(SCHAAPCOMMON_ALL_MODULES facets fitters h5parm math ducc0)
 
 # Define a couple of dedicated variables, to prevent that standard CMAKE
 # variables get screwed up when schaapcommon included in another package
diff --git a/include/schaapcommon/fitters/gaussianfitter.h b/include/schaapcommon/fitters/gaussianfitter.h
index 160f2511d5c87cddb6e0ba88b724daf77d4bf511..aa744c11af9d3a65219a0c5bcfaa41f2524244b0 100644
--- a/include/schaapcommon/fitters/gaussianfitter.h
+++ b/include/schaapcommon/fitters/gaussianfitter.h
@@ -6,15 +6,11 @@
 
 #include <cstring>
 
+#include "../math/ellipse.h"
+
 namespace schaapcommon {
 namespace fitters {
 
-struct Ellipse {
-  double major;
-  double minor;
-  double position_angle;
-};
-
 /**
  * Deconvolve two Gaussian kernels. The result convolved with @p smallest
  * results again in @p largest. If @p smallest is larger than @p largest,
@@ -23,16 +19,18 @@ struct Ellipse {
  * As the output axes are returned with the same units as the input axes,
  * the input may be both in FWHM or both in sigma.
  */
-Ellipse DeconvolveGaussian(const Ellipse& largest, const Ellipse& smallest);
+math::Ellipse DeconvolveGaussian(const math::Ellipse& largest,
+                                 const math::Ellipse& smallest);
 
 /**
  * Fit Gaussian parameters to the shape in the center of the image.
  * @param beam_est provides the initial value for the beam's size. It must be
  * somewhat close and must be non-zero for the algorithm to converge (quickly).
  */
-Ellipse Fit2DGaussianCentred(const float* image, size_t width, size_t height,
-                             double beam_est, double box_scale_factor = 10.0,
-                             bool verbose = false);
+math::Ellipse Fit2DGaussianCentred(const float* image, size_t width,
+                                   size_t height, double beam_est,
+                                   double box_scale_factor = 10.0,
+                                   bool verbose = false);
 
 /**
  * Fits a circular (one parameter) Gaussian to the shape in the centre of the
diff --git a/include/schaapcommon/fft/compositefft.h b/include/schaapcommon/math/compositefft.h
similarity index 88%
rename from include/schaapcommon/fft/compositefft.h
rename to include/schaapcommon/math/compositefft.h
index 4d9a361f740bc192a14e00a64de3d1cd79493df2..d45708271097cebed7e6fa066cea8910deb67120 100644
--- a/include/schaapcommon/fft/compositefft.h
+++ b/include/schaapcommon/math/compositefft.h
@@ -7,8 +7,7 @@
 #include <fftw3.h>
 #include <aocommon/staticfor.h>
 
-namespace schaapcommon {
-namespace fft {
+namespace schaapcommon::math {
 void FftR2CComposite(fftwf_plan plan_r2c, fftwf_plan plan_c2c,
                      size_t image_height, size_t image_width, const float* in,
                      fftwf_complex* out, aocommon::StaticFor<size_t>& loop);
@@ -17,7 +16,6 @@ void FftC2RComposite(fftwf_plan plan_c2c, fftwf_plan plan_c2r,
                      size_t image_height, size_t image_width,
                      const fftwf_complex* in, float* out,
                      aocommon::StaticFor<size_t>& loop);
-}  // namespace fft
-}  // namespace schaapcommon
+}  // namespace schaapcommon::math
 
-#endif
\ No newline at end of file
+#endif
diff --git a/include/schaapcommon/fft/convolution.h b/include/schaapcommon/math/convolution.h
similarity index 96%
rename from include/schaapcommon/fft/convolution.h
rename to include/schaapcommon/math/convolution.h
index 00cca6d20839743dc3f4246442b4d32165e72818..d88add1e8afb727defb537a1e9b05331ef96f525 100644
--- a/include/schaapcommon/fft/convolution.h
+++ b/include/schaapcommon/math/convolution.h
@@ -6,8 +6,7 @@
 
 #include <cstring>
 
-namespace schaapcommon {
-namespace fft {
+namespace schaapcommon::math {
 
 /**
  * @brief Make the FFTW float planner thread safe.
@@ -56,7 +55,6 @@ void PrepareConvolutionKernel(float* dest, const float* source,
  */
 void Convolve(float* image, const float* kernel, size_t image_width,
               size_t image_height);
-}  // namespace fft
-}  // namespace schaapcommon
+}  // namespace schaapcommon::math
 
 #endif
diff --git a/include/schaapcommon/math/drawgaussian.h b/include/schaapcommon/math/drawgaussian.h
new file mode 100644
index 0000000000000000000000000000000000000000..8ce6d45c2dcc27b2ca9988e7469f50b2b3153313
--- /dev/null
+++ b/include/schaapcommon/math/drawgaussian.h
@@ -0,0 +1,22 @@
+#ifndef SCHAAPCOMMON_MATH_DRAW_GAUSSIAN_H_
+#define SCHAAPCOMMON_MATH_DRAW_GAUSSIAN_H_
+
+#include <cstring>
+
+#include "ellipse.h"
+
+namespace schaapcommon::math {
+
+/**
+ * Draw a two-dimensional Gaussian onto an lm-space image.
+ */
+void DrawGaussianToLm(float* image_data, size_t image_width,
+                      size_t image_height, long double ra, long double dec,
+                      long double pixel_scale_l, long double pixel_scale_m,
+                      long double l_shift, long double m_shift,
+                      long double position_ra, long double position_dec,
+                      const Ellipse& ellipse, long double integrated_flux);
+
+}  // namespace schaapcommon::math
+
+#endif
diff --git a/include/schaapcommon/math/ellipse.h b/include/schaapcommon/math/ellipse.h
new file mode 100644
index 0000000000000000000000000000000000000000..620cc3690d8f92fcf4e54fcf7f7a42b1f8a61c05
--- /dev/null
+++ b/include/schaapcommon/math/ellipse.h
@@ -0,0 +1,14 @@
+#ifndef SCHAAPCOMMON_MATH_ELLIPSE_H_
+#define SCHAAPCOMMON_MATH_ELLIPSE_H_
+
+namespace schaapcommon::math {
+
+struct Ellipse {
+  double major;
+  double minor;
+  double position_angle;
+};
+
+}  // namespace schaapcommon::math
+
+#endif
diff --git a/include/schaapcommon/fft/resampler.h b/include/schaapcommon/math/resampler.h
similarity index 96%
rename from include/schaapcommon/fft/resampler.h
rename to include/schaapcommon/math/resampler.h
index 3f8684a31c904154f04b33eccdff9a5a19eca98b..03cf7fce8345f3b5d084f56ff4682f9d773d13d8 100644
--- a/include/schaapcommon/fft/resampler.h
+++ b/include/schaapcommon/math/resampler.h
@@ -14,8 +14,7 @@
 
 #include <fftw3.h>
 
-namespace schaapcommon {
-namespace fft {
+namespace schaapcommon::math {
 
 class Resampler {
  private:
@@ -107,7 +106,7 @@ class Resampler {
   aocommon::Lane<Task> tasks_;
   std::vector<std::thread> threads_;
 };
-}  // namespace fft
-}  // namespace schaapcommon
+
+}  // namespace schaapcommon::math
 
 #endif
diff --git a/include/schaapcommon/fft/restoreimage.h b/include/schaapcommon/math/restoreimage.h
similarity index 94%
rename from include/schaapcommon/fft/restoreimage.h
rename to include/schaapcommon/math/restoreimage.h
index d896bb884504bf5dd6dd7018228236488ae476b6..434fa029ff8265b6df6c25d573495928cce07a09 100644
--- a/include/schaapcommon/fft/restoreimage.h
+++ b/include/schaapcommon/math/restoreimage.h
@@ -6,8 +6,7 @@
 
 #include <cstddef>
 
-namespace schaapcommon {
-namespace fft {
+namespace schaapcommon::math {
 
 /**
  * @brief Restore a diffuse image (e.g. produced with a multi-scale clean) using
@@ -34,6 +33,5 @@ void RestoreImage(float* image_data, const float* model_data,
                   long double beam_position_angle, long double pixel_scale_l,
                   long double pixel_scale_m);
 
-}  // namespace fft
-}  // namespace schaapcommon
+}  // namespace schaapcommon::math
 #endif
diff --git a/src/fitters/gaussianfitter.cc b/src/fitters/gaussianfitter.cc
index 1a6750c6c9275f34e3d668ac0bd88ac31a0e8506..5720bca32a082b8057774e4ee4a8742eefe49a39 100644
--- a/src/fitters/gaussianfitter.cc
+++ b/src/fitters/gaussianfitter.cc
@@ -16,6 +16,8 @@
 
 using aocommon::Matrix2x2;
 
+using schaapcommon::math::Ellipse;
+
 namespace schaapcommon::fitters {
 namespace {
 
diff --git a/src/fitters/test/tgaussianfitter.cc b/src/fitters/test/tgaussianfitter.cc
index 52ebb200d3cd50da866f7dcfe4b9c0be9ca165c1..8dc7c73fda686258e33a9b280adaaea5699276fc 100644
--- a/src/fitters/test/tgaussianfitter.cc
+++ b/src/fitters/test/tgaussianfitter.cc
@@ -8,8 +8,8 @@
 #include <aocommon/image.h>
 #include <aocommon/threadpool.h>
 
-#include "../../../include/schaapcommon/fft/convolution.h"
-#include "../../../include/schaapcommon/fft/restoreimage.h"
+#include "../../../include/schaapcommon/math/convolution.h"
+#include "../../../include/schaapcommon/math/restoreimage.h"
 
 namespace {
 constexpr size_t kThreadCount = 2;
@@ -18,7 +18,7 @@ constexpr size_t kHeight = 64;
 constexpr long double kPixelSize = 1 /*amin*/ * (M_PI / 180.0 / 60.0);
 }  // namespace
 
-using schaapcommon::fitters::Ellipse;
+using schaapcommon::math::Ellipse;
 
 BOOST_AUTO_TEST_SUITE(gaussian_fitter)
 
@@ -36,8 +36,8 @@ BOOST_AUTO_TEST_CASE(fit) {
     const long double beam_position_angle =
         beam_position_angle_index * M_PI / 10.0;
 
-    schaapcommon::fft::MakeFftwfPlannerThreadSafe();
-    schaapcommon::fft::RestoreImage(
+    schaapcommon::math::MakeFftwfPlannerThreadSafe();
+    schaapcommon::math::RestoreImage(
         restored.Data(), model.Data(), width, height, beam_major, beam_minor,
         beam_position_angle, kPixelSize, kPixelSize);
 
@@ -68,8 +68,8 @@ BOOST_AUTO_TEST_CASE(fit_full) {
     long double beam_minor = 5.0L * kPixelSize;
     long double beam_position_angle = beam_position_angle_index * M_PI / 10.0;
 
-    schaapcommon::fft::MakeFftwfPlannerThreadSafe();
-    schaapcommon::fft::RestoreImage(
+    schaapcommon::math::MakeFftwfPlannerThreadSafe();
+    schaapcommon::math::RestoreImage(
         restored.Data(), model.Data(), width, height, beam_major, beam_minor,
         beam_position_angle, kPixelSize, kPixelSize);
 
@@ -114,10 +114,10 @@ BOOST_AUTO_TEST_CASE(fit_circular) {
   const long double beam_minor = 4.0L * kPixelSize;
   const long double beam_position_angle = 0.0;
   const long double estimated_beam_pixel = 1.0;  // this is on purpose way off
-  schaapcommon::fft::MakeFftwfPlannerThreadSafe();
-  schaapcommon::fft::RestoreImage(restored.Data(), model.Data(), kWidth,
-                                  kHeight, beam_major, beam_minor,
-                                  beam_position_angle, kPixelSize, kPixelSize);
+  schaapcommon::math::MakeFftwfPlannerThreadSafe();
+  schaapcommon::math::RestoreImage(restored.Data(), model.Data(), kWidth,
+                                   kHeight, beam_major, beam_minor,
+                                   beam_position_angle, kPixelSize, kPixelSize);
 
   const Ellipse result = schaapcommon::fitters::Fit2DGaussianCentred(
       restored.Data(), restored.Width(), restored.Height(),
@@ -152,10 +152,10 @@ BOOST_AUTO_TEST_CASE(fit_small_beam) {
   const long double beam_position_angle = 0.0;
   const long double estimated_beam_pixel = 1.0;  // this is on purpose way off
 
-  schaapcommon::fft::MakeFftwfPlannerThreadSafe();
-  schaapcommon::fft::RestoreImage(restored.Data(), model.Data(), kWidth,
-                                  kHeight, beam_major, beam_minor,
-                                  beam_position_angle, kPixelSize, kPixelSize);
+  schaapcommon::math::MakeFftwfPlannerThreadSafe();
+  schaapcommon::math::RestoreImage(restored.Data(), model.Data(), kWidth,
+                                   kHeight, beam_major, beam_minor,
+                                   beam_position_angle, kPixelSize, kPixelSize);
 
   Ellipse ellipse = schaapcommon::fitters::Fit2DGaussianCentred(
       restored.Data(), restored.Width(), restored.Height(),
diff --git a/src/fft/CMakeLists.txt b/src/math/CMakeLists.txt
similarity index 88%
rename from src/fft/CMakeLists.txt
rename to src/math/CMakeLists.txt
index 8286cca90da0e2675aadf257ac4423c48ba16b73..99a95f1a8b69d876a2e48372cb60af45c11bf30f 100644
--- a/src/fft/CMakeLists.txt
+++ b/src/math/CMakeLists.txt
@@ -19,13 +19,15 @@ get_filename_component(MODULE ${CMAKE_CURRENT_SOURCE_DIR} NAME)
 set(PUBLIC_HEADER_DIR ${SCHAAPCOMMON_SOURCE_DIR}/include/schaapcommon/${MODULE})
 
 set(PUBLIC_HEADERS
-    ${PUBLIC_HEADER_DIR}/convolution.h ${PUBLIC_HEADER_DIR}/resampler.h
+    ${PUBLIC_HEADER_DIR}/convolution.h ${PUBLIC_HEADER_DIR}/drawgaussian.h
+    ${PUBLIC_HEADER_DIR}/ellipse.h ${PUBLIC_HEADER_DIR}/resampler.h
     ${PUBLIC_HEADER_DIR}/restoreimage.h)
 
 target_sources(
   ${SCHAAPCOMMON_PROJECT_NAME}
   PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/convolution.cc
           ${CMAKE_CURRENT_SOURCE_DIR}/compositefft.cc
+          ${CMAKE_CURRENT_SOURCE_DIR}/drawgaussian.cc
           ${CMAKE_CURRENT_SOURCE_DIR}/resampler.cc
           ${CMAKE_CURRENT_SOURCE_DIR}/restoreimage.cc)
 
diff --git a/src/fft/compositefft.cc b/src/math/compositefft.cc
similarity index 98%
rename from src/fft/compositefft.cc
rename to src/math/compositefft.cc
index 50efbd38f6b1057b3a9d46b23226894a21a9b895..e0a27b58219e807e66e5b1adf0cb357623eca19f 100644
--- a/src/fft/compositefft.cc
+++ b/src/math/compositefft.cc
@@ -19,8 +19,7 @@ size_t RoundUp(size_t a, size_t b) { return ((a + b) / b) * b; }
 
 }  // namespace
 
-namespace schaapcommon {
-namespace fft {
+namespace schaapcommon::math {
 void FftR2CComposite(fftwf_plan plan_r2c, fftwf_plan plan_c2c,
                      size_t image_height, size_t image_width, const float* in,
                      fftwf_complex* out, aocommon::StaticFor<size_t>& loop) {
@@ -162,5 +161,4 @@ void FftC2RComposite(fftwf_plan plan_c2c, fftwf_plan plan_c2r,
 
   fftwf_free(temp1);
 }
-}  // namespace fft
-}  // namespace schaapcommon
\ No newline at end of file
+}  // namespace schaapcommon::math
diff --git a/src/fft/convolution.cc b/src/math/convolution.cc
similarity index 98%
rename from src/fft/convolution.cc
rename to src/math/convolution.cc
index cae161e391e8496d76d1e772631952b8a6373c78..2a6e02716e46e3a95d9c41f7c25c1b3e9212cea0 100644
--- a/src/fft/convolution.cc
+++ b/src/math/convolution.cc
@@ -15,8 +15,7 @@
 
 #include <iostream>
 
-namespace schaapcommon {
-namespace fft {
+namespace schaapcommon::math {
 
 void MakeFftwfPlannerThreadSafe() { fftwf_make_planner_thread_safe(); }
 
@@ -167,5 +166,4 @@ void Convolve(float* image, const float* kernel, size_t image_width,
   fftwf_destroy_plan(plan_c2r);
 }
 
-}  // namespace fft
-}  // namespace schaapcommon
+}  // namespace schaapcommon::math
diff --git a/src/math/drawgaussian.cc b/src/math/drawgaussian.cc
new file mode 100644
index 0000000000000000000000000000000000000000..345d6a707d6516d2acf7c0843592cfdd84b6d609
--- /dev/null
+++ b/src/math/drawgaussian.cc
@@ -0,0 +1,104 @@
+#include "drawgaussian.h"
+
+#include <cmath>
+
+#include <aocommon/imagecoordinates.h>
+
+using aocommon::ImageCoordinates;
+
+namespace {
+inline long double Gaussian(long double x, long double sigma) {
+  // Evaluate unnormalized Gaussian (unit valued peak, but not a unit integral.
+  const long double xi = x / sigma;
+  return std::exp(-0.5L * xi * xi);
+}
+}  // namespace
+
+namespace schaapcommon::math {
+
+void DrawGaussianToLm(float* image_data, size_t image_width,
+                      size_t image_height, long double ra, long double dec,
+                      long double pixel_scale_l, long double pixel_scale_m,
+                      long double l_shift, long double m_shift,
+                      long double position_ra, long double position_dec,
+                      const Ellipse& ellipse, long double integrated_flux) {
+  // Using the FWHM formula for a Gaussian:
+  const long double fwhm_constant = 2.0L * std::sqrt(2.0L * M_LN2);
+  const long double sigma_major_axis = ellipse.major / fwhm_constant;
+  const long double sigma_minor_axis = ellipse.minor / fwhm_constant;
+  // TODO this won't work for non-equally spaced dimensions
+  const long double min_pixel_scale = std::min(pixel_scale_l, pixel_scale_m);
+
+  // Position angle is angle from North:
+  const long double angle = ellipse.position_angle + M_PI_2;
+  const long double cos_angle = std::cos(angle);
+  const long double sin_angle = std::sin(angle);
+
+  // Make rotation matrix
+  long double transf[4];
+  transf[0] = cos_angle;
+  transf[1] = -sin_angle;
+  transf[2] = sin_angle;
+  transf[3] = cos_angle;
+
+  const double sigma_max = std::max(std::fabs(sigma_major_axis * transf[0]),
+                                    std::fabs(sigma_major_axis * transf[1]));
+  // Multiply with scaling matrix to make variance 1.
+  transf[0] /= sigma_major_axis;
+  transf[1] /= sigma_major_axis;
+  transf[2] /= sigma_minor_axis;
+  transf[3] /= sigma_minor_axis;
+  const int bounding_box_size = std::ceil(sigma_max * 20.0 / min_pixel_scale);
+  long double source_l;
+  long double source_m;
+  ImageCoordinates::RaDecToLM(position_ra, position_dec, ra, dec, source_l,
+                              source_m);
+
+  // Calculate the bounding box
+  int source_x;
+  int source_y;
+  ImageCoordinates::LMToXY<long double>(
+      source_l - l_shift, source_m - m_shift, pixel_scale_l, pixel_scale_m,
+      image_width, image_height, source_x, source_y);
+  const int x_left =
+      std::clamp(source_x - bounding_box_size, 0, int(image_width));
+  const int x_right =
+      std::clamp(source_x + bounding_box_size, x_left, int(image_width));
+  const int y_top =
+      std::clamp(source_y - bounding_box_size, 0, int(image_height));
+  const int y_bottom =
+      std::clamp(source_y + bounding_box_size, y_top, int(image_height));
+
+  std::vector<double> values;
+  double flux_sum = 0.0;
+  for (int y = y_top; y != y_bottom; ++y) {
+    for (int x = x_left; x != x_right; ++x) {
+      long double l, m;
+      ImageCoordinates::XYToLM<long double>(x, y, pixel_scale_l, pixel_scale_m,
+                                            image_width, image_height, l, m);
+      l += l_shift;
+      m += m_shift;
+      const long double l_transf =
+          (l - source_l) * transf[0] + (m - source_m) * transf[1];
+      const long double m_transf =
+          (l - source_l) * transf[2] + (m - source_m) * transf[3];
+      const long double dist =
+          std::sqrt(l_transf * l_transf + m_transf * m_transf);
+      long double v = Gaussian(dist, 1.0L);
+      flux_sum += static_cast<double>(v);
+      values.emplace_back(v);
+    }
+  }
+  const double* iter = values.data();
+  const double factor = integrated_flux / flux_sum;
+  for (int y = y_top; y != y_bottom; ++y) {
+    float* image_data_ptr = image_data + y * image_width + x_left;
+    for (int x = x_left; x != x_right; ++x) {
+      (*image_data_ptr) += *iter * factor;
+      ++image_data_ptr;
+      ++iter;
+    }
+  }
+}
+
+}  // namespace schaapcommon::math
diff --git a/src/fft/resampler.cc b/src/math/resampler.cc
similarity index 98%
rename from src/fft/resampler.cc
rename to src/math/resampler.cc
index c5fbf288303e298b8293b965e62032946f4a7c39..f57dc0d3393926ea7a43c8edaa93fbe8ed16601c 100644
--- a/src/fft/resampler.cc
+++ b/src/math/resampler.cc
@@ -3,8 +3,7 @@
 #include <complex>
 
 using aocommon::WindowFunction;
-namespace schaapcommon {
-namespace fft {
+namespace schaapcommon::math {
 
 Resampler::Resampler(size_t input_width, size_t input_height,
                      size_t output_width, size_t output_height,
@@ -248,5 +247,4 @@ void Resampler::UnapplyWindow(float* data) const {
     data[i] /= window_out_[i];
   }
 }
-}  // namespace fft
-}  // namespace schaapcommon
+}  // namespace schaapcommon::math
diff --git a/src/fft/restoreimage.cc b/src/math/restoreimage.cc
similarity index 89%
rename from src/fft/restoreimage.cc
rename to src/math/restoreimage.cc
index ea7312bcb03a064c6a8826731e4c999f0c807be4..4920465ab657d0da0d2683761ab31e838bbf5c5d 100644
--- a/src/fft/restoreimage.cc
+++ b/src/math/restoreimage.cc
@@ -12,14 +12,13 @@
 
 #include "convolution.h"
 
-namespace schaapcommon {
-namespace fft {
+namespace schaapcommon::math {
 void RestoreImage(float* image_data, const float* model_data,
                   size_t image_width, size_t image_height,
                   long double beam_major_axis, long double beam_minor_axis,
                   long double beam_position_angle, long double pixel_scale_l,
                   long double pixel_scale_m) {
-  if (beam_major_axis == 0.0 && beam_minor_axis == 0.0) {
+  if (beam_major_axis == 0.0L && beam_minor_axis == 0.0L) {
     for (size_t j = 0; j != image_width * image_height; ++j) {
       image_data[j] += model_data[j];
     }
@@ -83,13 +82,12 @@ void RestoreImage(float* image_data, const float* model_data,
     aocommon::UVector<float> convolved_model(
         model_data, model_data + image_width * image_height);
 
-    schaapcommon::fft::ResizeAndConvolve(convolved_model.data(), image_width,
-                                         image_height, kernel.data(),
-                                         bounding_box_size);
+    schaapcommon::math::ResizeAndConvolve(convolved_model.data(), image_width,
+                                          image_height, kernel.data(),
+                                          bounding_box_size);
     for (size_t j = 0; j != image_width * image_height; ++j) {
       image_data[j] += convolved_model[j];
     }
   }
 }
-}  // namespace fft
-}  // namespace schaapcommon
+}  // namespace schaapcommon::math
diff --git a/src/fft/test/CMakeLists.txt b/src/math/test/CMakeLists.txt
similarity index 69%
rename from src/fft/test/CMakeLists.txt
rename to src/math/test/CMakeLists.txt
index 9f8ccfe376ee87b75d31208a11d7e7965e4c3de5..1fd519713d5a352f185c4239cdffc47d5dfbee66 100644
--- a/src/fft/test/CMakeLists.txt
+++ b/src/math/test/CMakeLists.txt
@@ -3,4 +3,4 @@
 
 include(${SCHAAPCOMMON_SOURCE_DIR}/cmake/unittest.cmake)
 
-add_unittest(fft runtests.cc tconvolution.cc tresampler.cc trestoreimage.cc)
+add_unittest(math runtests.cc tconvolution.cc tresampler.cc trestoreimage.cc)
diff --git a/src/fft/test/runtests.cc b/src/math/test/runtests.cc
similarity index 100%
rename from src/fft/test/runtests.cc
rename to src/math/test/runtests.cc
diff --git a/src/fft/test/tconvolution.cc b/src/math/test/tconvolution.cc
similarity index 84%
rename from src/fft/test/tconvolution.cc
rename to src/math/test/tconvolution.cc
index 490b962472f87d1a2cd689b48a6e935b556595e4..288d0ec0052ce340b00185703151ea6ae020f466 100644
--- a/src/fft/test/tconvolution.cc
+++ b/src/math/test/tconvolution.cc
@@ -25,7 +25,7 @@ BOOST_AUTO_TEST_CASE(prepare_kernel) {
                             {10.0, 11.0, 8.0, 9.0, 14.0, 15.0, 12.0, 13.0, 2.0,
                              3.0, 0.0, 1.0, 6.0, 7.0, 4.0, 5.0});
   aocommon::Image kernel_out(kWidth, kHeight);
-  schaapcommon::fft::PrepareConvolutionKernel(
+  schaapcommon::math::PrepareConvolutionKernel(
       kernel_out.Data(), kernel_in.Data(), kWidth, kHeight);
   for (size_t i = 0; i != kernel_out.Size(); ++i) {
     BOOST_CHECK_CLOSE(kernel_out[i], static_cast<float>(i), 1e-4);
@@ -38,11 +38,11 @@ BOOST_AUTO_TEST_CASE(prepare_small_kernel) {
   aocommon::Image kernel_in(kWidth / 2, kHeight / 2, {16.0, 13.0, 4.0, 1.0});
   aocommon::Image kernel_out(kWidth, kHeight, 0.0);
   BOOST_CHECK_THROW(
-      schaapcommon::fft::PrepareSmallConvolutionKernel(
+      schaapcommon::math::PrepareSmallConvolutionKernel(
           kernel_out.Data(), kWidth, kHeight, kernel_in.Data(), kWidth * 2),
       std::runtime_error);
 
-  schaapcommon::fft::PrepareSmallConvolutionKernel(
+  schaapcommon::math::PrepareSmallConvolutionKernel(
       kernel_out.Data(), kWidth, kHeight, kernel_in.Data(), kWidth / 2);
 
   // kernel_out should have non-zero corner values only
@@ -68,14 +68,14 @@ BOOST_AUTO_TEST_CASE(convolve) {
 
   kernel[kWidth * kHeight / 2 + kWidth / 2] = dirac_scale * 1.0f;
 
-  schaapcommon::fft::MakeFftwfPlannerThreadSafe();
+  schaapcommon::math::MakeFftwfPlannerThreadSafe();
   BOOST_CHECK_THROW(
-      schaapcommon::fft::ResizeAndConvolve(image.Data(), kWidth, kHeight,
-                                           kernel.Data(), kWidth * 2),
+      schaapcommon::math::ResizeAndConvolve(image.Data(), kWidth, kHeight,
+                                            kernel.Data(), kWidth * 2),
       std::runtime_error);
 
-  schaapcommon::fft::ResizeAndConvolve(image.Data(), kWidth, kHeight,
-                                       kernel.Data(), kWidth);
+  schaapcommon::math::ResizeAndConvolve(image.Data(), kWidth, kHeight,
+                                        kernel.Data(), kWidth);
 
   // Convolution with (scaled) dirac kernel should return same signal,
   // scaled by a constant factor.
diff --git a/src/fft/test/tresampler.cc b/src/math/test/tresampler.cc
similarity index 98%
rename from src/fft/test/tresampler.cc
rename to src/math/test/tresampler.cc
index c553ba6f46e09bfac8edb753d0ed355654ccb577..2680871f114e09d83a480c00fb988de9e76194e0 100644
--- a/src/fft/test/tresampler.cc
+++ b/src/math/test/tresampler.cc
@@ -9,7 +9,7 @@
 
 #include <aocommon/image.h>
 
-using schaapcommon::fft::Resampler;
+using schaapcommon::math::Resampler;
 
 namespace {
 constexpr size_t kInputWidth = 16;
@@ -137,4 +137,4 @@ BOOST_AUTO_TEST_CASE(upsample_regular_window) {
   CheckResampledImage(input_image, output_image, scale_output);
 }
 
-BOOST_AUTO_TEST_SUITE_END()
\ No newline at end of file
+BOOST_AUTO_TEST_SUITE_END()
diff --git a/src/fft/test/trestoreimage.cc b/src/math/test/trestoreimage.cc
similarity index 97%
rename from src/fft/test/trestoreimage.cc
rename to src/math/test/trestoreimage.cc
index e90d5a1054a4cfcd531251e7573bfaa471acc2db..1f3d326a96e64e0db5adbe6eb5c5b95a46fafd81 100644
--- a/src/fft/test/trestoreimage.cc
+++ b/src/math/test/trestoreimage.cc
@@ -15,8 +15,8 @@
 
 #include <iostream>
 
-using schaapcommon::fft::MakeFftwfPlannerThreadSafe;
-using schaapcommon::fft::RestoreImage;
+using schaapcommon::math::MakeFftwfPlannerThreadSafe;
+using schaapcommon::math::RestoreImage;
 
 namespace {
 constexpr size_t kWidth = 16;