# Arithmetic operations with complex number

Hello,

I am trying to work with complex number in Igor Pro for the first time. After reading some posts and the manual, I was able to do some computation but the results look strange, it would be great if you can give me a hint.

I have two complex waves, SWave and ScWave, defined as follow: (Iwave and Qwave are real)
make/o/c/n=(numRow) SWave, ScWave, momentWave
SWave = cmplx(Iwave, Qwave)
ScWave = cmplx(Iwave, -Qwave)

Then I try to first raise each wave to some power, SWave to m and ScWave to n, where m,n are from 0 to 4.
for (n=0;n<5;n=n+1)
for (m=0;m<5;m=m+1)
//compute the (S*)^n x (S)^m
momentWave=(ScWave^n)*(SWave^m)
print mean(momentWave)
momentMatOneTrig[n][m]=mean(momentWave)
endfor
endfor

Then it turns out, when m or n is zero, I get "nan" as the result.
Also, the first row and first column are both "nan".
I was expecting at least for the case where m=n=0, it should give me cmplx(1,0).

In the case of m=n=0,
following the debugger, the momentWave are all (1,0), but then mean(momentWave) becomes "nan".

Sorry if I am missing something obvious, thanks for reading.
Just tried, it seems I can get the expected result if I use Wavestats instead of mean.
My guess is mean doesn't work well in computing the average for complex number, would you have a better way than using wavestats (as it may be slow)?
johnweeks wrote:
You can use `MatrixOP theMean = mean(complexWave)`

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

Thanks John,
maybe I made a mistake but the evaluation given by this commend seems to be incorrect.
For example,
I created a wave of length 2e5, and assigned all elements with cmplx(1,0).
Then when I use the command to compute the mean, it gives (nan,nan).

Is there some sorts of limitation over length of wave?
johnweeks wrote:
Using Igor 8? Seems there was a bug, fixed just a few minutes ago. Try the nightly build tomorrow.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

I tried Igor 6 and 7, they were also affected.
I am still on Igor 6 mainly, but moving to Igor 7 lately.
Spending the whole day pinning down the issue, it was about the length of the complex wave.
Somehow as long as it is longer than an unknown value, it won't average correctly, even I am trying to use running average that didn't use the mean/sum or other internal functions.

The only way that works for me for now is going back to wavestats/q:

function/c computeCPXAvg(inputCPXwave)
wave/d/c inputCPXwave

variable numRow = dimsize(inputCPXwave,0)

variable n
variable/c m

//breaker for cmplx number
make/d/o/n=(numRow) inputReal, inputImag

inputReal = real(inputCPXwave)
inputImag = imag(inputCPXwave)

variable avgReal, avgImag
wavestats/q inputReal
avgReal = V_avg
wavestats/q inputImag
avgImag = V_avg

variable/c avgCPXnum = cmplx(avgReal,avgImag)

//  for (n=0;n<numRow;n=n+1)
//      //sumCPXnum = sumCPXnum + inputCPXwave[m]
////        avgCPXnum=(avgCPXnum*(m)+inputCPXwave[m])/(m+1)
////        avgReal=real((avgCPXnum*m+inputCPXwave[n])/(m+1))
////        avgImag=imag(avgCPXnum*m+inputCPXwave[n])/(m+1)
//
//      avgCPXnum=cmplx(avgReal,avgImag)
//      endfor

//variable/c avgCPXnum = sumCPXnum/numRow

return avgCPXnum

end
Sandbo wrote:
Spending the whole day pinning down the issue, it was about the length of the complex wave.
Somehow as long as it is longer than an unknown value, it won't average correctly, even I am trying to use running average that didn't use the mean/sum or other internal functions.

I unfortunately introduced a bug into the MatrixOP computation of the mean that kicks in when the number of points in the wave triggers automatic multithreading. This was fixed for the nightly build of IP8. There are several ways to get around the issue. For example, in the version you are using you can get around this by executing

or use the granularity of the operation to affect only MatrixOP.

