FFT Programming (Magnitude and Phase), Please Help

Share icon

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