Arithmetic operations with complex number
Sandbo
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