AG
Sandbo wrote:
I tried Igor 6 and 7, they were also affected.

The following works correctly for me on Macintosh and Windows with Igor6, Igor7, and Igor8 with the 2018-06-02 nightly build (choose Help->Igor Pro Nightly Builds):
Function TestMean()
Make/O/C/N=(2E5) ctest = cmplx(1,0)

// Prints "(1,0)" in Igor6, Igor7, Igor8
Print/C mean(ctest)

// Prints (1,0) in Igor6, Igor7, and Igor8 with 2018-06-02 nightly build or later
MatrixOp/O theMean = mean(ctest)
Print theMean

KillWaves ctest
End

If you would like to pursue this, try to create a simplified, self-contained example like the one above that shows the problem and post it.
hrodstein wrote:
Sandbo wrote:
I tried Igor 6 and 7, they were also affected.

The following works correctly for me on Macintosh and Windows with Igor6, Igor7, and Igor8 with the 2018-06-02 nightly build (choose Help->Igor Pro Nightly Builds):
Function TestMean()
Make/O/C/N=(2E5) ctest = cmplx(1,0)

// Prints "(1,0)" in Igor6, Igor7, Igor8
Print/C mean(ctest)

// Prints (1,0) in Igor6, Igor7, and Igor8 with 2018-06-02 nightly build or later
MatrixOp/O theMean = mean(ctest)
Print theMean

KillWaves ctest
End

If you would like to pursue this, try to create a simplified, self-contained example like the one above that shows the problem and post it.

Thanks for your code, I tried and it indeed worked, falsifying the length being a cause.
And I tried to reproduce the result using the closest possible code I made:
function deBuggingIgorCmplx(whichcase)
variable whichcase

//some constants
variable h=6.62606979e-34
variable freq=4.2e9
variable BW=2e5
variable Gain=4.7777e6
variable LC = -13.9

//parameters
variable size=4

//container matrix for output
make/D/c/o/n=(size,size) momentMatOneTrig

//make/D temporary containers for I and Q waves, note that they are now double precision
make/D/o/n=2e5 Iwave, Qwave

//To continue, I tried two cases, one with fake gnoise with similar amplitude as my data;
//The second case is my real data which should also be gaussian, acquired from a digitizer.
//The data wave are supposed to be 32-bit floating points.
switch (whichcase)
case 1:
//make fake input data
make/o/n=2e5 Ifake=(gnoise(1,1))*0.001, Qfake=(gnoise(1,1))*0.001
Iwave= Ifake[p]
Qwave= Qfake[p]
break
case 2:
//Get real data
wave testIQdata
//wave assignment
Iwave= testIQdata[p][0]
Qwave=testIQdata[p][1]
break
endswitch

variable m,n,k,i,j

Iwave = Iwave*sqrt(1e-3*10^(LC/10))/sqrt(h*freq*BW*Gain)
Qwave = Qwave*sqrt(1e-3*10^(LC/10))/sqrt(h*freq*BW*Gain)

//We now have the scaled voltage wave I and Q. We can form complex number wave S, as S = I+iQ
make/D/o/c/n=2e5 SWave, ScWave, momentWave
SWave = cmplx(Iwave, Qwave)
ScWave = cmplx(Iwave, -Qwave)

for (n=0;n<size;n=n+1)
for (m=0;m<size;m=m+1)
momentWave=(ScWave^n)*(SWave^m)

MatrixOp/o/c theMean=mean(momentWave)
momentMatOneTrig[n][m]=theMean[0]
endfor
endfor

end

The results are contained in the wave "momentMatOneTrig". The problem is, in case 2, the first row/column will be empty, caused by some nan generated in the mean calculation.
As it now turns out, it's my input data having something weird. I attached it so if you have time you can try and let me know if you spotted anything from the data.
Thanks.
In your double-nested loop, you have this:
momentWave=(ScWave^n)*(SWave^m)
MatrixOp/o/c theMean=mean(momentWave)

I ran your function, like this:
deBuggingIgorCmplx(2)

