I'm not sure I totally understand what would be best, but I believe that looking at Dynamic Linear Models, State Space modeling, and the Kalman Filter would be helpful. In R, the package dlm is fairly accessible.
Edit: Describe DLM's simply, eh? One explanation I've seen is that DLMs are regression where your coefficients are allowed to change with time. Doesn't really give me much intuition.
So, here is my I-realize-I-don't-understand-DLMs-well-enough-to-do-this-well-but-you-asked answer, which I hope others will correct as necessary. Let me use a situation that's similar to how they were invented...
Say you were controlling a remotely-piloted vehicle. You could summarize the vehicle's actual state (speed, direction, altitude, fuel, etc, etc) as a vector in a state space. You can't directly observe the actual state, but you do have sensors that observe linear combinations of the actual state, with some (gaussian) noise added in. That is, your sensors are not perfect.
Further, the vehicle has a program that determines how it changes states, as a linear combination of the current state, also with some (gaussian) noise. That is, your controls are not perfect.
Say you want to know the most likely and most accurate path possible for the vehicle -- that is, the states it actually went through: how it actually moved. What you observe is noisy, but you can use a Kalman filter on this system you model and the filter models the variances in the system and outputs the most likely actual states. At each time step, it forecasts what it believes the state will be, then observes what the sensors say, calculates the variances of the different parts, and comes up with a final estimate of the state as a weighted sum of the forecast and the sensor readings. Then the process is repeated for the next time step.
The whole concept is basically a continuous Hidden Markov Model, if you're familiar with those. The model is several equations which involve several matrices that are multiplied together, and you can add columns to the matrix to reflect trend, seasonal, and other types of components of the time series.
The R package dlm makes it particularly easy to define and combine components.
I believe that calculating the coherence would be of use to you.
The coherence at a frequency, $C_{XY}(f)$, between two series, $x(n)$ and $y(n)$, is defined as:
$C_{XY}(f) = \frac{|P_{XY}(f)|^2}{P_{X}(f)P_{Y}(f)}$
Where $P_X(f)$ and $P_Y(f)$ are the power spectrum of the series $x(n)$ and $y(n)$ respectively. $P_{XY}(f)$ is called the cross-spectrum.
Let's say that we calculate the spectrum of $x(n) = x_n$ by first Fourier transforming the windowed data:
$X(f) = \sum_{n=0}^{N-1}h_n x_n e^{-i2\pi fn}$
where $h_n = h(n)$ is the $n^{\text{th}}$ value of the window we are using. This implies that the power spectrum of $x(n)$ is:
$P_{X}(f) = |X(f)|^2$
The cross spectrum is then:
$P_{XY}(f) = |X(f) \cdot Y^{*}(f)|^{1/2}$, where $*$ denotes complex conjugate.
Thus, the coherence is:
$C_{XY}(f) = \frac{|X(f) \cdot Y^{*}(f)|}{|X(f)|^2 |Y(f)|^2}$.
After this, you could calculate the average coherence I suppose - I don't think averaging is the correct approach honestly, I will do some more checking. This is hopefully a start though.
Using a multitaper method for estimating your power spectra is a good idea (Spectrum Estimation and Harmonic Analysis, DJ Thomson, 1982 (paper) or Spectral Analysis for Physical Applications, Pervcival and Walden, 1993 (book)) as you are able to get a much reduced variance for your estimate.
If you decide to go that route, I would suggest checking out: http://www.spectraworks.com/Help/mtmtheory.html for more info.
Best Answer
Try using recurrence plot or perform a Fourier transform.
In the first case, the periodicity is translated to equally spaced lines with the slope 1 (and $R(t,t'+T) = R(t,t')$). This method is more robust, and easier to visualize, than the one below.
Then you can try to detect lines from the plot. E.g.:
Approach 1: Subtract the diagonal and calculate fractal dimension by box counting. The result, which is in the range of $[0,2]$, is 1 for a periodic signal, but not only for it.
Approach 2: Rotate by $45^\circ$ (or rather - take subimage rotated by $45^\circ$) and perform the Singluar Value Decomposition - then $\lambda_0/\sum_i \lambda_i$ (where $\lambda_0$ is the largest singular value) is an easy (but I guess not that robust) measure of periodicity (closer to one, the more periodic).
There are other methods for detection of lines (e.g. Radon transform) but here they seem to be using a canon to kill a fly.
In the second case, for a periodic function its Fourier transform $\tilde{f}$ is a sum of equally spaced Dirac delta peaks: $$\tilde{f}(\omega) = \sum_{k=-\infty}^\infty a_k \delta(\omega - k \omega_0).$$ Of course for every practical data it will be smeared. If you want a quantification then $\int |\tilde{f}(\omega)|^4 d\omega$ may work (or use Shannon/Renyi entropy) - the larger it is, the more periodic is the function.
EDIT: Edited considerably to give numbers characterizing periodicity.