complex numbers multiplication in double precision

조회 수: 127(최근 30일)
FastCar 2017년 7월 1일
댓글: James Tursa 2017년 7월 3일
I have to multiply couple of complex numbers and then I have to add all the product.
Let's say the two complex numbers (with unitary module) are c1 and c2. What I need is prod=c1 x c2* (but I can substitute c2 with its conjiugate in our example). c1 has phase phi1, c2 has phase phi2. I can do it in chartesian coordinates
c1 x c2* = real(c1) x real (c2) + imag(c1) x imag(c2) + j x (imag(c1) x real (c2) - real(c1) x imag(c2)
(if I have used the right formula)
or I can use
c1 x c2* = cos(phi1-phi2) + j x sin(phi1-phi2)
Mind that the phases are of the order of magnitude of 10^10 radians.
The second way is of course faster, but my question is which one is the more precise given that the computation is done in double precision?
Mind also that I have to sum 10^7 products and I can tell you that there is a difference, a not neglectable difference.
Many thanks in advance


Walter Roberson
Walter Roberson 2017년 7월 1일
편집: Walter Roberson 2017년 7월 1일
On my system, even excluding the time it takes to calculate the angle, the cos/sin version is about 1.8 to 2.4 times slower than the other version. This is to be expected as cos and sin require a many more calculations. There are hardware instructions, yes, but they are slower than plain multiply or addition.
The second version appears to have marginally better numerical properties eps(1) compared to 2*eps(1). That is not really a meaningful difference compared to the rough-off error in the rest of the calculations.
  댓글 수: 2
Walter Roberson
Walter Roberson 2017년 7월 2일
I need the formulae you used.

댓글을 달려면 로그인하십시오.

James Tursa
James Tursa 2017년 7월 2일
편집: James Tursa 2017년 7월 2일
"... Mind that the phases are of the order of magnitude of 10^10 radians. ..."
When you have angles that large, I have to start questioning where they came from and how exactly you are using them in your code. Realize that 10^10 is more than 9 orders of magnitude larger than 2pi, so this will come into play in the "accuracy" you expect for calculations. E.g., a naive example:
>> a = 10^10
a =
>> cos(a)+sin(a)*i
ans =
0.873119622676856 - 0.487506025087511i
>> cos(mod(a,2*pi))+sin(mod(a,2*pi))*i
ans =
0.873119809656583 - 0.487505690208075i
So you can see a difference in the 6th digit. What is going on? It turns out that the trig functions have a very sophisticated method of dealing with the argument. I don't know the exact method that is used behind the scenes, but effectively it treats the IEEE double precision argument as if it has the same absolute precision as a IEEE double precision pi. E.g.,
>> va = vpa(a)
va =
>> vp = vpa(pi)
vp =
>> n = floor(va/(2*vp))
n =
>> mod(a,2*pi)
ans =
>> double(va-n*2*vp)
ans =
>> cos(ans)+sin(ans)*i
ans =
0.873119622676856 - 0.487506025087511i
So you can see by treating the angle as if it had the same absolute precision as a double precision pi value we can do the mod in such a way that we can uncover how cos() and sin() are behaving behind the scenes. It is effectively a different angle that is used compared to the naive mod(a,2*pi) angle result.
Bottom line is that with 10^10 radian angles you need to be very careful that your problem has real meaning at these values, and then be extra careful how you are doing calculations with them to come up with meaningful results. For your phi1-phi2 method, realize that the cos() and sin() functions are going to treat that phi1-phi2 double precision result as if it had absolute precision the same as a double precision pi. This effect may be creeping into your results.

FastCar 2017년 7월 2일
편집: FastCar 2017년 7월 2일
I try to clarify my question: the phase is given by the distance of a satellite from the Earth multiplied by
4 x pi / lambda
The distance is big (around 40000 km) and thus gives a huge phase. So far I used these expressions:
distance = norm ( position vector )
phi = ( 4 x pi x distance ) / lambda
method one
c1= cos(phi1) + j x sin(phi1)
c2= cos(phi2) + j x sin(phi2)
c1 x c2* = cos(phi1) x cos(phi2) + sin(phi1) x sin(phi2) + j x (sin(phi1) x cos (phi2) - cos(phi2) x sin(c2)
method 2
c1 x c2* = cos(phi1-phi2) + j x sin(phi1-phi2)
For what I know every computation in matlab should be done in double precision and without the symbolic packed pi should have only 15 figures (don't blame me if I am wrong). Mind that 15 figures (and thus 14 sure figures) is a good output in my point of view.
I hope I have clarified the question.
I haven't answered to the single question because I think both of you needed the full picture.
  댓글 수: 1
James Tursa
James Tursa 2017년 7월 3일
I don't quite understand your physical application yet, but keep in mind what I have already written in my Answer. Because phi is so large, you have effectively already lost 9 significant digits from your answer. So the c1 x c2* result will only have about 6 meaningful digits of accuracy at most.

댓글을 달려면 로그인하십시오.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by