1
$\begingroup$

I am trying to identify the characteristics of an FIR filter in an ideal case as a sanity check. I observe that there is no strange behaviour, but there is an unexplainable issue in the phase plot of the original and identified filters. To make it simple, I've provided the transfer function coefficients below:

a1 = 1; b1 = [-0.0228015616127155 0.445214656159413 0.776226769351428 0.445214656159413 -0.0228015616127155]; b2 = [-0.0228015616127150 0.445214656159413 0.776226769351428 0.445214656159413 -0.0228015616127150]; G1 = filt(b1,a1); G2 = filt(b2,a1); fmin = 1e-1; fmax = 0.5; Bodeoptions = bodeoptions; Bodeoptions.FreqUnits = 'Hz'; Bodeoptions.Xlim = [fmin,fmax]; Bodeoptions.Ylim = {[-100,10],[0,720]}; Bodeoptions.Grid = 'on'; bode(G1,'-b',G2,':g',Bodeoptions) h = findobj(gcf,'type','line'); set(h,'LineWidth',2); wmin = 2*pi*fmin; H1 = freqresp(G1, wmin); H2 = freqresp(G2, wmin); gain1_dB = 20*log10(abs(H1)); gain2_dB = 20*log10(abs(H2)); phase1_deg = angle(H1)*180/pi; phase2_deg = angle(H2)*180/pi; fprintf('\nAt f = %.4f Hz:\n', fmin); fprintf('G1: Gain = %.16f dB, Phase = %.4f degrees\n', gain1_dB, phase1_deg); fprintf('G2: Gain = %.16f dB, Phase = %.4f degrees\n', gain2_dB, phase2_deg); [mag1, phase1, wout] = bode(G1, wmin); [mag2, phase2, ~] = bode(G2, wmin); gain1_dB = 20*log10(squeeze(mag1)); gain2_dB = 20*log10(squeeze(mag2)); phase1_deg = squeeze(phase1); phase2_deg = squeeze(phase2); fprintf('\nFrom bode():\n'); fprintf('G1: Gain = %.16f dB, Phase = %.4f degrees\n', gain1_dB, phase1_deg); fprintf('G2: Gain = %.16f dB, Phase = %.4f degrees\n', gain2_dB, phase2_deg);

Output:

At f = 0.1000 Hz: G1: Gain = 3.4199354940498718 dB, Phase = -72.0000 degrees G2: Gain = 3.4199354940498772 dB, Phase = -72.0000 degrees From bode(): G1: Gain = 3.4199354940498852 dB, Phase = 648.0000 degrees G2: Gain = 3.4199354940498634 dB, Phase = 288.0000 degrees

G1 is true, G2 is estimation

The difference between b1 and b2 is only at the level of machine precision, yet I observe 360$^\circ$ phase shift in the phase plot. I understand that a 360$^\circ$ phase shift has no physical significance, but I would like to understand why this happens with the bode function in MATLAB. Also, it would be great if you could tell me how to correct this when plotting?

Improve this question
$\endgroup$
3
  • $\begingroup$ Attaching the plot generated by your bode statement would be helpful, as well as the exact gain and phase of the two systems at the lowest frequency being plotted. $\endgroup$ Commented Apr 11 at 20:10
  • 2
    $\begingroup$ Probably what's happening is that Matlab's bode function is unwrapping the phase, so you don't get hard-to-interpret jumps in phase from $2 \pi - \epsilon$ and $0$, and that you're starting at slightly different phases, leading to your phase difference. $\endgroup$ Commented Apr 11 at 20:11
  • $\begingroup$ Thank you, I updated the bode plot, gain and phase margins. Is there a manual way to correct it? Like tweaking the TF? $\endgroup$ Commented Apr 12 at 7:43

1 Answer 1

Reset to default
3
$\begingroup$

but I would like to understand why this happens with the bode function in MATLAB.

That's difficult to answer since you basically asking about the inner workings of a commercial piece of software.

Also, it would be great if you could tell me how to correct this when plotting?

Both answers are correct. So it's not about "correctness" but about "consistent visual experience". This partially depends on personal preference and there is no single right or wrong answer.

Here is what I do:

  1. Find a good "anchor frequency": a frequency that's in the band of interest, has good SNR and doesn't have some crazy features in the magnitude response. Often, that's simply the location of maximum of the magnitude.
  2. Calculate the phase at this frequency in the range of $[-\pi,+\pi]$
  3. Us this is a starting point and then unwrap downwards to DC and unwrap upwards to Nyquist.

In 95% of all cases that gets me what I want but occasionally I will add or subtract integer multiples of $2\pi$ to make it look prettier.

Improve this answer
$\endgroup$
2
  • $\begingroup$ I'm pretty sure that Matlab both has an "unwrap" function (I think by that name, even), and that it will take an amplitude-phase vector and plot it. Dig into the documentation; there should be help somewhere for the sort of manual intervention Hilmar is suggesting. $\endgroup$ Commented Apr 13 at 5:08
  • 1
    $\begingroup$ @TimWescott Matlab has indeed an unrwap function but it always starts unwrapping at DC which isn't a great choice for, say, a high pass filter. $\endgroup$ Commented Apr 13 at 12:03

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.