Find the relation between two noisy periodic signals of equal length [Thank you, @leftaroundabout for your helpful comments to my related question. I'm trying to ask a less confusing question here.]

$0\leq\phi\leq 2\pi$ is a parameter that theoretically describes a system. $\alpha=\sin(\phi)$ and $\alpha' = \cos(\phi)$ are observables.

In a measurement, $\phi(t)$ is a continuous periodic function of the time $t$. $\alpha$ can be directly measured and $A_i$ is a noisy over-sampled digital representation of such a measurement. $\beta$ is a different observable, the relation between $\beta$ and $\phi$ shall be studied. $B_i$ is a sampled synchronously with $A_i$ and noisy.

Three questions are interesting and the answer shall be found from $A_i$ and $B_i$:


What is $\delta = \beta(\phi=-\pi/2) - \beta(\phi=+\pi/2)$? This the difference of the $\beta$ values corresponding to the extrema of $\alpha$.
Is $\beta = a * \alpha + c$ a valid model and what is the coefficient $a$?
Can we give a confidence interval for $a$ or $\delta$?


Inspection of the data shows that a better model for $B_i$ is $\beta' = a * \alpha + a' * \alpha' + c$. Does this affect the recipe to find $\delta$ or $a$?

A third alternative model is $\beta'' = f(\phi) + c$. How to reliably find $\delta$?

The nature of the noise is unknown. The noise spectrum has both frequencies less then the modulation frequency of $\phi$ and higher frequencies. A Fourier filter (e.g. lowpass) renders $A_i$ and $B_i$ reasonably smooth.
All right, we can cope with this one.

So first, we should like to have some knowledge about the relation of  $i$ and $\phi$,
$$
  \phi\colon\ \mathbb{Z} \to S^1,\qquad i\ \mapsto\ \phi_0+\frac1{\nu_{\!s}\cdot\omega_0}i
  =: \varphi + \lambda i\quad {}_{^{(\operatorname{mod}2\pi)}}
$$
Where we can essentially forget about the $\operatorname{mod}2\pi$, just view it as $S^1=\mathbb{R}/\!{\sim}$.
There are multiple ways to determine the most likely values of $\varphi$ and $\lambda$; what I'd do is a least-square fit in those variables and $A_M$, perhaps also an offset $A_O$ such that
$$
  \zeta := \sum_i \bigl(A_{\!M}\cdot\sin\,\phi(i) + A_O - A_i\bigr)^2
$$
is minimized. A fitting algorithm like this can be somewhat slowish if you want it to work reliable, I don't know if that's a problem.

Once we have this relation $\phi(i)$, we can just read out the $B_i$ for which $i$ is close to $\tfrac\pi2$ resp. close to $-\tfrac\pi2$. We can simply average over all next-neighbours to these real values, this way we're quite immune from nonlinear response, but rather more subject to noise than with using some kind of interpolation. You want to do it this way if you have a very large amount of data, spanning many periods of $\phi$. If you don't, better use some interpolation between the $B_i$ with $\phi(i)\approx\pm\tfrac\pi2$. I'd try a Friedrichs function as convolution kernel.