FFT Programming (Magnitude and Phase), Please Help

Share icon

Fri, 04/19/2013 - 12:25 am

Hi everyone,

I've been using Igor Pro for about five years now, however I have only just ventured into the art of programming. Currently I am attempting to modify a prexisting experiment that performs FFT on interferometric fringe patterns. The original function listed below works well but only outputs the Magnitude of the Fourier Transform:

Function RealTimeFFT2(w,windowing,resolution,limits)    //designed to do a fast FFT of wave named "w"
//w is considered to be an optical spectrum from the ocean optics CCD.  The routine takes the wavelength x-axis from
//"Wavelength" and converts it to the proper frequency scaling before FFT.
    string w
    variable windowing  //This is the windowing: 1= None; 2=Hanning
    variable resolution //Resolution enhancement: 1, 2, 3, 4, 5, 6 are: none;2;4;8;16;32 respectively.  Actually just determines how many data points are in the FFT
    variable limits
//Modified from Wavemetrics Library
//Given an input wave ( doesn't have to be a power of two), creates a new wave
//with the suffix "_Mag" that contains the normalized frequency response
//Several levels of resolution enhancement (really just sin x/x interpolation) are provided.
//Options include windowing (Hann) vs no windowing.
    string destw=w+"_Mag"
    Variable n      //= numpnts($w)
    variable start_point, end_point
    if (limits) //want to do FFT only on portion of spectrum selected between the cursors in the spectrum window
        NVAR GCursorAPoint_spectrum=GCursorAPoint_spectrum
        NVAR GCursorBPoint_spectrum=GCursorBPoint_spectrum
        GCursorAPoint_spectrum=pcsr(A,"Graph4")
        GCursorBPoint_spectrum=pcsr(B,"Graph4")
        start_point=GCursorAPoint_spectrum
        end_point=GCursorBPoint_spectrum
    else    //do FFT on entire spectrum
        start_point=0
        end_point=numpnts($w)-1
    endif
    if( (resolution<1) %| (resolution>6) )
        Abort "resolution out of range"
    endif
    resolution=2^resolution     //convert this resolution value to the appropriate binary number (if file length is 2^n, it speeds up the FFT algorithm)
    //resolution -= 1
    Duplicate/O/R=[start_point,end_point] Wavelength, xData
    xData=1/xData       //X axis will be 1/nm so FT will be in nm
    Duplicate/O/R=[start_point,end_point] $w, yData
    Sort xData, xData, yData    // sort by increasing X
    Duplicate/O yData, tempFFT  // make tempFFT file to contain properly scaled waveform
    SetScale/I x xData[0], xData[(numpnts(xData)-1)], tempFFT   // set the range of X values
    tempFFT = interp(x, xData, yData)
    if(windowing==2)
        Hanning tempFFT; tempFFT *= 2           //assumes continuous rather than pulsed data
    endif
    n= numpnts(yData)
    Redimension/N=(CeilPwr2(n)*2^resolution) tempFFT        //pad with zeros to power of 2
    //print "resolution is " +num2str(resolution)
    //print "Ceilpwr2(n) is " +num2str(CeilPwr2(n))
    //print "number of new data points is " +num2str(CeilPwr2(n)*2^resolution)
    FFT tempFFT
    Wave/C tempFFTC= tempFFT        //this resets the complex part of FTspectrum, which was generated in the fft
    tempFFTC=r2polar(tempFFT)
    Redimension/R tempFFTC
    //NOTE: depending on your application you may want to un-comment the next line
//  $destw[0] /= 2                              //dc is special
    tempFFT /= n/2
    SetScale y,0,0,"Cts",tempFFT
    Redimension/R tempFFT
    Duplicate/O tempFFT $destw
end


I would like to also calculate the phase angle from the FFT. I am basing my modification on the example experiement "MagPhase Demo". As I have previously mentioned this is my first attempt at igor pro programming and while I am slowly learning how it all works, I keep getting an error code near the end of the function, "Function not available for this numbertype", relating to the line ----- tempFFT_Phase = imag(r2polar(tempFFT2)). I believe that this indicates that the wave has not been identified as complex, however I thought I had addressed this previously with ----- Wave/C tempFFT_Phase= tempFFT2.

Any help would be greatly appreciated and as I am still a beginner, example code would be fantastic.

Function RealTimeFFTMagPhase(w,windowing,resolution,limits) //designed to do a fast FFT of wave named "w"
//w is considered to be an optical spectrum from the ocean optics CCD.  The routine takes the wavelength x-axis from
//"Wavelength" and converts it to the proper frequency scaling before FFT.
    string w
    variable windowing  //This is the windowing: 1= None; 2=Hanning
    variable resolution //Resolution enhancement: 1, 2, 3, 4, 5, 6 are: none;2;4;8;16;32 respectively.  Actually just determines how many data points are in the FFT
    variable limits
