ZGainsforth/ReadSPA.py

## ReadSPA.py
# Created 2015, Zack Gainsforth
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import struct
from numpy.fft import fft, fftfreq

def LoadSPAInterferogram(FileName):

    # Open the SPA file.
    with open(FileName, 'rb') as f:

        # Go to offset 294 in the header which tells us the number of sections in the file.
        f.seek(294)
        NumSections = struct.unpack('h', f.read(2))[0]
        print('This spa file has ' + str(NumSections) + ' sections.')

        # We will process each section in the spa file.
        for n in range(NumSections):
            # Go to the section start.  Each section is 16 bytes, and starts after offset 304.
            f.seek(304+16*n)
            SectionType = struct.unpack('h', f.read(2))[0]
            SectionOffset = struct.unpack('i', f.read(3)+b'\x00')[0]
            SectionLength = struct.unpack('i', f.read(3)+b'\x00')[0]

            print('Section #%d, Type: %d, Offset: %d, Length %d' % (n, SectionType, SectionOffset, SectionLength))

            # If this is section type 3, then it contains the result spectrum (interferogram, single beam, etc.) Note where it is in the file.
            if SectionType==3:
                print('Found raw spectrum section.')
                Section3Offset = SectionOffset
                Section3Length = SectionLength
            # If this is section type 102, then it contains a processed spectrum.
            if SectionType==102:
                print('Found processed spectrum section.')
                Section102Offset = SectionOffset
                Section102Length = SectionLength

        # We will default to using section 102 (processed data) if it exists.  Otherwise, we read section 3 (raw data).
        if 'Section102Length' in locals():
            print('Using processed spectrum.')
            SectionOffset = Section102Offset
            SectionLength = Section102Length
        else:
            print('Using raw spectrum.')
            SectionOffset = Section3Offset
            SectionLength = Section3Length

        # Read in the data!
        f.seek(SectionOffset)
        DataText = f.read(SectionLength)

        Data = np.fromstring(DataText, dtype='float32')[:-5]

        return Data

def PlotInterferogram(Interferogram):

    plt.figure()
    plt.plot(Interferogram)
    plt.title('Interferogram')

    # Apply the hamming window to smooth out the FT (we don't want ripples.)
    H = np.hamming(len(Interferogram))
    I = H*Interferogram
    plt.figure()
    plt.plot(I)
    plt.plot(H)
    plt.title('Hamminged Interferogram')

    # Compute the FFT
    F = fft(I)

    # Make the magnitude spectrum.
    Mag = np.abs(F)
    Mag = Mag[:int(len(Mag)/2)]

    LaserFreq = 15798.3 # cm-1, this is the sampling distance for the interferogram.
    #dE = 1/(2*LaserFreq)
    E = fftfreq(len(Interferogram), 1/(2*LaserFreq))
    E=E[:int(len(E)/2)]

    plt.figure()
    plt.plot(E, Mag)
    plt.title('Magnitude spectrum')

    # Make the phase spectrum.
    Phase = np.unwrap(np.angle(F))
    Phase = Phase[:int(len(Phase)/2)]
    plt.figure()
    plt.plot(E, Phase)
    plt.title('Phase spectrum')

    return

def GetMagAndPhaseFromInterferogram(Interferogram, SamplingLaserFreq=15798.3):
    # Returns:
    # Wavenumbers is the energy axis in cm-1
    # Mag is the magnitude spectrum
    # Phase is the phase spectrum.

    # Apply the hamming window to smooth out the FT (we don't want ripples.)
    H = np.hamming(len(Interferogram))
    I = H*Interferogram

    # Compute the FFT
    F = fft(I)

    # Make the magnitude spectrum.
    Mag = np.abs(F)
    Mag = Mag[:int(len(Mag)/2)]

    LaserFreq = SamplingLaserFreq #15798.3 # cm-1, this is the sampling distance for the interferogram.
    #dE = 1/(2*LaserFreq)
    E = fftfreq(len(Interferogram), 1/(2*LaserFreq))
    E=E[:int(len(E)/2)]

    # Make the phase spectrum.
    Phase = np.unwrap(np.angle(F))
    Phase = Phase[:int(len(Phase)/2)]

    return E, Mag, Phase

