diff --git a/fourier/06_direct_fourier.py b/fourier/06_direct_fourier.py
index d29c758..8d3d441 100644
--- a/fourier/06_direct_fourier.py
+++ b/fourier/06_direct_fourier.py
@@ -13,7 +13,7 @@ if __name__ == "__main__":
     t = np.linspace(0, 100, n_samples+1)
 
     phi_in = np.linspace(0, 2*np.pi, nphi)
-    test_freqs = f_beacon + np.linspace(-0.1, 0.1, 300+1)
+    test_freqs = f_beacon + 0.1 * np.linspace(-1, 1, 100+1)
 
 
     phi_out = np.zeros( (len(phi_in), len(test_freqs)) )
@@ -80,13 +80,13 @@ if __name__ == "__main__":
         # Single Amplitude / Frequency plot showing frequency fitting
         freq_out = None
         if True:
-            from numpy.polynomial import Polynomial
+            from numpy.polynomial import Polynomial as P
         
             freq_out = np.zeros(len(phi_in))
             amp_cut = 0.8
 
             fig, ax = plt.subplots()
-            ax.set_title("Frequency estimation by parabola fitting.")
+            ax.set_title("Frequency estimation by parabola fitting.\nStars are used for the parabola fit, vertical line is where $\\partial_f = 0 $")
             ax.set_xlabel("Frequency")
             ax.set_ylabel("Amplitude")
         
@@ -100,28 +100,33 @@ if __name__ == "__main__":
                 # filter amplitudes below amp_cut*max_amp
                 valid_idx = amp >= amp_cut*max_amp
         
-                p_fit = Polynomial.fit(test_freqs[valid_idx], amp[valid_idx], 2)
+                p_fit = P.fit(test_freqs[valid_idx], amp[valid_idx], 2)
                 func  = p_fit.convert()
-                
-                tmp_test_freqs = test_freqs[max_amp_idx] + 0.05*np.linspace(-1,1,101, endpoint=True)
-                func_amps = func(tmp_test_freqs)
-                freq_id = np.argmax(func_amps)
 
-                freq_out[j] = tmp_test_freqs[freq_id]
+                # Find frequency of derivative == 0
+                deriv = func.deriv(1)
+                freq = deriv.roots()[0]
+                freq_out[j] = freq
 
-                # plot tmp_test_freqs and freq_id
-                func_amps_idx = func_amps > 0
-                func_amps = func_amps[func_amps_idx]
-                tmp_test_freqs = tmp_test_freqs[func_amps_idx]
 
-                ax.plot(test_freqs, amp, marker='o')
-                l = ax.plot(tmp_test_freqs, func_amps)
+                l = ax.plot(test_freqs, amp, marker='.')
+                ax.plot(test_freqs[valid_idx], amp[valid_idx], marker='*', color=l[0].get_color())
                 ax.axvline(freq_out[j], color=l[0].get_color())
 
+                if True: # plot the fit
+                    tmp_test_freqs = test_freqs[max_amp_idx] + 0.05*np.linspace(-1,1,101, endpoint=True)
+                    func_amps = func(tmp_test_freqs)
+
+                    func_amps_idx = func_amps > 0
+                    func_amps = func_amps[func_amps_idx]
+                    tmp_test_freqs = tmp_test_freqs[func_amps_idx]
+
+                    ax.plot(tmp_test_freqs, func_amps, ls='dotted', color=l[0].get_color())
+
         # Amplitudes figure
-        if True:
+        if not True:
             fig, ax = plt.subplots()
-            ax.set_ylabel("Amplitude")
+            ax.set_ylabels("Amplitude")
             ax.set_xlabel("Frequency")
             if True:
                 for j, amp in enumerate(amp_out.T):