//Modified from Wavemetrics Library
//Given an input wave ( doesn't have to be a power of two), creates a new wave
//with the suffix "_Mag" that contains the normalized frequency response
//Several levels of resolution enhancement (really just sin x/x interpolation) are provided.
//Options include windowing (Hann) vs no windowing.
    string destw_Mag, destw_Phase
    destw_Mag=w+"_Mag"      //Destination file for Magnitude
    destw_Phase=w+"_Phase"  //Destingation file for Phase
   
    Variable n      //= numpnts($w)
    variable start_point, end_point
    if (limits) //want to do FFT only on portion of spectrum selected between the cursors in the spectrum window
        NVAR GCursorAPoint_spectrum=GCursorAPoint_spectrum
        NVAR GCursorBPoint_spectrum=GCursorBPoint_spectrum
        GCursorAPoint_spectrum=pcsr(A,"Graph4")
        GCursorBPoint_spectrum=pcsr(B,"Graph4")
        start_point=GCursorAPoint_spectrum
        end_point=GCursorBPoint_spectrum
    else    //do FFT on entire spectrum
        start_point=0
        end_point=numpnts($w)-1
    endif
    if( (resolution<1) %| (resolution>6) )
        Abort "resolution out of range"
    endif
    resolution=2^resolution     //convert this resolution value to the appropriate binary number (if file length is 2^n, it speeds up the FFT algorithm)
    //resolution -= 1
    Duplicate/O/R=[start_point,end_point] Wavelength, xData
    xData=1/xData       //X axis will be 1/nm so FT will be in nm
    Duplicate/O/R=[start_point,end_point] $w, yData
    Sort xData, xData, yData    // sort by increasing X
    Duplicate/O yData, tempFFT  // make tempFFT file to contain properly scaled waveform
    SetScale/I x xData[0], xData[(numpnts(xData)-1)], tempFFT   // set the range of X values
    tempFFT = interp(x, xData, yData)
    if(windowing==2)
        Hanning tempFFT; tempFFT *= 2           //assumes continuous rather than pulsed data
    endif
    n= numpnts(yData)
    Redimension/N=(CeilPwr2(n)*2^resolution) tempFFT        //pad with zeros to power of 2
    //print "resolution is " +num2str(resolution)
    //print "Ceilpwr2(n) is " +num2str(CeilPwr2(n))
    //print "number of new data points is " +num2str(CeilPwr2(n)*2^resolution)
    FFT tempFFT
    Duplicate/O/C tempFFT, tempFFT2     //Duplicated Fourier Transformed wave for Mag and Phase analysis
    Wave/C tempFFT_Mag= tempFFT     //this resets the complex part of FTspectrum, which was generated in the fft
    Wave/C tempFFT_Phase= tempFFT2
    tempFFT_Mag=r2polar(tempFFT)
    Redimension/R tempFFT_Mag
    //NOTE: depending on your application you may want to un-comment the next line
//  $destw_Mag[0] /= 2                              //dc is special
    tempFFT /= n/2
    SetScale y,0,0,"Cts",tempFFT
    Redimension/R tempFFT
    Duplicate/O tempFFT $destw_Mag
    tempFFT_Phase = imag(r2polar(tempFFT2))  //**************ERROR**********  "Function not available for this numbertype"
    tempFFT_Phase *= 180/PI                     // Convert phase to degrees
    Redimension/R tempFFT_Phase
    tempFFT2 /= n/2
    SetScale y,0,0,"Cts",tempFFT2
    Redimension/R tempFFT2
    Duplicate/O tempFFT2 $destw_Phase                  
end


Thanks

Andy
imag() returns a real number, not a complex one. So Igor throws an error because you declared tempFFT_Phase to be complex.

April 19, 2013 at 01:41 am - Permalink

Thanks 741, changing

Wave/C tempFFT_Phase= tempFFT2
to
Wave tempFFT_Phase= tempFFT2

cleared up the error associated with the imag() tag, however now I receive an error on the very next line:

tempFFT_Phase *= 180/PI wave assignment error: "Complex wave used in real expression."

I thought that the imag() tag converted the imaginary component of the complex wave into a real number, as shown in the example experiment MagPhase Demo.

Thanks again for your fast reply

April 19, 2013 at 07:38 am - Permalink

baitzer wrote:

I thought that the imag() tag converted the imaginary component of the complex wave into a real number, as shown in the example experiment MagPhase Demo.


There is no conversion here. imag() is a function that returns a real value. If you assign it to a complex wave the compiler will find it objectionable.

The solution to many of these problems is to use MatrixOP. MatrixOP does not care if your input is real or complex. It attempts to give you the correct output as implied by the return values of all the tokens that appear on the right hand side of the MatrixOP expression.

Also note that you were apparently working from a very old example. FFT has a /DEST flag and /OUT flag that would make your code a bit more compact.

A.G.
WaveMetrics, Inc.

April 19, 2013 at 10:08 am - Permalink

Thanks A.G.

Modified my function to include MatrixOP and now it works great.

April 20, 2013 at 04:31 am - Permalink

Hi all,

Thanks again for your earlier help. I just have a few additional questions relating to the phase angle. I modified my original code according to A.G. and 741's suggestions which solved the associated error's I was receiving.

New Code:

