diff --git a/NOTES.txt b/NOTES.txt
index dead55b..30ad49b 100644
--- a/NOTES.txt
+++ b/NOTES.txt
@@ -8,8 +8,9 @@ Valid sample rates
 
 Sample rates between 300001 Hz and 900000 Hz (inclusive) are not supported.
 They cause an invalid configuration of the RTL chip.
-  rsamp_ratio = 28.8 MHz * 2**22 / sample_rate
-  If bit 27 and bit 28 of rsamp_ratio are different, the RTL chip malfunctions.
+
+rsamp_ratio = 28.8 MHz * 2**22 / sample_rate
+If bit 27 and bit 28 of rsamp_ratio are different, the RTL chip malfunctions.
 
 
 Behaviour of RTL and Elonics tuner
@@ -54,6 +55,36 @@ Elonics IF filters: matched to sample rate (note this may not be optimal)
 RTL AGC mode off
 
 
+Effect of IF signal filtering
+-----------------------------
+
+Carson bandwidth rule:
+  IF_half_bandwidth = peak_freq_devation + modulating_freq
+
+In case of broadcast FM, this is
+  75 kHz + 53 kHz = 128 kHz  (worst case)
+  19 kHz + 53 kHz =  72 kHz  (typical case)
+
+Simulations of IF filtering show:
+ * narrow IF filter reduces noise in the baseband
+ * narrow IF filter causes gain roll-off for high modulating frequencies
+ * narrow IF filter causes harmonic distortion at high modulating deviation
+
+IF filter with 100 kHz half-bandwidth:
+ * baseband gain >= -1 dB up to 75 kHz
+ * less than 0.1% distortion of modulating signal at 19 kHz peak deviation
+ * ~ 2% distortion of modulating signal at 75 kHz peak devation
+
+IF filter with 75 kHz half-bandwidth:
+ * baseband gain ~ -3 dB at 60 kHz, ~ -8 dB at 75 kHz
+ * ~ 1% distortion of modulating signal at 19 kHz peak deviation
+
+Optimal IF bandwidth is probably somewhere between 75 and 100 kHz, with
+roll-off not too steep.
+Weak stations benefit from a narrow IF filter to reduce noise.
+Strong stations benefit from a wider IF filter to reduce harmonics.
+
+
 Effect of settings on baseband SNR
 ----------------------------------
 
diff --git a/TODO.txt b/TODO.txt
index 1d56db0..d9e0683 100644
--- a/TODO.txt
+++ b/TODO.txt
@@ -1,4 +1,4 @@
-* (experiment) make nice plot of baseband distortion due to IF filtering
+* (experiment) consider reducing IF filter bandwidth to ~ 80 kHz
 * (experiment) consider downsampling IF signal before FM detection
 * (experiment) measure effect of IF gain on baseband SNR
 * (experiment) measure effect of IF gain linearity on baseband SNR
@@ -6,7 +6,6 @@
 * (experiment) try if RTL AGC mode improves FM decoding
 
 * (feature) support 'M' 'k' suffixes for sample rates and tuning frequency
-* (feature) implement off-line FM decoder in Python for experimentation
 * (feature) implement stereo pilot pulse-per-second
 * (speedup) maybe replace high-order FIR downsampling filter with 2nd order butterworth followed by lower order FIR filter
 * figure out why we sometimes lose stereo lock
diff --git a/pyfm.py b/pyfm.py
index dfa2069..8ac305c 100644
--- a/pyfm.py
+++ b/pyfm.py
@@ -295,7 +295,7 @@ def pilotLevel(d, fs, freqshift, nfft=None, bw=150.0e3):
     return (p19db, guarddb, guarddb - p19db)
 
 
-def modulateAndReconstruct(sigfreq, sigampl, nsampl, fs, noisebw=None, ifbw=None, ifnoise=0):
+def modulateAndReconstruct(sigfreq, sigampl, nsampl, fs, noisebw=None, ifbw=None, ifnoise=0, ifdownsamp=1):
     """Create a pure sine wave, modulate to FM, add noise, filter, demodulate.
 
     sigfreq     :: frequency of sine wave in Hz
@@ -305,6 +305,7 @@ def modulateAndReconstruct(sigfreq, sigampl, nsampl, fs, noisebw=None, ifbw=None
     noisebw     :: calculate noise after demodulation over this bandwidth
     ifbw        :: IF filter bandwidth in Hz, or None for no filtering
     ifnoise     :: IF noise level
+    ifdownsamp  :: downsample factor before demodulation
 
     Return (ampl, phase, noise)
     where ampl  is the amplitude of the reconstructed sine wave (~ sigampl)
@@ -325,18 +326,24 @@ def modulateAndReconstruct(sigfreq, sigampl, nsampl, fs, noisebw=None, ifbw=None
 
     # Filter IF.
     if ifbw is not None:
-        b  = scipy.signal.firwin(61, 2.0 * ifbw / fs, window='nuttall')
+        b  = scipy.signal.firwin(101, 2.0 * ifbw / fs, window='nuttall')
         fm = scipy.signal.lfilter(b, 1, fm)
         fm = fm[61:]
 
+    # Downsample IF.
+    fs1 = fs
+    if ifdownsamp != 1:
+        fm = fm[::ifdownsamp]
+        fs1 = fs / ifdownsamp
+
     # Demodulate.
-    sig1 = quadratureDetector(fm, fs=fs)
+    sig1 = quadratureDetector(fm, fs=fs1)
 
     # Fit original sine wave.
     k = len(sig1)
     m = numpy.zeros((k, 3))
-    m[:,0] = numpy.sin(2*numpy.pi*sigfreq/fs * (numpy.arange(k) + nsampl - k))
-    m[:,1] = numpy.cos(2*numpy.pi*sigfreq/fs * (numpy.arange(k) + nsampl - k))
+    m[:,0] = numpy.sin(2*numpy.pi*sigfreq/fs1 * (numpy.arange(k) + nsampl - k))
+    m[:,1] = numpy.cos(2*numpy.pi*sigfreq/fs1 * (numpy.arange(k) + nsampl - k))
     m[:,2] = 1
     fit = numpy.linalg.lstsq(m, sig1)
     csin, ccos, coffset = fit[0]
@@ -350,7 +357,7 @@ def modulateAndReconstruct(sigfreq, sigampl, nsampl, fs, noisebw=None, ifbw=None
     res1   = sig1 - m[:,0] * csin - m[:,1] * ccos
 
     if noisebw is not None:
-        b  = scipy.signal.firwin(61, 2.0 * noisebw / fs, window='nuttall')
+        b  = scipy.signal.firwin(101, 2.0 * noisebw / fs1, window='nuttall')
         res1 = scipy.signal.lfilter(b, 1, res1)
 
     noise1 = numpy.sqrt(numpy.mean(res1 ** 2))