Arithmetic operations with complex number

Share icon

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

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

May 29, 2018 at 04:05 pm - Permalink

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?

June 1, 2018 at 10:31 am - Permalink

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.

June 1, 2018 at 06:40 pm - Permalink

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

June 1, 2018 at 08:09 pm - Permalink

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
MultiThreadingControl setMode=0

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

AG

June 2, 2018 at 04:22 am - Permalink

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.

June 2, 2018 at 06:22 am - Permalink

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.
testIQdata.ibw

June 2, 2018 at 03:44 pm - Permalink

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.

June 2, 2018 at 07:56 pm - Permalink

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.

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:

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

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

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