diff --git a/doc/docs/Python_User_Interface.md b/doc/docs/Python_User_Interface.md
index 0c16a5849..0b3323806 100644
--- a/doc/docs/Python_User_Interface.md
+++ b/doc/docs/Python_User_Interface.md
@@ -6506,6 +6506,7 @@ def __init__(self,
              start_time=-1e+20,
              end_time=1e+20,
              center_frequency=0,
+             fwidth=0,
              **kwargs):
 ```
 
@@ -6537,6 +6538,11 @@ Construct a `CustomSource`.
 + **`center_frequency` [`number`]** — Optional center frequency so that the
   `CustomSource` can be used within an `EigenModeSource`. Defaults to 0.
 
++ **`fwidth` [`number`]** — Optional bandwidth in frequency units.
+  Default is 0. For bandwidth-limited sources, this parameter is used to
+  automatically determine the decimation factor of the time-series updates
+  of the DFT fields monitors (if any).
+
 </div>
 
 </div>
diff --git a/python/adjoint/filter_source.py b/python/adjoint/filter_source.py
index cce5fc881..6c82cff27 100644
--- a/python/adjoint/filter_source.py
+++ b/python/adjoint/filter_source.py
@@ -26,6 +26,10 @@ def __init__(
         self.t = np.arange(0, dt * (self.N), dt)
         self.n = np.arange(self.N)
         f = self.func()
+
+        # frequency bandwidth of the Nuttall window function
+        fwidth = self.nuttall_bandwidth()
+
         self.bf = [
             lambda t, i=i: 0
             if t > self.T else (self.nuttall(t, self.center_frequencies) /
@@ -36,7 +40,8 @@ def __init__(
             CustomSource(src_func=bfi,
                          center_frequency=center_frequency,
                          is_integrated=False,
-                         end_time=self.T) for bfi in self.bf
+                         end_time=self.T,
+                         fwidth=fwidth) for bfi in self.bf
         ]
 
         if time_src:
@@ -58,7 +63,8 @@ def __init__(
         super(FilteredSource, self).__init__(src_func=f,
                                              center_frequency=center_frequency,
                                              is_integrated=False,
-                                             end_time=self.T)
+                                             end_time=self.T,
+                                             fwidth=fwidth)
 
     def cos_window_td(self, a, t, f0):
         cos_sum = np.sum([(-1)**k * a[k] * np.cos(2 * np.pi * t * k / self.T)
@@ -103,6 +109,20 @@ def nuttall_dtft(self, f, f0):
         a = [0.355768, 0.4873960, 0.144232, 0.012604]
         return self.cos_window_fd(a, f, f0)
 
+    ## compute the bandwidth of the DTFT of the Nuttall window function
+    ## (magnitude) assuming it has decayed from its peak value by some
+    ## tolerance by fitting it to an asymptotic power law of the form
+    ## C / f^3 where C is a constant and f is the frequency
+    def nuttall_bandwidth(self):
+        tol = 1e-7
+        fwidth = 1/(self.N * self.dt)
+        frq_inf = 10000*fwidth
+        na_dtft = self.nuttall_dtft(frq_inf, 0)
+        coeff = frq_inf**3 * np.abs(na_dtft)
+        na_dtft_max = self.nuttall_dtft(0, 0)
+        bw = 2 * np.power(coeff / (tol * na_dtft_max), 1/3)
+        return bw.real
+
     def dtft(self, y, f):
         return np.matmul(
             np.exp(1j * 2 * np.pi * f[:, np.newaxis] * np.arange(y.size) *
diff --git a/python/meep-python.hpp b/python/meep-python.hpp
index e254f8f26..0e6c6e278 100644
--- a/python/meep-python.hpp
+++ b/python/meep-python.hpp
@@ -8,8 +8,8 @@ namespace meep {
 class custom_py_src_time : public src_time {
 public:
   custom_py_src_time(PyObject *fun, double st = -infinity, double et = infinity,
-                     std::complex<double> f = 0)
-      : func(fun), freq(f), start_time(float(st)), end_time(float(et)) {
+                     std::complex<double> f = 0, double fw = 0)
+      : func(fun), freq(f), start_time(float(st)), end_time(float(et)), fwidth(fw) {
     SWIG_PYTHON_THREAD_SCOPED_BLOCK;
     Py_INCREF(func);
   }
@@ -50,17 +50,19 @@ class custom_py_src_time : public src_time {
     const custom_py_src_time *tp = dynamic_cast<const custom_py_src_time *>(&t);
     if (tp)
       return (tp->start_time == start_time && tp->end_time == end_time && tp->func == func &&
-              tp->freq == freq);
+              tp->freq == freq && tp->fwidth == fwidth);
     else
       return 0;
   }
   virtual std::complex<double> frequency() const { return freq; }
   virtual void set_frequency(std::complex<double> f) { freq = f; }
+  virtual double get_fwidth() const { return fwidth; };
+  virtual void set_fwidth(double fw) { fwidth = fw; }
 
 private:
   PyObject *func;
   std::complex<double> freq;
-  double start_time, end_time;
+  double start_time, end_time, fwidth;
 };
 
 } // namespace meep
diff --git a/python/source.py b/python/source.py
index 205eab2d7..52d3e3fc6 100644
--- a/python/source.py
+++ b/python/source.py
@@ -269,7 +269,7 @@ class CustomSource(SourceTime):
     [`examples/chirped_pulse.py`](https://github.com/NanoComp/meep/blob/master/python/examples/chirped_pulse.py).
     """
 
-    def __init__(self, src_func, start_time=-1.0e20, end_time=1.0e20, center_frequency=0, **kwargs):
+    def __init__(self, src_func, start_time=-1.0e20, end_time=1.0e20, center_frequency=0, fwidth=0, **kwargs):
         """
         Construct a `CustomSource`.
 
@@ -296,13 +296,20 @@ def __init__(self, src_func, start_time=-1.0e20, end_time=1.0e20, center_frequen
 
         + **`center_frequency` [`number`]** — Optional center frequency so that the
           `CustomSource` can be used within an `EigenModeSource`. Defaults to 0.
+
+        + **`fwidth` [`number`]** — Optional bandwidth in frequency units.
+          Default is 0. For bandwidth-limited sources, this parameter is used to
+          automatically determine the decimation factor of the time-series updates
+          of the DFT fields monitors (if any).
         """
         super(CustomSource, self).__init__(**kwargs)
         self.src_func = src_func
         self.start_time = start_time
         self.end_time = end_time
+        self.fwidth = fwidth
         self.center_frequency = center_frequency
-        self.swigobj = mp.custom_py_src_time(src_func, start_time, end_time, center_frequency)
+        self.swigobj = mp.custom_py_src_time(src_func, start_time, end_time,
+                                             center_frequency, fwidth)
         self.swigobj.is_integrated = self.is_integrated
 
 
diff --git a/python/tests/test_adjoint_cyl.py b/python/tests/test_adjoint_cyl.py
index 6f16b0db1..f02727f0c 100644
--- a/python/tests/test_adjoint_cyl.py
+++ b/python/tests/test_adjoint_cyl.py
@@ -26,8 +26,8 @@
 design_region_resolution = int(2*resolution)
 design_r = 4.8
 design_z = 2
-Nx = int(design_region_resolution*design_r)
-Nz = int(design_region_resolution*design_z)
+Nr = int(design_region_resolution*design_r) + 1
+Nz = int(design_region_resolution*design_z) + 1
 
 fcen = 1/1.55
 width = 0.2
@@ -37,20 +37,20 @@
 src = mp.GaussianSource(frequency=fcen,fwidth=fwidth)
 source = [mp.Source(src,component=mp.Er,
                     center=source_center,
-                   size=source_size)]
+                    size=source_size)]
 
 ## random design region
-p = np.random.rand(Nx*Nz)
+p = np.random.rand(Nr*Nz)
 ## random epsilon perturbation for design region
 deps = 1e-5
-dp = deps*np.random.rand(Nx*Nz)
+dp = deps*np.random.rand(Nr*Nz)
 
 
 def forward_simulation(design_params):
-    matgrid = mp.MaterialGrid(mp.Vector3(Nx,0,Nz),
+    matgrid = mp.MaterialGrid(mp.Vector3(Nr,0,Nz),
                               SiO2,
                               Si,
-                              weights=design_params.reshape(Nx,1,Nz))
+                              weights=design_params.reshape(Nr,1,Nz))
 
     geometry = [mp.Block(center=mp.Vector3(0.1+design_r/2,0,0),
                                  size=mp.Vector3(design_r,0,design_z),
@@ -68,9 +68,8 @@ def forward_simulation(design_params):
     far_x = [mp.Vector3(5,0,20)]
     mode = sim.add_near2far(
         frequencies,
-        mp.Near2FarRegion(center=mp.Vector3(0.1+design_r/2,0 ,(sz/2-dpml+design_z/2)/2),size=mp.Vector3(design_r,0,0), weight=+1),
-        decimation_factor=10
-    )
+        mp.Near2FarRegion(center=mp.Vector3(0.1+design_r/2,0,(sz/2-dpml+design_z/2)/2),
+                          size=mp.Vector3(design_r,0,0),weight=+1))
 
     sim.run(until_after_sources=1200)
     Er = sim.get_farfield(mode, far_x[0])
@@ -81,9 +80,13 @@ def forward_simulation(design_params):
 
 def adjoint_solver(design_params):
 
-    design_variables = mp.MaterialGrid(mp.Vector3(Nx,0,Nz),SiO2,Si)
-    design_region = mpa.DesignRegion(design_variables,volume=mp.Volume(center=mp.Vector3(0.1+design_r/2,0,0), size=mp.Vector3(design_r, 0,design_z)))
-    geometry = [mp.Block(center=design_region.center, size=design_region.size, material=design_variables)]
+    design_variables = mp.MaterialGrid(mp.Vector3(Nr,0,Nz),SiO2,Si)
+    design_region = mpa.DesignRegion(design_variables,
+                                     volume=mp.Volume(center=mp.Vector3(0.1+design_r/2,0,0),
+                                                      size=mp.Vector3(design_r,0,design_z)))
+    geometry = [mp.Block(center=design_region.center,
+                         size=design_region.size,
+                         material=design_variables)]
 
     sim = mp.Simulation(cell_size=cell_size,
         boundary_layers=boundary_layers,
@@ -94,8 +97,10 @@ def adjoint_solver(design_params):
         m=m)
 
     far_x = [mp.Vector3(5,0,20)]
-    NearRegions = [mp.Near2FarRegion(center=mp.Vector3(0.1+design_r/2,0 ,(sz/2-dpml+design_z/2)/2),size=mp.Vector3(design_r,0,0), weight=+1)]
-    FarFields = mpa.Near2FarFields(sim, NearRegions ,far_x, decimation_factor=5)
+    NearRegions = [mp.Near2FarRegion(center=mp.Vector3(0.1+design_r/2,0,(sz/2-dpml+design_z/2)/2),
+                                     size=mp.Vector3(design_r,0,0),
+                                     weight=+1)]
+    FarFields = mpa.Near2FarFields(sim, NearRegions ,far_x)
     ob_list = [FarFields]
 
     def J(alpha):
diff --git a/python/tests/test_adjoint_jax.py b/python/tests/test_adjoint_jax.py
index e3bcad587..768d7ac82 100644
--- a/python/tests/test_adjoint_jax.py
+++ b/python/tests/test_adjoint_jax.py
@@ -127,8 +127,7 @@ def build_straight_wg_simulation(
         mpa.EigenmodeCoefficient(simulation,
                                  mp.Volume(center=center, size=monitor_size),
                                  mode=1,
-                                 forward=forward,
-                                 decimation_factor=5)
+                                 forward=forward)
         for center in monitor_centers for forward in [True, False]
     ]
     return simulation, sources, monitors, design_regions, frequencies
diff --git a/python/tests/test_adjoint_solver.py b/python/tests/test_adjoint_solver.py
index 481366cc3..a918372e7 100644
--- a/python/tests/test_adjoint_solver.py
+++ b/python/tests/test_adjoint_solver.py
@@ -81,7 +81,6 @@ def forward_simulation(design_params,mon_type, frequencies=None, use_complex=Fal
                                     mp.ModeRegion(center=mp.Vector3(0.5*sxy-dpml-0.1),
                                                   size=mp.Vector3(0,sxy-2*dpml,0)),
                                     yee_grid=True,
-                                    decimation_factor=10,
                                     eig_parity=eig_parity)
 
     elif mon_type.name == 'DFT':
@@ -89,8 +88,7 @@ def forward_simulation(design_params,mon_type, frequencies=None, use_complex=Fal
                                   frequencies,
                                   center=mp.Vector3(1.25),
                                   size=mp.Vector3(0.25,1,0),
-                                  yee_grid=False,
-                                  decimation_factor=10)
+                                  yee_grid=False)
 
     sim.run(until_after_sources=mp.stop_when_dft_decayed())
 
@@ -145,7 +143,6 @@ def adjoint_solver(design_params, mon_type, frequencies=None, use_complex=False,
                                              mp.Volume(center=mp.Vector3(0.5*sxy-dpml-0.1),
                                                        size=mp.Vector3(0,sxy-2*dpml,0)),
                                              1,
-                                             decimation_factor=5,
                                              eig_parity=eig_parity)]
 
         def J(mode_mon):
@@ -155,8 +152,7 @@ def J(mode_mon):
         obj_list = [mpa.FourierFields(sim,
                                       mp.Volume(center=mp.Vector3(1.25),
                                                 size=mp.Vector3(0.25,1,0)),
-                                      mp.Ez,
-                                      decimation_factor=5)]
+                                      mp.Ez)]
 
         def J(mode_mon):
             return npa.power(npa.abs(mode_mon[:,4,10]),2)
@@ -166,8 +162,7 @@ def J(mode_mon):
         objective_functions=J,
         objective_arguments=obj_list,
         design_regions=[matgrid_region],
-        frequencies=frequencies,
-        decimation_factor=10)
+        frequencies=frequencies)
 
     f, dJ_du = opt([design_params])
 
@@ -213,7 +208,7 @@ def test_adjoint_solver_DFT_fields(self):
             adj_scale = (dp[None,:]@adjsol_grad).flatten()
             fd_grad = S12_perturbed-S12_unperturbed
             print("Directional derivative -- adjoint solver: {}, FD: {}".format(adj_scale,fd_grad))
-            tol = 0.04 if mp.is_single_precision() else 0.005
+            tol = 0.0461 if mp.is_single_precision() else 0.005
             self.assertClose(adj_scale,fd_grad,epsilon=tol)
 
 
@@ -267,14 +262,14 @@ def test_gradient_backpropagation(self):
                 bp_adjsol_grad = tensor_jacobian_product(mapping,0)(p,filter_radius,eta,beta,adjsol_grad)
 
             ## compute unperturbed S12
-            S12_unperturbed = forward_simulation(mapped_p, MonitorObject.EIGENMODE,frequencies)
+            S12_unperturbed = forward_simulation(mapped_p,MonitorObject.EIGENMODE,frequencies)
 
             ## compare objective results
             print("S12 -- adjoint solver: {}, traditional simulation: {}".format(adjsol_obj,S12_unperturbed))
             self.assertClose(adjsol_obj,S12_unperturbed,epsilon=1e-6)
 
             ## compute perturbed S12
-            S12_perturbed = forward_simulation(mapping(p+dp,filter_radius,eta,beta), MonitorObject.EIGENMODE,frequencies)
+            S12_perturbed = forward_simulation(mapping(p+dp,filter_radius,eta,beta),MonitorObject.EIGENMODE,frequencies)
 
             if bp_adjsol_grad.ndim < 2:
                 bp_adjsol_grad = np.expand_dims(bp_adjsol_grad,axis=1)
diff --git a/src/dft.cpp b/src/dft.cpp
index 1c7e438af..efaf130b0 100644
--- a/src/dft.cpp
+++ b/src/dft.cpp
@@ -176,13 +176,12 @@ dft_chunk *fields::add_dft(component c, const volume &where, const double *freq,
   data.vc = vc;
 
   if (decimation_factor == 0) {
-    double tol = 1e-7;
     double src_freq_max = 0;
     for (src_time *s = sources; s; s = s->next) {
-      if (s->get_fwidth(tol) == 0)
+      if (s->get_fwidth() == 0)
         decimation_factor = 1;
       else
-        src_freq_max = std::max(src_freq_max, std::abs(s->frequency().real())+0.5*s->get_fwidth(tol));
+        src_freq_max = std::max(src_freq_max, std::abs(s->frequency().real())+0.5*s->get_fwidth());
     }
     double freq_max = 0;
     for (size_t i = 0; i < Nfreq; ++i)
@@ -1376,4 +1375,4 @@ void fields::get_mode_mode_overlap(void *mode1_data, void *mode2_data, dft_flux
   get_overlap(mode1_data, mode2_data, flux, 0, overlaps);
 }
 
-} // namespace meep
\ No newline at end of file
+} // namespace meep
diff --git a/src/meep.hpp b/src/meep.hpp
index b8b5aaa07..45d45bb7c 100644
--- a/src/meep.hpp
+++ b/src/meep.hpp
@@ -988,7 +988,8 @@ class src_time {
     return 1;
   }
   virtual std::complex<double> frequency() const { return 0.0; }
-  virtual double get_fwidth(double tol) const { (void)tol; return 0.0; }
+  virtual double get_fwidth() const { return 0.0; }
+  virtual void set_fwidth(double fw) { (void)fw; }
   virtual void set_frequency(std::complex<double> f) { (void)f; }
 
 private:
@@ -1010,12 +1011,13 @@ class gaussian_src_time : public src_time {
   virtual src_time *clone() const { return new gaussian_src_time(*this); }
   virtual bool is_equal(const src_time &t) const;
   virtual std::complex<double> frequency() const { return freq; }
-  virtual double get_fwidth(double tol) const;
+  virtual double get_fwidth() const { return fwidth; };
+  virtual void set_fwidth(double fw) { fwidth = fw; };
   virtual void set_frequency(std::complex<double> f) { freq = real(f); }
   std::complex<double> fourier_transform(const double f);
 
 private:
-  double freq, width, peak_time, cutoff;
+  double freq, fwidth, width, peak_time, cutoff;
 };
 
 // Continuous (CW) source with (optional) slow turn-on and/or turn-off.
@@ -1031,7 +1033,7 @@ class continuous_src_time : public src_time {
   virtual src_time *clone() const { return new continuous_src_time(*this); }
   virtual bool is_equal(const src_time &t) const;
   virtual std::complex<double> frequency() const { return freq; }
-  virtual double get_fwidth(double tol) const { (void)tol; return 0.0; };
+  virtual double get_fwidth() const { return 0.0; };
   virtual void set_frequency(std::complex<double> f) { freq = f; }
 
 private:
@@ -1043,8 +1045,8 @@ class continuous_src_time : public src_time {
 class custom_src_time : public src_time {
 public:
   custom_src_time(std::complex<double> (*func)(double t, void *), void *data, double st = -infinity,
-                  double et = infinity, std::complex<double> f = 0)
-      : func(func), data(data), freq(f), start_time(float(st)), end_time(float(et)) {}
+                  double et = infinity, std::complex<double> f = 0, double fw = 0)
+      : func(func), data(data), freq(f), start_time(float(st)), end_time(float(et)), fwidth(fw) {}
   virtual ~custom_src_time() {}
 
   virtual std::complex<double> current(double time, double dt) const {
@@ -1064,14 +1066,15 @@ class custom_src_time : public src_time {
   virtual src_time *clone() const { return new custom_src_time(*this); }
   virtual bool is_equal(const src_time &t) const;
   virtual std::complex<double> frequency() const { return freq; }
-  virtual double get_fwidth(double tol) const { (void)tol; return 0.0; };
   virtual void set_frequency(std::complex<double> f) { freq = f; }
+  virtual double get_fwidth() const { return fwidth; };
+  virtual void set_fwidth(double fw) { fwidth = fw; }
 
 private:
   std::complex<double> (*func)(double t, void *);
   void *data;
   std::complex<double> freq;
-  double start_time, end_time;
+  double start_time, end_time, fwidth;
 };
 
 class monitor_point {
diff --git a/src/sources.cpp b/src/sources.cpp
index 3c814bd47..5476dc1e6 100644
--- a/src/sources.cpp
+++ b/src/sources.cpp
@@ -60,11 +60,20 @@ double src_time::last_time_max(double after) {
     return after;
 }
 
-gaussian_src_time::gaussian_src_time(double f, double fwidth, double s) {
+// bandwidth (in frequency units, not angular frequency) of the
+// continuous Fourier transform of the Gaussian source function
+// when it has decayed by a tolerance tol below its peak value
+static double gaussian_bandwidth(double width) {
+  double tol = 1e-7;
+  return sqrt(-2.0 * log(tol)) / (width * pi);
+}
+
+gaussian_src_time::gaussian_src_time(double f, double my_fwidth, double s) {
   freq = f;
-  width = 1.0 / fwidth;
+  width = 1.0 / my_fwidth;
   peak_time = width * s;
   cutoff = width * s;
+  fwidth = gaussian_bandwidth(width);
 
   // this is to make last_source_time as small as possible
   while (exp(-cutoff * cutoff / (2 * width * width)) < 1e-100)
@@ -77,6 +86,7 @@ gaussian_src_time::gaussian_src_time(double f, double w, double st, double et) {
   width = w;
   peak_time = 0.5 * (st + et);
   cutoff = (et - st) * 0.5;
+  fwidth = gaussian_bandwidth(width);
 
   // this is to make last_source_time as small as possible
   while (exp(-cutoff * cutoff / (2 * width * width)) < 1e-100)
@@ -105,13 +115,6 @@ std::complex<double> gaussian_src_time::fourier_transform(const double f) {
   return width * polar(1.0, omega * peak_time) * exp(-0.5 * delta * delta);
 }
 
-// bandwidth (in frequency units, not angular frequency) of the
-// continuous Fourier transform of the Gaussian source function
-// when it has decayed by a tolerance tol below its peak value
-double gaussian_src_time::get_fwidth(double tol) const {
-  return sqrt(-2.0 * log(tol)) / (width * pi);
-}
-
 bool gaussian_src_time::is_equal(const src_time &t) const {
   const gaussian_src_time *tp = dynamic_cast<const gaussian_src_time *>(&t);
   if (tp)
@@ -154,7 +157,7 @@ bool custom_src_time::is_equal(const src_time &t) const {
   const custom_src_time *tp = dynamic_cast<const custom_src_time *>(&t);
   if (tp)
     return (tp->start_time == start_time && tp->end_time == end_time && tp->func == func &&
-            tp->data == data && tp->freq == freq);
+            tp->data == data && tp->freq == freq && tp->fwidth == fwidth);
   else
     return 0;
 }