if __name__ == '__main__':

    # FileName = 'Nuth_FeO_smoke_DriftSpot2_Au_0006.spa'
    FileName = 'background/bkg_au_harmonic2.spa'
    # FileName = 'section5/scan2/spec0040.spa'
    Interferogram = LoadSPAInterferogram(FileName)

    PlotInterferogram(Interferogram)

    plt.show()
	# Created 2015, Zack Gainsforth
	import matplotlib
	import matplotlib.pyplot as plt
	import numpy as np
	import struct
	from numpy.fft import fft, fftfreq

	def LoadSPAInterferogram(FileName):

	# Open the SPA file.
	with open(FileName, 'rb') as f:

	# Go to offset 294 in the header which tells us the number of sections in the file.
	f.seek(294)
	NumSections = struct.unpack('h', f.read(2))[0]
	print('This spa file has ' + str(NumSections) + ' sections.')

	# We will process each section in the spa file.
	for n in range(NumSections):
	# Go to the section start. Each section is 16 bytes, and starts after offset 304.
	f.seek(304+16*n)
	SectionType = struct.unpack('h', f.read(2))[0]
	SectionOffset = struct.unpack('i', f.read(3)+b'\x00')[0]
	SectionLength = struct.unpack('i', f.read(3)+b'\x00')[0]

	print('Section #%d, Type: %d, Offset: %d, Length %d' % (n, SectionType, SectionOffset, SectionLength))

	# If this is section type 3, then it contains the result spectrum (interferogram, single beam, etc.) Note where it is in the file.
	if SectionType==3:
	print('Found raw spectrum section.')
	Section3Offset = SectionOffset
	Section3Length = SectionLength
	# If this is section type 102, then it contains a processed spectrum.
	if SectionType==102:
	print('Found processed spectrum section.')
	Section102Offset = SectionOffset
	Section102Length = SectionLength

	# We will default to using section 102 (processed data) if it exists. Otherwise, we read section 3 (raw data).
	if 'Section102Length' in locals():
	print('Using processed spectrum.')
	SectionOffset = Section102Offset
	SectionLength = Section102Length
	else:
	print('Using raw spectrum.')
	SectionOffset = Section3Offset
	SectionLength = Section3Length

	# Read in the data!
	f.seek(SectionOffset)
	DataText = f.read(SectionLength)

	Data = np.fromstring(DataText, dtype='float32')[:-5]

	return Data

	def PlotInterferogram(Interferogram):

	plt.figure()
	plt.plot(Interferogram)
	plt.title('Interferogram')

	# Apply the hamming window to smooth out the FT (we don't want ripples.)
	H = np.hamming(len(Interferogram))
	I = H*Interferogram
	plt.figure()
	plt.plot(I)
	plt.plot(H)
	plt.title('Hamminged Interferogram')

	# Compute the FFT
	F = fft(I)

	# Make the magnitude spectrum.
	Mag = np.abs(F)
	Mag = Mag[:int(len(Mag)/2)]

	LaserFreq = 15798.3 # cm-1, this is the sampling distance for the interferogram.
	#dE = 1/(2*LaserFreq)
	E = fftfreq(len(Interferogram), 1/(2*LaserFreq))
	E=E[:int(len(E)/2)]

	plt.figure()
	plt.plot(E, Mag)
	plt.title('Magnitude spectrum')

	# Make the phase spectrum.
	Phase = np.unwrap(np.angle(F))
	Phase = Phase[:int(len(Phase)/2)]
	plt.figure()
	plt.plot(E, Phase)
	plt.title('Phase spectrum')

	return

	def GetMagAndPhaseFromInterferogram(Interferogram, SamplingLaserFreq=15798.3):
	# Returns:
	# Wavenumbers is the energy axis in cm-1
	# Mag is the magnitude spectrum
	# Phase is the phase spectrum.

	# Apply the hamming window to smooth out the FT (we don't want ripples.)
	H = np.hamming(len(Interferogram))
	I = H*Interferogram

	# Compute the FFT
	F = fft(I)

	# Make the magnitude spectrum.
	Mag = np.abs(F)
	Mag = Mag[:int(len(Mag)/2)]

	LaserFreq = SamplingLaserFreq #15798.3 # cm-1, this is the sampling distance for the interferogram.
	#dE = 1/(2*LaserFreq)
	E = fftfreq(len(Interferogram), 1/(2*LaserFreq))
	E=E[:int(len(E)/2)]

	# Make the phase spectrum.
	Phase = np.unwrap(np.angle(F))
	Phase = Phase[:int(len(Phase)/2)]

	return E, Mag, Phase

	if __name__ == '__main__':

	# FileName = 'Nuth_FeO_smoke_DriftSpot2_Au_0006.spa'
	FileName = 'background/bkg_au_harmonic2.spa'
	# FileName = 'section5/scan2/spec0040.spa'
	Interferogram = LoadSPAInterferogram(FileName)

	PlotInterferogram(Interferogram)

	plt.show()