Error in FFT example from manual. "Complex wave used in real expression"

Share icon

Wed, 10/07/2020 - 02:29 pm

The Manual gives the following example in II-5: Waves

    Make/O/N=1024 wave0
    SetScale x 0,1,"s", wave0
    wave0 = sin(2*PI*x) + sin(6*PI*x)/3 + sin(10*PI*X)/5 + sin(14*PI*x)/7
    Display wave0 as "Time Domain"
   
    // Do FFT
    Duplicate/O wave0, cwave0
    FFT cwave0
    cwave0 /= 512; cwave0[0] /= 2
    Display cwave0 as "Frequency Domain";SetAxis bottom, 0, 25
   
    // Calcualate magnitude and phase
    Make/O/N=513 mag0, phase0, power0
    CopyScales cwave0, mag0, phase0, power0
    mag0 = real(r2polar(cwave0))
    phase0 = imag(r2polar(cwave0))
    phase0 *= 180/PI
    Display mag0 as "Magnitude and Phase"; AppendToGraph/R phase0
    SetAxis bottom, 0, 25
    Label left, "Magnitude";Label right, "Phase"
   
    // Calculate power spectrum
    power0 = magsqr(cwave0)
    Display power0 as "Power Sepctrum";SetAxis bottom, 0, 25

But the line "cwave /= 512; cwave0[0] /= 2" gives an error. "Complex wave used in a real expression." What would be the proper way to renormalize the amplitude of cwave0? Sort of discouraging that the manual wouldn't do it correctly.

The code was designed to be run from a Macro or command line many years ago.

Putting the code into a function wasn't tried after Igor added more stringent checking of wave reference types.

That checking is placated by using a WAVE/C wave reference on the FFT result:

Function fComplexFromManual()

    // Make a time domain waveform
    Make/O/N=1024 wave0
    SetScale x 0, 1, "s", wave0   // goes from 0 to 1 second
    wave0=sin(2*PI*x)+sin(6*PI*x)/3+sin(10*PI*x)/5+sin(14*PI*x)/7
    Display wave0 as "Time Domain"

    // Do the FFT
    Duplicate/O wave0, cwave0               // get copy to do FFT on
    FFT cwave0                              // cwave0 is now complex...
    WAVE/C cw0 = cwave0                     // ... so we need a complex wave ref
    cw0 /= 512;cw0[0] /= 2                  // normalize amplitude
    Display cw0 as "Frequency Domain";SetAxis bottom, 0, 25

    // Calculate magnitude and phase
    Make/O/N=513 mag0, phase0, power0           // these are real waves
    CopyScales cw0, mag0, phase0, power0
    mag0 = real(r2polar(cw0))
    phase0 = imag(r2polar(cw0))
    phase0 *= 180/PI                                    // convert to degrees
    Display mag0 as "Magnitude and Phase";AppendToGraph/R phase0
    SetAxis bottom, 0, 25
    Label left, "Magnitude";Label right, "Phase"

    // Calculate power spectrum
    power0 = magsqr(cw0)
    Display power0 as "Power Spectrum";SetAxis bottom, 0, 25
End

October 7, 2020 at 03:28 pm - Permalink

This version runs in both a function and in a macro (macros don't permit or need WAVE statements):

// Make a time domain waveform
Make/O/N=1024 wave0
SetScale x 0, 1, "s", wave0         // goes from 0 to 1 second
wave0=sin(2*PI*x)+sin(6*PI*x)/3+sin(10*PI*x)/5+sin(14*PI*x)/7
Display wave0 as "Time Domain"

// Do the FFT
FFT/DEST=cwave0 wave0              // cwave0 is created complex, 513 points
cwave0 /= 512;cwave0[0] /= 2       // normalize amplitude
Display cwave0 as "Frequency Domain";SetAxis bottom, 0, 25

// Calculate magnitude and phase
Make/O/N=513 mag0, phase0, power0  // these are real waves
CopyScales cwave0, mag0, phase0, power0
mag0 = real(r2polar(cwave0))
phase0 = imag(r2polar(cwave0))
phase0 *= 180/PI                   // convert to degrees
Display mag0 as "Magnitude and Phase";AppendToGraph/R phase0
SetAxis bottom, 0, 25
Label left, "Magnitude";Label right, "Phase"

// Calculate power spectrum
power0 = magsqr(cwave0)
Display power0 as "Power Spectrum";SetAxis bottom, 0, 25

October 7, 2020 at 03:50 pm - Permalink