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

  • $\begingroup$ 1 - Please post your "Bode plot". 2 - You can filter to improve the quality of the measurement, but you should use the "filtfilt" method to remove the phase shift caused by the filter. 3 - You probably need to unwrap the phase so that your plot looks more like a Bode plot. $\endgroup$
    – Ben
    Mar 23, 2021 at 16:58
  • $\begingroup$ You show the plots of the input and output signal but you didn’t plot your actual Bode plot, can you include that? To compute the phase, hard limit the input and output and then XOR the result (multiply the two square waves). Filter this and the result will be linearly proportional to the phase shift crossing 0 at 90 degrees. This will take advantage of the entire sig to give you an avg result. Be careful if you are filtering the input waveforms as the filter will introduce its own phase shift, so you need to use the exact same filter— best to just do this on the product output of the two $\endgroup$ Mar 23, 2021 at 17:01
  • $\begingroup$ The phase detector approach I gave you is useful over a range of 0 to 180 degrees after which it repeats (aliases) $\endgroup$ Mar 23, 2021 at 17:03
  • $\begingroup$ @Ben I edited my question to include the phase plot. Thanks for the suggestion, currently I'm only using a lowpass filter on the output. If I were to use filtfilt instead, would I again just use it on the output? Moreover, what do you mean by unwrap the phase? Thanks! $\endgroup$
    – migeve3630
    Mar 23, 2021 at 17:31
  • $\begingroup$ @DanBoschen Thanks for the suggestion! I edited my question to include the phase diagram. For XORing the signals, do I square the input and output signals, multiply them together and plot that against time? What exactly would I be looking for on the plot? $\endgroup$
    – migeve3630
    Mar 23, 2021 at 17:33

1 Answer 1


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.

  • $\begingroup$ I really like this approach! Just to make sure I'm understanding this correctly, I should calculate $x$ along the same time domain that I use for $y$, and then point-wise multiply $x$ and $y$ together. I should do this many times and take the mean of all the $z$'s. Once I've done this, I will be left with $z_{mean}$ which is still a vector. How do I go from here to obtaining the angle? Do I take the average of the $z_{mean}$ vector, then compute the argument? $\endgroup$
    – migeve3630
    Mar 23, 2021 at 18:02
  • $\begingroup$ $z_{mean}$ is a complex number, so the angle of the complex number will be the phase shift. $\endgroup$
    – Ben
    Mar 23, 2021 at 18:03
  • $\begingroup$ What I'm confused about is the input signal and output signal will be a vector covering some time domain (ie. 0 to 10 seconds). So wouldn't $z$ also be a vector? $\endgroup$
    – migeve3630
    Mar 23, 2021 at 18:05
  • $\begingroup$ Yep, but the mean of z(t) over 10 seconds, $z_{mean}$ will be a complex number. $\endgroup$
    – Ben
    Mar 23, 2021 at 18:06
  • $\begingroup$ Ah I see my misunderstanding, thank you! Approximately how many cycles do you think would be appropriate? A concern I have is that when I'm dealing with really small frequencies, such as $0.0001$, the period is very large. Running the system for even a couple cycles at that frequency can take a very long time. $\endgroup$
    – migeve3630
    Mar 23, 2021 at 18:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.