This created ScWave and SWave.

I then executed this in the command line:
momentWave=(ScWave^0)*(SWave^0)
MatrixOp/o/c theMean=mean(momentWave); Print theMean

This prints "cmplx(NaN,NaN)".

I then did WaveStats on your momentWave, I see that it has 29 NaNs.

I displayed momentWave in a table and used Edit->Find to search for blanks (NaNs). The first NaN occurs at point 1114. The second at 3043.

So the problem starts with
momentWave=(ScWave^0)*(SWave^0)

which stores NaN for point 1114.

Looking at SWave and ScWave for point 1114, I see that both are cmplx(0,0).

I then executed this:
Print cmplx(0,0)^0

This prints "(NaN,NaN)". That explains why point 1114 in momentWave is NaN.

The question is: "Is it correct that cmplx(0,0)^0 is NaN?

According to https://www.quora.com/What-is-the-result-of-a-complex-number-raised-to-…
even in the complex plane, 0^0 is still an indeterminate form

https://en.wikipedia.org/wiki/Zero_to_the_power_of_zero says:
Zero to the power of zero, denoted by 0^0, is a mathematical expression with no necessarily obvious value. Possibilities include 0, 1, or leaving the expression undefined altogether, and there is no consensus as to which approach is best.

This does not address the case of raising a complex zero to the zeroth power, but it lends credence to the idea that Igor's result for "cmplx(0,0)^0" is reasonable.

I wonder what Mathematica returns for cmplx(0,0)^0.

At any rate, that explains why you are getting NaNs in row 0 and column 0 of momentMatOneTrig.

Thanks a lot for the check, I really appreciate it.
That explains it, these data was from a digitizer, so it is safe to say there are no true zeros, I may just manually set them to very small nonzero values for the script to work then.

PS: I finally decided to just drop them, they are around a few hundreds from the 2e5 data points.

I have a related problem with complex waves, in that I can't figure out how to directly assign values to the real and complex parts of a complex wave independently.  The following code (which converts a constructed complex power spectral density to a complex amplitude wave with random phase noise to allow for ifft back to the time domain) compiles because of the for loop, but I actually want the commented out line just before the for loop to work. I've tried a few variations and keep getting a "function not available for this number type" error. I don't see any relevant documentation in the manual or help files:

Function PSDtoTimeDomain(Wi,Rg,trelax,N,b)
//Wave PSD
Variable Wi,Rg,trelax,N,b

Variable viscS,nu,kb,T, tau1,i

//make wave for PSD, dimensionless for now
Make/N=10000 PSD
SetScale/I x, 0, 10, "",PSD     //Wi = 0 overlaps with Wi>0 for f>~10 dimensionless

viscS = 1e-3            // water solvent viscosity Pa.s
nu = 3/5                    // good solvent, nu=3/5  (1/2 for theta solvent
kb = 1.380649e-23       // Boltzmann constant, J/k
T = 300                 // room temperature, K

//tau1 = viscS*N^(3*nu)*b^3/(sqrt(3*Pi)*kb*T)           //relaxation time in s
tau1 = 1
print tau1

PSD=0
for(i=1;i<51;i+=2)  // 25 terms
//Do summation from Hur eq50
PSD[] = PSD[p] + 1/( i^(2/5)*(i^(18/5)+(2*Pi*x*tau1)^2)) + Wi^2/( i^(2/5)*(i^(18/5)+(2*Pi*x*tau1)^2)^2)
endfor                      // Execute body code until continue test is FALSE
//Leave dimensionless for now, for validation with Hur paper dimensionless PSDs
PSD[] = PSD[p]*96/(Pi^2)
//PSD[] = PSD[p]*96*Rg^2*tau1/(Pi^2)                //add shared coefficients

//Convert PSD to time domain stretch signal
Duplicate PSD, PSD_complex          //Create wave to store complex PSD
Redimension/C PSD_complex           // Make it complex
Duplicate PSD_complex, sqrtPSD_complex
sqrtPSD_complex = sqrt(sqrtPSD_complex) //Convert from power to amplitude