Function RealTimeFFTMagPhase(w,windowing,resolution,limits) //designed to do a fast FFT of wave named "w"
//w is considered to be an optical spectrum from the ocean optics CCD.  The routine takes the wavelength x-axis from
//"Wavelength" and converts it to the proper frequency scaling before FFT.
    string w
    variable windowing  //This is the windowing: 1= None; 2=Hanning
    variable resolution //Resolution enhancement: 1, 2, 3, 4, 5, 6 are: none;2;4;8;16;32 respectively.  Actually just determines how many data points are in the FFT
    variable limits
//Modified from Wavemetrics Library
//Given an input wave ( doesn't have to be a power of two), creates a new wave
//with the suffix "_Mag" that contains the normalized frequency response
//Several levels of resolution enhancement (really just sin x/x interpolation) are provided.
//Options include windowing (Hann) vs no windowing.
    string destw_Mag, destw_Phase
    destw_Mag=w+"_Mag"
    destw_Phase=w+"_Phase"
   
    Variable n      //= numpnts($w)
    variable start_point, end_point
    if (limits) //want to do FFT only on portion of spectrum selected between the cursors in the spectrum window
        NVAR GCursorAPoint_spectrum=GCursorAPoint_spectrum
        NVAR GCursorBPoint_spectrum=GCursorBPoint_spectrum
        GCursorAPoint_spectrum=pcsr(A,"Graph4")
        GCursorBPoint_spectrum=pcsr(B,"Graph4")
        start_point=GCursorAPoint_spectrum
        end_point=GCursorBPoint_spectrum
    else    //do FFT on entire spectrum
        start_point=0
        end_point=numpnts($w)-1
    endif
    if( (resolution<1) %| (resolution>6) )
        Abort "resolution out of range"
    endif
    resolution=2^resolution     //convert this resolution value to the appropriate binary number (if file length is 2^n, it speeds up the FFT algorithm)
    //resolution -= 1
    Duplicate/O/R=[start_point,end_point] Wavelength, xData
    xData=1/xData       //X axis will be 1/nm so FT will be in nm
    Duplicate/O/R=[start_point,end_point] $w, yData
    Sort xData, xData, yData    // sort by increasing X
    Duplicate/O yData, tempFFT  // make tempFFT file to contain properly scaled waveform
    SetScale/I x xData[0], xData[(numpnts(xData)-1)], tempFFT   // set the range of X values
    tempFFT = interp(x, xData, yData)
    if(windowing==2)
        Hanning tempFFT; tempFFT *= 2           //assumes continuous rather than pulsed data
    endif
    n= numpnts(yData)
    Redimension/N=(CeilPwr2(n)*2^resolution) tempFFT        //pad with zeros to power of 2
    //print "resolution is " +num2str(resolution)
    //print "Ceilpwr2(n) is " +num2str(CeilPwr2(n))
    //print "number of new data points is " +num2str(CeilPwr2(n)*2^resolution)
    FFT tempFFT
    Duplicate/O/C tempFFT, tempFFT2
    Wave/C tempFFT_Mag= tempFFT     //this resets the complex part of FTspectrum, which was generated in the fft
    Wave tempFFT_Phase= tempFFT2
    tempFFT_Mag=(r2polar(tempFFT))
    Redimension/R tempFFT_Mag
    //NOTE: depending on your application you may want to un-comment the next line
//  $destw_Mag[0] /= 2                              //dc is special
    tempFFT /= n/2
    SetScale y,0,0,"Cts",tempFFT
    Redimension/R tempFFT
    Duplicate/O tempFFT $destw_Mag
    Variable degree=(180/pi)
    MatrixOP/O tempFFT_Phase=degree*imag(tempFFT2)  // Extracts the imaginary component of the complex wave (Phase) and then Convert phase to degrees via (180/PI)
    Print "test", atan(tempFFT2)
    Redimension/R tempFFT_Phase
    tempFFT2 /= n/2    
    SetScale y,0,0,"Deg",tempFFT2
    Redimension/R tempFFT2     
    Duplicate/O tempFFT2 $destw_Phase
    Print "Phase", tempFFT2
end


I have now managed to plot both the Magnitude and Phase on a single graph as shown in the attached image (Magnitude = red, phase = blue), however the phase is massive >20 kDegrees. I thought that the phase extracted from the FFT was between -Pi and Pi or -180 to 180 degrees if converted.

http://www.igorexchange.com/files/images/FFT%20Magnitude%20&%20Phase%20(Zoom).JPG

Does anyone have a suggestion as to why my code is creating such a large phase?

Also the project I am modifying the code for requires the phase angle at the magnitude. Now that I have extracted the phase information from the FFT, how would I determine the phase angle. I can easily find the magnitude via peak fitting the red trace but I am unsure how to determine the related phase angle.

Cheers

April 26, 2013 at 10:26 pm - Permalink

You can convert between rectangular (real, imag) and polar (mag, phase) representations of complex numbers using the functions r2polar() and p2rect(). With phase variables sometimes the Unwrap operation comes in handy.

April 27, 2013 at 07:22 am - Permalink