Harmonic analysis

[swhite2401]

The harmonic_analysis module is rewritten and cleaned up. I have done test to evaluate the effect windowing and padding for the 2 methods and came to the conclusion that they are not useful for the interpolated FFT case. A warning is now issue when users try to activate windowing or padding using interpolated fft.

The method name is now interp_fft to be more consistent with what it really does, the old keyword laskar is kept for backward compatibility.

Returning only the tune was found to be inappropriate for tbt analysis, the amplitudes and phases are now returned in addition to the tune in an additional function get_main_harmonic and the existing get_spectrum_harmonic. No change to other internal AT function, except for the help, only get_spectrum_harmonic is not backward compatible.

[swhite2401]

I don;t understand why tests are failing

[swhite2401]

I started this in the hope of improving the resolution after the discussion in #556 but the windowing did not work very well for the interpolated FFT. We may want to add the NAFF algorithm as a 3rd method in the future.

[swhite2401]

Looking only at the first harmonic for get_tune (which is what should be done) solves the test issue

[swhite2401]

Looking only at the first harmonic for get_tune (which is what should be done) solves the test issue

In fact using 2 harmonics is safer for real input tbt input data

[lfarv]

We may want to add the NAFF algorithm as a 3rd method in the future.

There is a C implementation of NAFF already in the repository, it was used in Matlab, though apparently it's not proposed anymore. It should be fast enough, however it requires full testing.

Apart from that, is this one ready and tested ?

[swhite2401]

I have tested it and believe it is ready, but it would be good if @oscarxblanco, @lcarver and @carmignani. I particular @carmignani complained about the lost information on the phase form the previous version so I would like to know whether this is satisfactory, if possible I would not want to come back to it

[swhite2401]

Wasn't there a license issue with the NAFF version present in the direction? I can't remember. If not we can use it in fact

[oscarxblanco]

@swhite2401 the licence was GPL-3.0 and should be ok. The issue was speed. The pyNAFF implementation is x20 slower and it is not longer maintained, therefore, some numpy calls are already deprecated.

[oscarxblanco]

@swhite2401

It seems you removed the parameter argument num_harmonics from get_tunes_harmonic.

  File "/home/sources/physmach/blanco-garcia/codes/pyenv/venv_test_harmonic/lib/python3.10/site-packages/at/physics/frequency_maps.py", line 230, in fmap_parallel_track
    xfreqfirst = get_tunes_harmonic(xfirst, num_harmonics=1)
TypeError: get_tunes_harmonic() got an unexpected keyword argument 'num_harmonics'

Is there a replacement ?

[lfarv]

@swhite2401:

Wasn't there a license issue with the NAFF version present in the direction? I can't remember. If not we can use it in fact

@lnadolski ensured that there was no licensing problem with this implementation.

@oscarxblanco: We are talking here of the NAFF implementation in C which is in <at>/atmat/atphysics/nafflib. It has never been compiled for PyAT, to my knowledge, and is very fast. And it has nothing to do with PyNAFF. However, it has not been maintained for long and it would be nice to have an expert's advice. Maybe @lnadolski ?

[swhite2401]

It seems you removed the parameter argument num_harmonics from get_tunes_harmonic.

Yes I did, since the tune is supposed to be the main harmonic I thought this argument was useless. I case many lines are needed one should use the function get_spectrum.

If you think this is an issue it is very easy to restore and set the default to 1

[swhite2401]

It is true it also breaks backward compatibility... I will restore it with a default to 1 it is safer

[oscarxblanco]

I removed the argument num_harmonics=1 in my local copy of frequency_maps.py and the calculation of a frequency map finishes without errors. The resulting diffusion plot appears to be the same wrt to a previous calculation (I attach an image).

With respect to precision, I tested get_tunes_harmonic with a sin function and tried to recover the frequency while varying the number of samples (no noise added, just numerical precision). I attach a figure, but, in general I see a similar of better result for larger number of samples with respect to what was discussed in #556 .

I think one could safely remove num_harmonics=1 from lines 230, 235, 240 and 245 in the file frequency_analysis.py if such modification is included in this PR.

However, it might be still a good idea to keep the backwards compatibility. It might be possible that num_harmonics is used in some other part of the code or related codes. A warning could be shown in case it is decided to remove it on another release.

[swhite2401]

Thanks @oscarxblanco , I have made the requested changes and restored backward compatibility

[swhite2401]

One thing is odd, though. The phase calculated with standard FFT is wrong and I can't figure out why. Any idea?

[oscarxblanco]

Thanks @oscarxblanco , I have made the requested changes and restored backward compatibility

I did a test and it works as expected.

[oscarxblanco]

One thing is odd, though. The phase calculated with standard FFT is wrong and I can't figure out why. Any idea?

Could you provide an example showing inputs, outputs and call to the function ?

[swhite2401]

Here is the code:

import at
import numpy
from at.physics import get_main_harmonic

x = numpy.arange(1024)
qx = 0.23
y = numpy.cos(2*numpy.pi*qx*x+numpy.pi)
tune1, a1, p1 = get_main_harmonic(y, method='fft')
tune2, a2, p2 = get_main_harmonic(y, method='interp_fft')
print('')
print(tune1, tune2)
print(a1, a2)
print(p1, p2)

and the output

[0.23046875] [0.23]
[0.33104111] [0.50001535]
[1.6365377] [-3.14148306]

[swhite2401]

While the amplitude can be different I was expecting a phase of ~pi in both case. Increasing the number does not bring the result to the correct value.

[oscarxblanco]

While the amplitude can be different I was expecting a phase of ~pi in both case. Increasing the number does not bring the result to the correct value.

@swhite2401 It seems to be just an offset on the phase. Are you using functions defined on the same interval ?

[swhite2401]

Ah thanks @oscarxblanco! In fact this was just stupid... taking the periodic signal everything is fine. Hanning window is making the signal periodic but also shifts the phase.

See result below for 2000 turns:

All is ok, and in any case on offset is not really detrimental as you generally look at phase differences.

So this is ready to merge for me, any other comments?

[oscarxblanco]

Ah thanks @oscarxblanco! In fact this was just stupid... taking the periodic signal everything is fine. Hanning window is making the signal periodic but also shifts the phase.

See result below for 2000 turns:

All is ok, and in any case on offset is not really detrimental as you generally look at phase differences.

So this is ready to merge for me, any other comments?

I still get a difference in phase. Is there anyway to avoid using hann ? It seems to be false by default ... `

exec(open('test_simon_white_2.py').read()) [0.23046875] [0.23] [0.33104111] [0.50001535] [1.6365377] [-3.14148306] `

[swhite2401]

Ah this is strange... hann is turned off by default. Did you change the number of turn to 2000? To get the correct phase you need to have nturns*qx=integer

[oscarxblanco]

Right, with 2000 points both functions return the same value.

[swhite2401]

Good! Do you agree to merge then? Can you approve the PR or was there anything else than needs to be changed?

[swhite2401]

Perfect I merge!