Creating Bode Plot from Experimental Data I have a blackbox system in which I can input a function and obtain the output signal (in MATLAB). I'm attempting to reconstruct a Bode diagram and have had success with the Magnitude plot, however, I'm having a lot of difficulty with the phase plot. What I've tried to do is collect input and output data for the input $\sin(ft)$ where $f$ is the frequency that I am changing. For each frequency I am determining the phase shift by finding the time difference (I'll denote by $td$) between the input and output signal on a plot and calculating the phase shift using
$$td \div\frac{2\pi}{f}$$
After collecting the data over a range of frequencies I plot this (phase shift on y-axis, frequency on a log scale on the x axis). The plot does not at all look like a phase plot, I don't really see much of a pattern, nor does the scale seem correct for the phase shifts. Below are some examples of data I have collected for different frequencies:
Frequency: $f=0.05$: https://ibb.co/swCGLw8
Frequency: $f=10$: https://ibb.co/jL8JhTv
Frequency: $f=1000$: https://ibb.co/rkZ8WLh
Note that a filter has been applied to try and reduce the noise
Any help would be great!
EDIT: Image of the phase plot: https://ibb.co/K2LcHqJ
Assuming you have an input signal $ u = A cos(2\pi ft) $ and you measure an output signal
$y = B cos(2\pi ft + \theta) $
The $\frac{B}{A}$ ratio is the gain and $\theta$ is the phase shift for frequency $f$.
You could simply demodulate $y$ by multiplying by $x = cos(2\pi ft) - jsin(2\pi ft)$
$$ z(t) = y*x = \frac{B}{2}(cos(4\pi ft) + cos(\theta) + j(sin(4\pi ft) + sin(\theta))$$
Now simply compute the  average of $z(t)$ over multiple cycles so that the $cos(4\pi ft)$ and $sin(4\pi ft)$ terms will average to 0.
Divide the magnitude of $z_{mean}$ by $A$ to get the gain and compute the angle of $z_{mean}$ to get the phase shift.
The advantage of this method is that by averaging over multiple cycles, you will significantly reduce the noise and improve the accuracy of your phase and gain estimates. Apply this technique for multiple frequencies, you will get your Bode plot.