Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

fix bugs in unwrap function #16977

Closed
wants to merge 1 commit into from
Closed

fix bugs in unwrap function #16977

wants to merge 1 commit into from

Conversation

tatsuoka999
Copy link

Fixed the following two bugs.
Bug1: Returned wrong result when the parameter "discont" is not pi.
Bug2: Wrapped when the diff value is a little larger than pi or discont.

Changes:
L1540, L1541 for Bug1.
L1543 for Bug2.

Example for Bug2:
Fig1

@eric-wieser
Copy link
Member

eric-wieser commented Jul 30, 2020

See also #14877, where some questions were raised about the correctness of this function.

@tatsuoka999
Copy link
Author

tatsuoka999 commented Jul 31, 2020

See also #14877, where some questions were raised about the correctness of this function.

@eric-wieser, thank you for your prompt reply.
I read #14877.
However, the current source would not be correct. I show an example as below.

The following code draws two normalized phase response, one is calculated in radian and another is in degree.
The result in radian is correct, it is not necessary to wrap 2*pi at the point shifted just pi though.
The result in degree is not correct. The phase is changed 180 degree at the point shifted pi.

Code (Python ver 3.7.7):

# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
import numpy as np
from scipy import signal

deg_180 = 180

# HPF frequency response
Fs = 5000
Fc = 20
win = 'hamming'

numtaps = 501
res = 65536

b = signal.firwin(numtaps, Fc, window=win, fs=Fs, pass_zero=False)
w, h = signal.freqz(b, fs=Fs, worN=res)

# Radian
p_wrap = np.angle(h)
p_unwrap = np.unwrap(p_wrap)

p_wrap_n = p_wrap / (2 * np.pi)
p_unwrap_n = p_unwrap / (2 * np.pi)

plt.figure()
plt.plot(w, p_wrap_n, label="Wrapped")
plt.plot(w, p_unwrap_n, label="Unwrapped")
plt.axis([0, 100, -6, 6])
plt.ylabel('Normalized phase [*2*PI rad]')
plt.xlabel('Frequency [Hz]')
plt.legend(loc='best', framealpha=.5, numpoints=1)
plt.title('Radian')
plt.show()

# Degree
p_wrap_deg = np.angle(h, deg=True)
p_unwrap_deg = np.unwrap(p_wrap_deg, discont=deg_180)

p_wrap_deg_n = p_wrap_deg / (2 * deg_180)
p_unwrap_deg_n = p_unwrap_deg / (2 * deg_180)

plt.figure()
plt.plot(w, p_wrap_deg_n, label="Wrapped")
plt.plot(w, p_unwrap_deg_n, label="Unwrapped")
plt.axis([0, 100, -6, 6])
plt.ylabel('Normalized phase [*2*180 deg]')
plt.xlabel('Frequency [Hz]')
plt.legend(loc='best', framealpha=.5, numpoints=1)
plt.title('Degree')
plt.show()

Result:
radian
degree

The correct result is obtained in degree by using the pull request code.

Result:
pullreq

@scimax also mentioned about the problem in the case discont=180 in #14877.
My pull request may not be good but could you please confirm the above?

@tatsuoka999
Copy link
Author

I withdraw about L1543. It was wrong. I am sorry to bother you.
Still, regarding L1540 and L1541, it looks correct.

@scimax
Copy link
Contributor

scimax commented Aug 1, 2020

I have to say I took a quick look yesterday and was very confused about the factor 1.5, but as I had no time to think this through I didn't want to ask yet.

@@ -1537,10 +1537,10 @@ def unwrap(p, discont=pi, axis=-1):
slice1 = [slice(None, None)]*nd # full slices
slice1[axis] = slice(1, None)
slice1 = tuple(slice1)
ddmod = mod(dd + pi, 2*pi) - pi
_nx.copyto(ddmod, pi, where=(ddmod == -pi) & (dd > 0))
ddmod = mod(dd + discont, 2 * discont) - discont
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This 2*pi change is definitely incorrect - the docs refer specifically to comparing discont to pi and 2*pi, but after this change the function does not mention pi at all.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, absolutely true, the meaning of discont was misused here. But the logic is correct. I will correct this. I would propose a name like ´interval_sizesincemin_valandmax_val` from #14877 don't make sense at all. It only corrects differences, so it will do the same correction for

unwrap([ 1, 2, 3, 4, 5, 6, 1, 2, 3], min_val=1, max_val=7)

as for

unwrap([ 1, 2, 3, 4, 5, 6, 1, 2, 3], min_val=0, max_val=6)

. Only the difference max_val - min_val matters.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

since min_valandmax_val from #14877 don't make sense at all.

I don't agree with that - they makes sense to me, and should satisfy unwrap(arr, min=-a, max=a) + b== unwrap(arr + b, min=-a + b, max=a + b)

@scimax
Copy link
Contributor

scimax commented Aug 1, 2020

I compared it to my solution, and, from a semantics point of view I came up with the same solution. But I couldn't solve the issue with integers returning floats, which are mentioned in #14877 .
And I think you misused the parameter discont here. It's purpose was to restrict the size of jumps to be corrected while still dealing with an interval size of 2 pi. But I would like to test your examples.

@eric-wieser , the tests you mentioned in #14877 which failed, are working with this. It's just a mistake in L1540.

@scimax
Copy link
Contributor

scimax commented Aug 1, 2020

I just realized, I can't propose corrections. Neither here nor in #14877 . Shall I open another PR for that?

@tatsuoka999
Copy link
Author

I misunderstood that unwrap() could also be used in degree by changing discont to 180.
Now, I know that this function is only for radian.
I apologize for taking your time and thank you all for your polite answer. @eric-wieser @scimax
This pull request will be deleted.

Also, I confirmed #16987.
This solves what I wanted to do. I agree with this change.

@tatsuoka999 tatsuoka999 closed this Aug 1, 2020
@eric-wieser
Copy link
Member

Thanks for the thorough summary @tatsuoka999, reports with graphs make a nice change :)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

3 participants