diff --git a/fourier/06_direct_fourier.py b/fourier/06_direct_fourier.py index 8d3d441..e160b1e 100644 --- a/fourier/06_direct_fourier.py +++ b/fourier/06_direct_fourier.py @@ -1,4 +1,12 @@ #!/usr/bin/env python3 +# vim: fdm=indent ts=4 + +__doc__ = \ +""" +Show how the fourier transform can be calculated +in a continuous fashion +""" + if __name__ == "__main__": import numpy as np import matplotlib.pyplot as plt @@ -83,12 +91,14 @@ if __name__ == "__main__": from numpy.polynomial import Polynomial as P freq_out = np.zeros(len(phi_in)) - amp_cut = 0.8 + amp_cut = 0.5 fig, ax = plt.subplots() 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") + + ax.axvline(f_beacon, lw=5, ls=(0,(5,5))) for j, amp in enumerate(amp_out): @@ -98,9 +108,20 @@ if __name__ == "__main__": max_amp_idx = np.argmax(amp) max_amp = amp[max_amp_idx] # filter amplitudes below amp_cut*max_amp - valid_idx = amp >= amp_cut*max_amp - - p_fit = P.fit(test_freqs[valid_idx], amp[valid_idx], 2) + valid_mask = amp >= amp_cut*max_amp + + if True: + # make sure not to use other peaks + lower_mask = valid_mask[0:max_amp_idx] + upper_mask = valid_mask[max_amp_idx:] + + lower_end = np.argmin(lower_mask[::-1]) + upper_end = np.argmin(upper_mask) + + valid_mask[0:(max_amp_idx - lower_end)] = False + valid_mask[(max_amp_idx + upper_end):] = False + + p_fit = P.fit(test_freqs[valid_mask], amp[valid_mask], 2) func = p_fit.convert() # Find frequency of derivative == 0 @@ -110,7 +131,7 @@ if __name__ == "__main__": l = ax.plot(test_freqs, amp, marker='.') - ax.plot(test_freqs[valid_idx], amp[valid_idx], marker='*', color=l[0].get_color()) + ax.plot(test_freqs[valid_mask], amp[valid_mask], marker='*', color=l[0].get_color()) ax.axvline(freq_out[j], color=l[0].get_color()) if True: # plot the fit @@ -126,7 +147,7 @@ if __name__ == "__main__": # Amplitudes figure if not True: fig, ax = plt.subplots() - ax.set_ylabels("Amplitude") + ax.set_ylabel("Amplitude") ax.set_xlabel("Frequency") if True: for j, amp in enumerate(amp_out.T):