# Arithmetic operations with complex number

Sandbo

Tue, 05/29/2018 - 03:46 pm

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)

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

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.

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)?

May 29, 2018 at 04:05 pm - Permalink

`MatrixOP theMean = mean(complexWave)`

John Weeks

WaveMetrics, Inc.

support@wavemetrics.com

May 29, 2018 at 05:07 pm - Permalink

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?

June 1, 2018 at 10:31 am - Permalink

John Weeks

WaveMetrics, Inc.

support@wavemetrics.com

June 1, 2018 at 02:09 pm - Permalink

I tried Igor 6 and 7, they were also affected.

I am still on Igor 6 mainly, but moving to Igor 7 lately.

June 1, 2018 at 06:40 pm - Permalink

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:

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

June 1, 2018 at 08:09 pm - Permalink

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

June 2, 2018 at 04:22 am - Permalink

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):

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.

June 2, 2018 at 06:22 am - Permalink

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:

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.

June 2, 2018 at 03:44 pm - Permalink

MatrixOp/o/c theMean=mean(momentWave)

I ran your function, like this:

This created ScWave and SWave.

I then executed this in the command line:

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

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:

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.

June 2, 2018 at 07:56 pm - Permalink

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.

June 4, 2018 at 09:31 am - Permalink

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:

//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

//Add phase noise

//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.

April 27, 2023 at 05:31 am - Permalink

you have a parenthesis problem:

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

April 27, 2023 at 08:19 am - Permalink

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)`

April 27, 2023 at 07:44 pm - Permalink

Here is a version of your original function that compiles:

//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

//Add phase noise

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)

April 28, 2023 at 05:38 am - Permalink

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.

April 28, 2023 at 05:42 am - Permalink

In reply to Thank you, that fixed it. I… by jmcclimo2003

Yeah, that's a common mistake to make. Unless told to do otherwise, the compiler will default to compiling a wave as a real numeric wave. So without the /C flag, the compiler assumes the source and destination waves of the Duplicate statement are both real numeric.

You'll run into the same kind of problem if you deal with text waves or the less frequently used wave/df reference waves.

April 28, 2023 at 06:04 am - Permalink