From 25fb4444d030694078bc9dbb6edfb101024bd6d2 Mon Sep 17 00:00:00 2001 From: Eric Teunis de Boone Date: Mon, 14 Nov 2022 11:41:35 +0100 Subject: [PATCH] Continuous DTFT: use derivative == 0 --- fourier/06_direct_fourier.py | 39 ++++++++++++++++++++---------------- 1 file changed, 22 insertions(+), 17 deletions(-) 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):