Continuous DTFT: only use one peak for fitting

This commit is contained in:
Eric Teunis de Boone 2022-11-14 15:38:32 +01:00
parent 25fb4444d0
commit feccf64293

View file

@ -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):