//sqrtPSD_complex[] = p2rect(cmplx(real(sqrtPSD_complex[p],enoise(Pi))))
Variable/C temp
for(i=0;i<DimSize(sqrtPSD_complex,0);i+=1)
temp = cmplx(sqrtPSD_complex[i],enoise(Pi))
sqrtPSD_complex[i] = temp
endfor                      // Execute body code until continue test is FALSE
Duplicate PSD, TimeDomain_PSD
ifft/R/DEST=TimeDomain_PSD sqrtPSD_complex

end

Is there a way to do the wave assignment without using the for loop?  The for loop works, but if there is a more concise way to do this, I'd rather do that.

you have a parenthesis problem:

`sqrtPSD_complex = p2rect(cmplx(real(sqrtPSD_complex),enoise(Pi)))`

Oops, you're right. Thank you for finding that, but your line still gives the "function not available for this number type" error.

This simpler construction does too:

`sqrtPSD_complex =     cmplx(real(sqrtPSD_complex),4)`

Here is a version of your original function that compiles:

Function PSDtoTimeDomain(Wi,Rg,trelax,N,b)
//Wave PSD
Variable Wi,Rg,trelax,N,b

Variable viscS,nu,kb,T, tau1,i

//make wave for PSD, dimensionless for now
Make/N=10000 PSD
SetScale/I x, 0, 10, "",PSD     //Wi = 0 overlaps with Wi>0 for f>~10 dimensionless

viscS = 1e-3            // water solvent viscosity Pa.s
nu = 3/5                    // good solvent, nu=3/5  (1/2 for theta solvent
kb = 1.380649e-23       // Boltzmann constant, J/k
T = 300                 // room temperature, K

//tau1 = viscS*N^(3*nu)*b^3/(sqrt(3*Pi)*kb*T)           //relaxation time in s
tau1 = 1
print tau1

PSD=0
for(i=1;i<51;i+=2)  // 25 terms
//Do summation from Hur eq50
PSD[] = PSD[p] + 1/( i^(2/5)*(i^(18/5)+(2*Pi*x*tau1)^2)) + Wi^2/( i^(2/5)*(i^(18/5)+(2*Pi*x*tau1)^2)^2)
endfor                      // Execute body code until continue test is FALSE
//Leave dimensionless for now, for validation with Hur paper dimensionless PSDs
PSD[] = PSD[p]*96/(Pi^2)
//PSD[] = PSD[p]*96*Rg^2*tau1/(Pi^2)                //add shared coefficients

//Convert PSD to time domain stretch signal
Duplicate PSD, PSD_complex          //Create wave to store complex PSD
Redimension/C PSD_complex           // Make it complex
Duplicate/C PSD_complex, sqrtPSD_complex
sqrtPSD_complex = sqrt(sqrtPSD_complex) //Convert from power to amplitude

sqrtPSD_complex[] = p2rect(cmplx(real(sqrtPSD_complex[p]),enoise(Pi)))
//  Variable/C temp
//  for(i=0;i<DimSize(sqrtPSD_complex,0);i+=1)
//      temp = cmplx(sqrtPSD_complex[i],enoise(Pi))
//      sqrtPSD_complex[i] = temp
//  endfor                      // Execute body code until continue test is FALSE
Duplicate PSD, TimeDomain_PSD
ifft/R/DEST=TimeDomain_PSD sqrtPSD_complex

end

In order for the sqrtPSD_complex[] wave assignment statement to compile, Igor's compiler needs to know that this is a complex wave. You can tell it by adding the /C flag to the Duplicate statement that creates that wave.

When dealing with complex waves, make sure to add the /C flag to Duplicate and when you create WAVE reference variables (WAVE/C myComplexWave=root:foo:someComplexWave)

Thank you, that fixed it.  I saw that the output wave was complex after the duplicate command and assumed the compiler would know that too.  I'll start seeing looking for missing /C flags when I run into this sort of error.