0

Being A and B two cosine data series (same Xvalues, different Yvalues), I want to identify a third cosine function (and its parameters) so that: A(x) = B(x) + C(x) [A is the blue curve, B is the red curve below]

I just want to discover which cosine function (or N cosine functions) I must add to B(x) in order to yield A(x) (or as close as possible).

inserir a descrição da imagem aqui

So far I have managed to read data points for A(x) and B(x) and to use scipy.optimize.curve_fit() to recover cosine functions parameters for each series. 

    def fit_cos1(angles,energy):

        ang = numpy.array(angles)
        energy = numpy.array(energy)
        ff = numpy.fft.fftfreq(len(ang), (ang[1]-ang[0]))   # assume uniform spacing
        Fenergy = abs(numpy.fft.fft(energy))

        guess_freq = abs(ff[numpy.argmax(Fenergy[1:])+1])   # excluding the zero frequency peak, which is related to offset
        guess_amp = numpy.std(energy) * 2.**0.5
        guess_offset = numpy.mean(energy)
        guess = numpy.array([guess_amp, 2.*numpy.pi*guess_freq, 0., guess_offset])

        def cosfunc(t, K, N, Y, c):
            return K*(1 - numpy.cos(N*t + Y)) + c

        popt, pcov = scipy.optimize.curve_fit(cosfunc, ang, energy, p0=guess)
        K, N, Y, c = popt


        fitfunc = lambda t: K*(1 - numpy.cos(N*t + Y)) + c

        return {"amp": round(K,2), "omega": round(N*57.2958,2), "phase": round(Y*57.2958,2), "coeff": round(c,2), "fitfunc": fitfunc, "maxcov": numpy.max(pcov), "rawres": (guess,popt,pcov)} # outputing in degrees


angles = [0,30,60,90,120,150,180,210,240,270,300,330,360]
energyA = [0,3.3908, 10.7669, 14.9809, 11.526, 3.9065, 0.0016, 3.9057, 11.5265, 14.9811, 10.7667, 3.39, 0]
energyB = [4.679575, 3.251112, 0.951731, 0.0, 1.047923, 3.352932, 4.688559, 3.31077, 1.047938, 0.00, 0.951859, 3.25869, 4.679554]

A = fit_cos1(angles, energyA)
print( "CurveA: Amplitude=%(amp)s, Periodicity=%(omega)s, Phase=%(phase)s, V-shift=%(coeff)s, Max. Cov.=%(maxcov)s" % A )

B = fit_cos1(angles, energyB)
print( "CurveB: Amplitude=%(amp)s, Periodicity=%(omega)s, Phase=%(phase)s, V-shift=%(coeff)s, Max. Cov.=%(maxcov)s" % B )

The cosine parameters obtained seems ok and within the error margin I would expect.

CurveA: Amplitude=7.48, Periodicity=2.03, Phase=-5.31, V-shift=0.04, Max. Cov.=0.00130465635936
CurveB: Amplitude=-2.35, Periodicity=1.98, Phase=4.31, V-shift=4.6, Max. Cov.=0.00214051830825

From there, does anyone has an idea on how to find the third cosine function so that A(x) = B(x) + C(x)?

eyllanesc
  • 235,170
  • 19
  • 170
  • 241
mdpoleto
  • 711
  • 2
  • 6
  • 10
  • What exactly is your problem? Are the parameters found not correct? You don't need `fitfunc` and `cosfunc`. Try to use only one of them. For debugging remove the FFT stuff and use some fixed values and see if the fit finds parameters that you started with. From simple to complex. – Joe Sep 06 '19 at 06:05
  • 1
    Why not just subtract one from the other? – Mad Physicist Sep 06 '19 at 12:32
  • Thanks for the suggestions, Joe. I've edited the question in order to make it clearer. I have debugged the code and everything seems to work fine for what I have done (finding the parameters for the functions I do have). My problem is on how to find the third (or N) cosine function that would satisfy A(x) = B(x) + C(x) – mdpoleto Sep 06 '19 at 12:32
  • Interestingly, from the graphs, I'd say you would need an offset and multiplication to go from A to B, not an added cosine function. But I assume that doesn't fit the (physical) model underneath it. – 9769953 Sep 06 '19 at 12:38
  • @MadPhysicist Indeed, that was my first thought, but the main idea is to derive N cosine functions that could solve the equation (A(x) = B(x) + C(x) ; or A(x) = B(x) + C(x) + D(x) and so on..). Is there an option to derive N functions from the residue of A(x) and B(x)? – mdpoleto Sep 06 '19 at 13:56
  • The dct of the difference is exactly that. Keep in mind that this only works if there is no phase shift. If there's a shift, you need a full fft – Mad Physicist Sep 06 '19 at 18:19
  • @MadPhysicist ,I see. Could you elaborate a bit more on how to design a complete FFT (including phase shift) for the difference? – mdpoleto Sep 07 '19 at 22:19
  • Subtract the two functions and take the FFT. It's really that simple. The FFT is the coefficients of the sine and cosine functions of different frequencies that add together to give the input. – Mad Physicist Sep 08 '19 at 04:33
  • @MadPhysicist but is `sin` allowed? If it is only `cos` it'll not work due to symmetry. – mikuszefski Sep 09 '19 at 10:52
  • You seem to be in need of an introductory signal processing textbook, not really programming help. That's out of scope for SO. – Mad Physicist Sep 09 '19 at 12:02

0 Answers0