diff --git a/applications/lofar2/model/rtdsp/fourier.py b/applications/lofar2/model/rtdsp/fourier.py index b49c8e79bcfe2e6227061e9ae135961722ad701f..1b319fa9722f67e1ff93a2e45064c080aafb8f8e 100644 --- a/applications/lofar2/model/rtdsp/fourier.py +++ b/applications/lofar2/model/rtdsp/fourier.py @@ -133,6 +133,45 @@ def dtft(coefs, Ndtft=None, zeroCenter=True, fftShift=True): return h, f, HF +def estimate_gain_at_frequency(f, HF, freq): + """Find gain = abs(HF) at frequency in f that is closest to freq. + + Input: + . f, HF: frequency axis and spectrum as from dtft() with fftShift=True + . freq : frequency to look for in f + Return: + . fIndex: index of frequency in f that is closest to freq + . fValue: frequency at fIndex + . fGain : gain = abs(HF) at fIndex + """ + fDiff = np.abs(f - freq) + fIndex = fDiff.argmin() + fValue = f[fIndex] + fGain = np.abs(HF[fIndex]) + return fIndex, fValue, fGain + + +def estimate_frequency_at_gain(f, HF, gain): + """Find positive frequency in f that has gain = abs(HF) closest to gain. + + Input: + . f, HF: frequency axis and spectrum as from dtft() with fftShift=True + . gain : gain to look for in HF + Return: + . fIndex: index of frequency in f that has abs(HF) closest to gain + . fValue: frequency at fIndex + . fGain : gain = abs(HF) at fIndex + """ + pos = len(f[f <= 0]) + HFpos = HF[pos:] + HFgain = np.abs(HFpos) + HFdiff = np.abs(HFgain - gain) + fIndex = pos + HFdiff.argmin() + fValue = f[fIndex] + fGain = np.abs(HF[fIndex]) + return fIndex, fValue, fGain + + ############################################################################### # Radix FFT ###############################################################################