In [859]:
%reset
%load_ext autoreload
%autoreload 2
%matplotlib qt
In [860]:
%qtconsole
In [ ]:
%matplotlib?
In [861]:
import matplotlib.pyplot as plt
import numpy as np
import os
from dataIO.span import span
from pyGeneralRoutines.fn_add_subfix import fn_add_subfix
from skimage import data as skidata
from skimage import img_as_float
#from plot_grains import *
#from radialProfiles import *
#from peak_finder import *
from astropy.io import fits
from pyProfile.psd import plot_sig_psd,psd2prof
from IPython.display import display
In [ ]:
"""#Interesting data in header
BITPIX = -64 / number of bits per data pixel
NAXIS = 2 / number of data axes
NAXIS1 = 2048 / length of data axis 1
NAXIS2 = 2048 / length of data axis 2
OBJECT = 'Talysurf CCI' / object_name
OPERATOR= 'TalysurfCCI L5xZ1B1S1F0Hps' / operator_name
POINTSIZ= 32 / pointsize
ZMIN = 112364268 / zmin
ZMAX = 112757161 / zmax
XRES = 2048 / xres
YRES = 2048 / yres
NOFPOINT= 4194304 / nofpoints
DX = 1.6171875 / dx
DY = 1.6171875 / dy
DZ = 9.99999974737875E-06 / dz
DX_UNIT = ' m ' / dx_unit
DY_UNIT = ' m ' / dy_unit
DZ_UNIT = ' m ' / dz_unit
XLENGTHU= 'mm ' / xlength_unit
YLENGTHU= 'mm ' / ylength_unit
ZLENGTHU= ' m ' / zlength_unit
XUNITRAT= 1000. / xunit_ratio
YUNITRAT= 1000. / yunit_ratio
ZUNITRAT= 1. / zunit_ratio
IMPRINT = 0 / imprint
INVERSIO= 0 / inversion
LEVELING= 0 / leveling
SECONDS = 27 / seconds
MINUTES = 46 / minutes
HOURS = 15 / hours
DAY = 20 / day
MONTH = 1 / month
YEAR = 2016 / year
DAYOF = 4 / dayof
XOFFSET = 84.8548202514648 / XOffset
YOFFSET = 70.7438201904297 / YOffset
ZOFFSET = 67264.828125 / ZOffset
"""
Meaning of frequency and PSD vectors¶
In [863]:
#make a test signal
from pyProfile.psd import psd
from pyProfile.profile import line,make_signal
N=500
noise=0
amp=10.
L=30.
nwaves=5.8 #32
ystartend=(3,20)
#ystartend=(0,0)
outfolder='PSDtest'
In [5]:
x,y=make_signal(amp=10.,L=30.,N=500,nwaves=5.8,ystartend=(0,0),noise=0)
f,p=psd(x,y)
print np.std(y)**2
print np.trapz(p/len(y),f)*1e6
print '#--------------------'
x,y=make_signal(amp=10.,L=30.,N=500,nwaves=5.8,ystartend=(3,20),noise=0)
f2,p2=psd(x,y)
print np.std(y)**2
print np.trapz(p/len(y),f)*1e6
plt.plot(f,p/len(x),label='p1')
plt.plot(f2,p2/len(x),label='p2')
plt.loglog()
plt.grid(1)
plt.legend(loc=0)
Out[5]:
In [6]:
print np.std(y)**2
In [7]:
a=np.trapz(np.abs(np.fft.rfft(y))**2/len(x)**2)*2
print a
In [8]:
b=np.trapz(np.abs(np.fft.rfft(y))**2/len(x)**2,np.fft.fftfreq(len(x),span(x,1)/len(x))[:251])*2
print b, b/(f[2]-f[1])
In [9]:
print a/b,b/a
In [10]:
print f[2]-f[1],x[2]-x[1]
Make a sin wave with a few points and calculate PSD. This is the easiest case, a single component (n=5) and no offset or tilt
In [867]:
%matplotlib qt
plt.ion()
x,y=make_signal(amp=10.,L=30.,N=20,nwaves=5.,ystartend=(0,0),noise=0)
f,p=psd(x,y)
fig=plt.figure(11)
plt.clf()
axes=plot_sig_psd(x,y,label='5.0')
plt.cla()
plt.plot(f,p,'o')
plt.grid()
plt.show()
#ax2.set_xscale('linear')
#ax2.set_yscale('linear
display(plt.gcf())
In [173]:
print "rms**2=",np.std(y)**2, '(rms=%f5.3)'%np.std(y)
print "integral of PSD is ",np.trapz(p/len(y),f)
In [49]:
#print all Fourier components
p=np.fft.fft(y,len(x))
print p
In [50]:
#With even number of points, these are symmetries:
print p[0],p[10]
print p[1],p[-1]
In [59]:
#xorresponding frequencies are:
ff=np.fft.fftfreq(len(x),span(x,1)/len(x))
print ff[0],ff[10]
print ff[1],ff[-1]
In [52]:
#With odd number of points it is
p=np.fft.fft(y[:-1],len(x[:-1]))
print p[0],p[10]
print p[1],p[-1]
In [53]:
#where symmetry is the same, but the central point is missing
print p
Fig. 11: Effects of leveling¶
In [291]:
#generic messy signal
x,y=make_signal(amp=10.,L=30.,N=200,nwaves=5.4,ystartend=(29,42),noise=5)
f,p=psd(x,y,scale=1)
print "rms**2=",np.std(y)**2, '(rms=%f5.3)'%np.std(y)
print "integral of PSD is ",p.sum()-p[0]
print "psd[0]=",p[0]
In [293]:
#usually all three give similar results, with the first two being exactly
# equivalent. The third one gives a different number, very different
# if component 0 is large.
print np.trapz(p[1:]*span(x,1),f[1:]) #this gives correct result, but formula miss a /L
#factor respect Press et al for normalization.
print np.trapz(p[1:]) #indeed the formula as in ref. works if
#integral with default dx=1.
print p.sum()-p[0] #this is because formula is defined on
#discrete sums, not on integrals.
The simplest sine wave has wavelength integer fraction of the scanned length L. Extremes are 0.
In [24]:
#simplest sine wave
x,y=make_signal(amp=10.,L=30.,N=200,nwaves=5.,ystartend=(0,0),noise=0)
f,p=psd(x,y,scale=1.)
fig=plt.figure(11)
plt.clf()
axes=plot_sig_psd(x,y,scale=1.)
plt.cla()
plt.plot(f,p,'o',label='5.0')
plt.semilogy()
plt.grid()
plt.show()
#ax2.set_xscale('linear')
#ax2.set_yscale('linear
plt.plot(f[0],p[0],'Db',ms=10)
plt.legend(loc=0)
plt.annotate('note the very small term for 0 frequency.',(f[0],p[0]),(1.,p[0]))
display(plt.gcf())
print "rms**2=",np.std(y)**2, '(rms=%f5.3)'%np.std(y)
print "integral of PSD is ",p.sum()-p[0]
print "psd[0]=",p[0]
In [25]:
#wave of non multiple frequency
#%matplotlib inline
x,y=make_signal(amp=10.,L=30.,N=200,nwaves=5.7,ystartend=(0,0),noise=0)
f,p=psd(x,y,scale=1)
plot_sig_psd(x,y,scale=1.,label='5.7',c='g')
#plt.cla()
#ax1.plot(x,y)
#ax2.plot(f,p,'o')
ax2.set_xscale('linear')
plt.plot(f[0:1],p[0:1],'Dg',ms=10)
plt.annotate('much higher for 5.7',(f[0],p[0]),(1.,p[0]))
plt.grid(1)
plt.show()
ax1.plot(x,line(x,y),'--g')
from IPython.display import display
display(plt.gcf())
print "rms**2=",np.std(y)**2, '(rms=%f5.3)'%np.std(y)
print "integral of PSD is ",p.sum()-p[0]
print "psd[0]=",p[0]
In [26]:
#subtract line from 5.7
y=y-line(x,y)
f,p=psd(x,y,scale=1)
plot_sig_psd(x,y,scale=1.,label='line sub.',c='r',linestyle='--')
ax2.set_xscale('linear')
plt.plot(f[0],p[0],'dr',ms=10)
plt.legend(loc=0)
display(plt.gcf())
print "rms**2=",np.std(y)**2, '(rms=%f5.3)'%np.std(y)
print "integral of PSD is ",p[1:].sum()
print "psd[0]=",p[0]
In [311]:
plt.savefig(os.path.join(outfolder,'leveling.png'))
Fig.12, 13: Effects of windowing¶
Fig. 12: Effects on the simplest sinusoid¶
The simplest sine wave has wavelength integer fraction of the scanned length L. Extremes are 0.
In [28]:
#simplest sine wave
x,y=make_signal(amp=10.,L=30.,N=200,nwaves=5.,ystartend=(0,0),noise=0)
f,p=psd(x,y,scale=1.)
fig=plt.figure(12)
plt.clf()
axes=plot_sig_psd(x,y,scale=1.)
plt.cla()
plt.plot(f,p,'o',label='5.0')
plt.semilogy()
plt.grid()
plt.show()
#ax2.set_xscale('linear')
#ax2.set_yscale('linear
plt.plot(f[0],p[0],'Db',ms=10)
plt.legend(loc=0)
plt.annotate('note the very small term for 0 frequency.',(f[0],p[0]),(1.,p[0]))
display(plt.gcf())
In [30]:
#%matplotlib inline
f,p=psd(x,y*np.hanning(len(y)),scale=1)
plot_sig_psd(x,y*np.hanning(len(y)),scale=1.,label='hanning',c='g')
#plt.cla()
#ax1.plot(x,y)
#ax2.plot(f,p,'o')
ax2.set_xscale('linear')
plt.plot(f[0:1],p[0:1],'Dg',ms=10)
plt.grid(1)
plt.show()
ax1.plot(x,line(x,y),'--g')
plt.legend(loc=0)
display(plt.gcf())
Effect on non fractional sinusoid¶
In [333]:
fig=plt.figure(13)
plt.clf()
x,y=make_signal(amp=10.,L=30.,N=200,nwaves=5.7,ystartend=(0,0),noise=0)
f,p=psd(x,y,scale=1)
plot_sig_psd(x,y,scale=1.,label='line sub.',c='b',linestyle='-')
ax2.set_xscale('linear')
plt.plot(f[0],p[0],'db',ms=10)
plt.legend(loc=0)
display(plt.gcf())
print "rms**2=",np.std(y)**2, '(rms=%f5.3)'%np.std(y)
print "integral of PSD is ",p[1:].sum()
print "psd[0]=",p[0]
In [334]:
#hanning window
f,p=psd(x,y*np.hanning(len(y)),scale=1)
plot_sig_psd(x,y*np.hanning(len(y)),scale=1.,label='hanning',c='g')
#plt.cla()
#ax1.plot(x,y)
#ax2.plot(f,p,'o')
ax2.set_xscale('linear')
plt.plot(f[0:1],p[0:1],'Dg',ms=10)
plt.grid(1)
plt.show()
ax1.plot(x,line(x,y),'--g')
plt.legend(loc=0)
display(plt.gcf())
In [335]:
#line subtraction
y=y-line(x,y)
f,p=psd(x,y,scale=1)
plot_sig_psd(x,y,scale=1.,label='line sub.',c='r',linestyle='-')
ax2.set_xscale('linear')
plt.plot(f[0],p[0],'dr',ms=10)
plt.legend(loc=0)
display(plt.gcf())
print "rms**2=",np.std(y)**2, '(rms=%f5.3)'%np.std(y)
print "integral of PSD is ",p[1:].sum()
print "psd[0]=",p[0]
In [336]:
#add window
plot_sig_psd(x,y*np.hanning(len(y)),scale=1.,label='line sub., hanning',c='m')
#plt.cla()
#ax1.plot(x,y)
#ax2.plot(f,p,'o')
ax2.set_xscale('linear')
plt.plot(f[0:1],p[0:1],'Dm',ms=10)
plt.grid(1)
plt.show()
plt.legend(loc=0)
display(plt.gcf())
In [ ]:
In [337]:
fig=plt.figure(13)
plt.clf()
x,y=make_signal(amp=10.,L=30.,N=200,nwaves=5.7,ystartend=(20,44),noise=0)
f,p=psd(x,y,scale=1)
plot_sig_psd(x,y,scale=1.,label='line sub.',c='b',linestyle='-')
ax2.set_xscale('linear')
plt.plot(f[0],p[0],'db',ms=10)
plt.legend(loc=0)
display(plt.gcf())
print "rms**2=",np.std(y)**2, '(rms=%f5.3)'%np.std(y)
print "integral of PSD is ",p[1:].sum()
print "psd[0]=",p[0]
In [338]:
#line subtraction
y=y-line(x,y)
f,p=psd(x,y,scale=1)
plot_sig_psd(x,y,scale=1.,label='line sub.',c='r',linestyle='-')
ax2.set_xscale('linear')
plt.plot(f[0],p[0],'dr',ms=10)
plt.legend(loc=0)
display(plt.gcf())
print "rms**2=",np.std(y)**2, '(rms=%f5.3)'%np.std(y)
print "integral of PSD is ",p[1:].sum()
print "psd[0]=",p[0]
In [339]:
#add window
plot_sig_psd(x,y*np.hanning(len(y)),scale=1.,label='line sub., hanning',c='m')
#plt.cla()
#ax1.plot(x,y)
#ax2.plot(f,p,'o')
ax2.set_xscale('linear')
plt.plot(f[0:1],p[0:1],'Dm',ms=10)
plt.grid(1)
plt.show()
plt.legend(loc=0)
display(plt.gcf())
Fine¶
In [42]:
plt.figure(12)
x,y=make_signal(amp=10.,L=30.,N=200,nwaves=5.7,ystartend=(0,0),noise=0)
f,p=psd(x,y,scale=1)
plot_sig_psd(x,y*np.hanning(len(y)),scale=1.,label='5.7',c='b',ls='-')
plot_sig_psd(x,y*np.hanning(len(y)),scale=1.,label='5.7 hanning',c='r',ls='--')
plt.legend(loc=0)
Out[42]:
In [ ]:
#same test with a real signal
In [40]:
#same wave, with a tilt and offset
#%matplotlib inline
x,y=make_signal(amp=10.,L=30.,N=200,nwaves=5.,ystartend=(13,45),noise=0)
f,p=psd(x,y,scale=1)
plot_sig_psd(x,y,scale=1.)
#plt.cla()
#ax1.plot(x,y)
#ax2.plot(f,p,'o')
plt.semilogy()
plt.grid()
In [ ]:
In [10]:
In [ ]:
print np.std(y)**2
print np.trapz(p/len(y),f)
print '#--------------------'
#x,y=make_signal(amp=10.,L=30.,N=500,nwaves=5.8,ystartend=(3,20),noise=0)
plt.plot(f,p/len(x),label='p1')
plt.loglog()
plt.grid(1)
plt.legend(loc=0)
In [ ]:
f2,p2=psd(x,y)
print np.std(y)**2
print np.trapz(p/len(y),f)
plt.plot(f2,p2/len(x),label='p2')
In [67]:
In [ ]:
plt.clf()
x,y=make_signal(amp=10.,L=30.,N=500,nwaves=5.7,ystartend=(0,0),noise=0)
f,p=psd(x,y,scale=1)
plot_sig_psd(x,y,label='5.7')
x,y=make_signal(amp=10.,L=30.,N=500,nwaves=5.,ystartend=(0,0),noise=0)
f,p=psd(x,y,scale=1)
plot_sig_psd(x,y,label='5.0')
plt.legend(loc=0)
In [ ]:
np.std(y)
In [ ]:
In [ ]:
from psd2d import *
In [ ]:
outfolder=r'defect_count_test\5x'
infolder=r'defect_count_test\input_data\fits'
datafile='4MMode_5x00003.fits' #file with metrology data
fitsfile=os.path.join(infolder,datafile)
#test load images
x,y,data=getdata(fitsfile) #x and y in mm, data in um
In [ ]:
freqs,p,pl=compare_2dpsd(data,data,x,y)
In [ ]:
#normalization
plt.clf()
plt.suptitle('Signal')
plt.plot(x,y)
f,PSD=psd(x,y,scale=1)
print np.trapz(PSD,f)
plt.clf()
plt.suptitle('Effects of windowing and terms subtraction')
ax1=plt.subplot(311)
plt.title('Signal=wave+line+noise')
plt.plot(x,y)
ax2=plt.subplot(312)
plt.title('VC and RA routines')
f,PSD=psd(x,y-line(x,y),scale=1)
plt.plot(f,PSD/N/L,label='my PSD/N/L (line subtraction)')
fw,PSDw=psd(x,y*np.hanning(len(y)),scale=1)
plt.plot(fw,PSDw/N/L,label='my PSD windowed/N/L')
fw,PSDw=psd(x,(y-line(x,y))*np.hanning(len(y)),scale=1)
plt.plot(fw,PSDw/N/L,label='my PSD windowed/N/L (line subtraction)')
plt.plot(*realPSD(y,dx=L/(N-1)),label='realPSD',linestyle='-.',color='r')
plt.plot(*realPSD(y*np.hanning(len(y)),dx=L/(N-1)),label='realPSD - w',linestyle='-.',color='k')
plt.loglog()
plt.legend(loc=0)
In [ ]:
plt.ion()
outfolder=r'defect_count_test\5x'
infolder=r'defect_count_test\input_data\fits'
datafile='4MMode_5x00003.fits' #file with metrology data
fitsfile=os.path.join(infolder,datafile)
#test load images
x,y,data=getdata(fitsfile) #x and y in mm, data in um
In [ ]:
#test outliers removal
ldata=remove_outliers2d(data,fignum=1)
In [ ]:
#plot comparison of psd for the two
freqs,p,pl=compare_2dpsd(data,ldata,x,y,fignum=3,titles=['data','data outliers corrected'],
vmin=1e-9,vmax=1e-3)
plt.suptitle('2D PSD - %s'%datafile).set_size('large')
In [ ]:
#test outliers removal
ldata=remove_outliers2d(data,fignum=1)
#plot comparison of psd for the two
freqs,p,pl=compare_2dpsd(data,ldata,x,y,fignum=3,titles=['data','data outliers corrected'],
vmin=1e-9,vmax=1e-3)
plt.suptitle('2D PSD - %s'%datafile).set_size('large')
In [ ]:
#do all processing on a fits file, return freq and psd.
def avgpsd(fitsfile,correct=False):
"""Wrapper on avgpsd2d. Extract average psd from a fits image,
includes reading and outliers removal."""
x,y,data=getdata(fitsfile) #x and y in mm, data in um
if correct:
data=remove_outliers2d(data,nsigma=3)
f,p=psd2d(x,y,data)
pm=avgpsd2d(p,axis=1,scale=1e-6)
return f,pm
fg,pg=avgpsd(os.path.join(infolder,'OAB01_S03_05x_00006.fits')) #slumped glass
In [ ]:
#get average profiles for different data sets and plot them together
fg,pg=avgpsd(os.path.join(infolder,'OAB01_S03_05x_00006.fits')) #slumped glass
fgc,pgc=avgpsd(os.path.join(infolder,'OAB01_S03_05x_00006.fits'),correct=True)
fm,pm=avgpsd(os.path.join(infolder,'4MMode_5x00003.fits' )) #mandrel
fmc,pmc=avgpsd(os.path.join(infolder,'4MMode_5x00003.fits' ),correct=True)
#fw,pw=avgpsd(os.path.join(infolder,'' )) #witness
#fwc,pwc=avgpsd(os.path.join(infolder,'' ),correct=True)
plt.figure(4)
plt.clf()
plt.plot(fg,pg,label='slumped glass S03')
plt.plot(fgc,pgc,'.',label='slumped glass S03 (corr.)')
plt.plot(fm,pm,label='mandrel')
plt.plot(fmc,pmc,label='mandrel (corr.)')
#plt.plot(freqs,pm,label='mandrel')
plt.loglog()
plt.ylabel('mm^3')
plt.xlabel('mm^-1')
plt.legend(loc=0)
#plt.ylim([2e-6,1])
#plt.xlim([2e-1,4e2])
plt.grid(1)
In [ ]:
#plot different psd data of same sample (OAB mandrel)
## load data
fm,pm=avgpsd(os.path.join(infolder,'4MMode_5x00003.fits' )) #mandrel
fmc,pmc=avgpsd(os.path.join(infolder,'4MMode_5x00003.fits' ),correct=True)
fsi,psi=avgpsd(os.path.join(infolder,'OAB01_50x_00001.fits'))
fsic,psic=avgpsd(os.path.join(infolder,'OAB01_50x_00001.fits'),correct=True)
fsi2,psi2=avgpsd(os.path.join(infolder,'OAB01_50x_00002.fits'))
fsi2c,psi2c=avgpsd(os.path.join(infolder,'OAB01_50x_00002.fits'),correct=True)
In [ ]:
fsi20,psi20=avgpsd(os.path.join(infolder,'XYMode_20x00002.fits'))
fsi20c,psi20c=avgpsd(os.path.join(infolder,'XYMode_20x00002.fits'),correct=True)
In [ ]:
##plot data
plt.figure(6)
plt.clf()
plt.plot(fm,pm,label='mandrel 5x')
plt.plot(fmc,pmc,label='mandrel 5x (corr.)')
plt.plot(fsi20,psi20,label='Mandrel 20x')
#plt.plot(fsi20c,psi20c,'.',label='Mandrel 20x (corr.)')
plt.plot(fsi,psi,label='Mandrel 50x')
#plt.plot(fsic,psic,'.',label='Mandrel 50x (corr.)')
plt.plot(fsi2,psi2,label='Mandrel 50x area 2')
#plt.plot(fsi2c,psi2c,'.',label='Mandrel 50x area 2 (corr.)')
plt.loglog()
plt.ylabel('mm^3')
plt.xlabel('mm^-1')
plt.legend(loc=0)
plt.grid(1)
In [ ]:
#effects of outliers correction on 20x mandrel surface data
ff='XYMode_20x00002.fits'
x,y,d=getdata(os.path.join(infolder,ff))
dl=remove_outliers2d(d)
compare_2dpsd(d,dl,x,y,fignum=2,titles=['raw','corrected'])
plt.suptitle(ff).set_size('large')
In [ ]:
#try to add a 50x data on a different sample
fsi,psi=avgpsd(os.path.join(infolder,'Si052TR_00002.fits'),correct=True)
plt.plot(fsi,psi,label='Coated Si wafer')
plt.legend(loc=0)
In [ ]:
# Si52TR (SR Pt)
# plot different psd data of same sample
## load data
f1,p1=avgpsd(os.path.join(infolder,'Si052TR_00002.fits'))
f2,p2=avgpsd(os.path.join(infolder,'Si052TR_00003.fits'))
f3,p3=avgpsd(os.path.join(infolder,'Si052TR_00004.fits'))
f3c,p3c=avgpsd(os.path.join(infolder,'Si052TR_00004.fits'),correct=True)
In [ ]:
##plot data
plt.figure(7)
plt.clf()
plt.plot(f1,p1,label='5x')
plt.plot(f2,p2,label='20x')
plt.plot(f3,p3,label='50x')
plt.plot(f3c,p3c,label='50x corr')
plt.loglog()
plt.ylabel('mm^3')
plt.xlabel('mm^-1')
plt.legend(loc=0)
plt.suptitle('Pt coated Si wafer').set_size('large')
plt.grid(1)
In [ ]:
"""
#add outliers corrected data
f1,p1=avgpsd(os.path.join(infolder,'Si052TR_00002.fits'),correct=True)
f2,p2=avgpsd(os.path.join(infolder,'Si052TR_00003.fits'),correct=True)
f3,p3=avgpsd(os.path.join(infolder,'Si052TR_00004.fits'),correct=True)
plt.plot(f1,p1,label='5x corr')
plt.plot(f2,p2,label='20x corr')
plt.plot(f3,p3,label='50x corr')
plt.legend(loc=0)
"""
In [ ]:
x,y,gd=getdata(os.path.join(infolder,'Si052TR_00002.fits'))
p1x=p1/np.trapz(p1,f1)*np.std(gd)**2
x,y,gd=getdata(os.path.join(infolder,'Si052TR_00003.fits'))
p2x=p2/np.trapz(p2,f2)*np.std(gd)**2
x,y,gd=getdata(os.path.join(infolder,'Si052TR_00004.fits'))
p3x=p3/np.trapz(p3,f3)*np.std(gd)**2
x,y,gd=getdata(os.path.join(infolder,'Si052TR_00004.fits'))
gdl=remove_outliers2d(gd,nsigma=1.5,fignum=2)
p3cx=p3c/np.trapz(p3c,f3c)*np.std(gd)**2
#compare_2dpsd(d,dl,x,y,fignum=2,titles=['raw','corrected'])
plt.figure(8)
plt.clf()
plt.plot(f1,p1x,label='5x')
plt.plot(f2,p2x,label='20x')
plt.plot(f3,p3x,label='50x')
plt.plot(f3c,p3cx,label='50x corr')
plt.loglog()
plt.ylabel('mm^3')
plt.xlabel('mm^-1')
plt.legend(loc=0)
plt.suptitle('Pt coated Si wafer').set_size('large')
plt.grid(1)
In [ ]:
#plot PSDs of selected profiles
from pyProfile.psd import psd
#x in um and y in nm -> PSD in um^3 with scale 1e-6
from pyProfile.profile import movingaverage
for i in [1000,1100,1200]:
plt.figure('line %i'%i)
plt.clf()
z=ldata[:,i]
plot_sig_psd(x,z)
p2=movingaverage(z,31)
plot_sig_psd(x,p2)
In [ ]:
#outliers removal details on 50x
x,y,data=getdata(os.path.join(infolder,'Si052TR_00004.fits'))
ldata=remove_outliers2d(data,fignum=1)
#plot comparison of psd for the two
freqs,p,pl=compare_2dpsd(data,ldata,x,y,fignum=2,titles=['data','data outliers corrected'],
vmin=1e-11,vmax=1e-8)
plt.suptitle('2D PSD - %s'%datafile).set_size('large')
In [ ]:
fitsfile=os.path.join(infolder,datafile)
In [ ]:
fitsfile=os.path.join(infolder,'003002.fits')
#test load images
x2,y2,data2=getdata(fitsfile) #x and y in mm, data in um
ldata2=remove_outliers2d(data2,nsigma=3,fignum=3)
In [ ]:
fsi,psi=avgpsd(os.path.join(infolder,'003002.fits'))
plt.plot(fsi,psi,label='Mandrel50x')
plt.legend(loc=0)
In [ ]:
#understanding scaling factors in psd
#avgpsd does:
x,y,data=getdata(fitsfile) #x and y in mm, data in um
i=1000
#freqs,pm from avgpsd2d
axis=1
scale=1e-6
N=data.shape[0]
L=span(x,True)
freqs = np.fft.fftfreq(N,L/N)
freqs = freqs[:np.ceil(N/2)]
p=(np.abs(np.fft.fftn(data)))**2/N*L*2*scale #psd = np.abs(yfft)**2/N*L*2*scale
pm=np.mean(p,axis=axis)
pm = pm[:np.ceil(N/2)] #yfft = yfft[:np.ceil(N/2)]
#using rfft2
p=(np.abs(np.fft.rfft2(data,axes=[0])))**2/N*L*2*scale #psd = np.abs(yfft)**2/N*L*2*scale
pm2=np.mean(p,axis=axis) [1:] #exclude constant part
#pm = pm[:np.ceil(N/2)] #yfft = yfft[:np.ceil(N/2)]
#for single profile instead:
y=data[:,i]
N = len(y)
L=span(x,True)
#y=y-line(x,y) #-y.mean()
yfft = np.fft.fft(y)
yfft = yfft[:np.ceil(N/2)]
p = np.abs(yfft)**2/N*L*2*scale
plt.clf()
plt.plot(freqs,p,label='single profile')
plt.plot(freqs,pm,label='average with fftn')
plt.plot(freqs,pm2,label='average with rfft2')
plt.loglog()
plt.ylabel('mm^3')
plt.xlabel('mm^-1')
plt.legend(loc=0)
plt.ylim([2e-6,1])
plt.xlim([2e-1,4e2])
plt.grid(1)
In [ ]:
#from example above, I want to understand why fftn and rfft2 give different
# results. The answer is in the keyword axes, but why?
i=1000
#this gives a rescaled result if axes is not set in fftn (it works if set to [0]
#freqs,pm from avgpsd2d
axis=1
scale=1e-6
N=data.shape[0]
L=span(x,True)
freqs = np.fft.fftfreq(N,L/N)
freqs = freqs[:np.ceil(N/2)]
p=(np.abs(np.fft.fftn(data)))**2/N*L*2*scale #psd = np.abs(yfft)**2/N*L*2*scale
pm=np.mean(p,axis=axis)
pm = pm[:np.ceil(N/2)] #yfft = yfft[:np.ceil(N/2)]
#using rfft2
p=(np.abs(np.fft.rfft2(data,axes=[0])))**2/N*L*2*scale #psd = np.abs(yfft)**2/N*L*2*scale
pm2=np.mean(p,axis=axis) [1:] #exclude constant term
plt.figure()
plt.clf()
plt.plot(freqs,pm,label='average with fftn')
plt.plot(freqs,pm2,label='average with rfft2')
plt.loglog()
plt.ylabel('mm^3')
plt.xlabel('mm^-1')
plt.legend(loc=0)
plt.ylim([2e-6,1])
plt.xlim([2e-1,4e2])
plt.grid(1)
math:: A{kl} = \sum{m=0}^{M-1} \sum{n=0}^{N-1} a{mn}\exp\left{-2\pi i \left({mk\over M}+{nl\over N}\right)\right} \qquad k = 0, \ldots, M-1;\quad l = 0, \ldots, N-1,
In [ ]:
#explore meaning of axis
xt=np.arange(100)*2*np.pi/10
yt=np.sin(xt)
ft,pt=psd(xt,yt,scale=1.)
plt.plot(ft,pt) #this is a single peak at 0.16 for some reason
plt.figure()
plot_sig_psd(xt,yt)
#now create 2d data
yt2d=np.repeat(yt,100).reshape(100,100)
pt2d=np.fft.fftn(yt2d,axes=[0])
print pt2d.shape
plt.plot(ft,np.mean(pt,axis=1))
Fine¶
In [76]:
from pyProfile.profile import rebin_profile
%matplotlib qt
#return frequencies and PSD of a profile in mm^3,
#with x in mm and y in um with default scale value 1e-6.
# 1e-6 --> psd in um^3 if x in um and y in nm
#
i=1200
plt.figure('line %i'%i)
plt.clf()
z=data[:,i].copy()
x=dx*np.arange(len(z))
plot_sig_psd(x,z,label='raw')
#remove outliers by interpolating
mask=np.abs(z-np.mean(z))<(1.5*np.std(z)) #points to keep
z[mask==False] = np.interp(x[mask==False], x[mask], z[mask])
#---------------------
plot_sig_psd(x,z,label='interp')
p2=movingaverage(z,31)
plot_sig_psd(x,p2,label='smooth')
x3,p3=rebin_profile(x,z,bins=200)
plot_sig_psd(x3,p3,label='rebin')
plt.subplot(212)
freqs,pm=avgpsd2d(data,axis=1)
plt.plot(freqs,pm,label='2D average')
plt.ylabel('mm^3')
plt.xlabel('mm^-1')
plt.legend(loc=0)
In [ ]:
plt.figure()
plt.plot(freqs,pm,label='2D average')
plt.ylabel('mm^3')
plt.xlabel('mm^-1')
plt.loglog()
plt.legend(loc=0)
Profile reconstruction¶
In [868]:
#make a test signal
from pyProfile.psd import psd
from pyProfile.profile import line,make_signal
N=500
noise=0
amp=10.
L=30.
nwaves=5.8 #32
ystartend=(3,20)
#ystartend=(0,0)
outfolder='PSDtest'
In [773]:
%qtconsole
In [870]:
#reconstruction work perfectly.
# normalization is perfectly accurate for odd N, it fails by a few .1%
# for even N. The discrepacy is not numerically relevant, but
# suspicious. Probably the way
def test_reconstruction(x1,y1):
"""Given a profile, calculate psd and phase, use them to reconstruct
the profile and overplot it."""
#plt.clf()
f1,p1,ph1=psd(x1,y1,retall=True)
axes=plot_sig_psd(x1,y1,marker='o',label='5.0',c='r')
#rebuild profile 1
x,y=psd2prof(f1,p1,ph1,N=len(x1))
plot_sig_psd(x,y,marker='x',label='5.0',c='b')
display(plt.gcf())
print "input profile:"
print "signal rms is:",np.std(y1)
print "rms**2=",np.std(y1)**2
print "psd sum",p1[1:].sum() #np.trapz(p/len(y),f)
print '#--------------------'
print 'rms of reconstruction error:',np.std(y-y1)
def test_normalization(x1,y1):
"""Given a profile, vweify."""
#plt.clf()
f1,p1,ph1=psd(x1,y1,retall=True)
axes=plot_sig_psd(x1,y1,marker='o',label='5.0')
display(plt.gcf())
print "input profile:"
print "signal rms is:",np.std(y1)
print "rms**2=",np.std(y1)**2
print "psd sum",p1[1:].sum() #np.trapz(p/len(y),f)
print '#--------------------'
return x,y,f1,p1,ph1
plt.figure()
plt.clf()
test_normalization(*make_signal(amp=10.,L=20.,N=20,
nwaves=5.7,ystartend=(20,30),noise=0))
test_normalization(*make_signal(amp=10.,L=20.,N=20,
nwaves=5.,ystartend=(20,30),noise=3))
test_normalization(*make_signal(amp=10.,L=20.,N=21,
nwaves=5.7,ystartend=(20,30),noise=0))
test_normalization(*make_signal(amp=10.,L=20.,N=21,
nwaves=5.,ystartend=(20,30),noise=3))
Out[870]:
In [870]:
#reconstruction work perfectly.
# normalization is perfectly accurate for odd N, it fails by a few .1%
# for even N. The discrepacy is not numerically relevant, but
# suspicious. Probably the way python includes the Nywuist freq.
# for even N and the way I handle the normalization don't match.
def test_reconstruction(x1,y1):
"""Given a profile, calculate psd and phase, use them to reconstruct
the profile and overplot it."""
#plt.clf()
f1,p1,ph1=psd(x1,y1,retall=True)
axes=plot_sig_psd(x1,y1,marker='o',label='5.0',c='r')
#rebuild profile 1
x,y=psd2prof(f1,p1,ph1,N=len(x1))
plot_sig_psd(x,y,marker='x',label='5.0',c='b')
display(plt.gcf())
print "input profile:"
print "signal rms is:",np.std(y1)
print "rms**2=",np.std(y1)**2
print "psd sum",p1[1:].sum() #np.trapz(p/len(y),f)
print '#--------------------'
print 'rms of reconstruction error:',np.std(y-y1)
def test_normalization(x1,y1):
"""Given a profile, vweify."""
#plt.clf()
f1,p1,ph1=psd(x1,y1,retall=True)
axes=plot_sig_psd(x1,y1,marker='o',label='5.0')
display(plt.gcf())
print "input profile:"
print "signal rms is:",np.std(y1)
print "rms**2=",np.std(y1)**2
print "psd sum",p1[1:].sum() #np.trapz(p/len(y),f)
print '#--------------------'
return x,y,f1,p1,ph1
plt.figure()
plt.clf()
test_normalization(*make_signal(amp=10.,L=20.,N=20,
nwaves=5.7,ystartend=(20,30),noise=0))
test_normalization(*make_signal(amp=10.,L=20.,N=20,
nwaves=5.,ystartend=(20,30),noise=3))
test_normalization(*make_signal(amp=10.,L=20.,N=21,
nwaves=5.7,ystartend=(20,30),noise=0))
test_normalization(*make_signal(amp=10.,L=20.,N=21,
nwaves=5.,ystartend=(20,30),noise=3))
Out[870]:
In [871]:
#indeed the discrepacy is much larger (and numerically relevant)
# when there is a signal with Nyquist frequency N/2
plt.figure()
plt.clf()
test_normalization(*make_signal(amp=10.,L=20.,N=20,
nwaves=10.,ystartend=(20,30),noise=0))
test_normalization(*make_signal(amp=10.,L=20.,N=21,
nwaves=10.,ystartend=(20,30),noise=3))
Out[871]:
In [875]:
#but it can also be a sampling effect (e.g. aliasing, compare these
# figures with the ones above). Note that for number of points much
# larger than max freq. (e.g. N=200 nwaves=10) all four tests give
# perfect accuracy.
plt.figure()
plt.clf()
test_normalization(*make_signal(amp=10.,L=20.,N=200,
nwaves=100,ystartend=(20,30),noise=0));
test_normalization(*make_signal(amp=10.,L=20.,N=200,
nwaves=5.,ystartend=(20,30),noise=3));
test_normalization(*make_signal(amp=10.,L=20.,N=201,
nwaves=100,ystartend=(20,30),noise=0));
test_normalization(*make_signal(amp=10.,L=20.,N=201,
nwaves=5.,ystartend=(20,30),noise=3));
In [876]:
#but it can also be a sampling effect (e.g. aliasing, compare these
# figures with the ones above). Note that for number of points much
# larger than max freq. (e.g. N=200 nwaves=10) all four tests give
# perfect accuracy.
plt.figure()
plt.clf()
test_normalization(*make_signal(amp=10.,L=20.,N=200,
nwaves=100,ystartend=(0,0),noise=0));
test_normalization(*make_signal(amp=10.,L=20.,N=200,
nwaves=5.,ystartend=(0,0),noise=3));
test_normalization(*make_signal(amp=10.,L=20.,N=201,
nwaves=100,ystartend=(0,0),noise=0));
test_normalization(*make_signal(amp=10.,L=20.,N=201,
nwaves=5.,ystartend=(0,0),noise=3));
In [794]:
#debug for development of psd2prof
#compare yfft frin direct transform of profile with the one from psd and phase.
N=len(x1)
#x1,y1=make_signal(amp=10.,L=30.,N=N,nwaves=5.,ystartend=(0,0),noise=0)
#f1,p1,ph1=psd(x1,y1,retall=True)
yfft=np.fft.rfft(y1)
p=p1/2
p[0]=p[0]*2
print np.abs(yfft)
print np.sqrt(p)*N
print yfft
a= np.abs(yfft)*np.exp(1j*np.angle(yfft))
b= np.sqrt(p)*np.exp(1j*ph1)*N
print a
print b
print np.std(np.abs(a-yfft))
print np.std(np.abs(a-b))
In [794]:
#debug for development of psd2prof
#compare yfft frin direct transform of profile with the one from psd and phase.
N=len(x1)
#x1,y1=make_signal(amp=10.,L=30.,N=N,nwaves=5.,ystartend=(0,0),noise=0)
#f1,p1,ph1=psd(x1,y1,retall=True)
yfft=np.fft.rfft(y1)
p=p1/2
p[0]=p[0]*2
print np.abs(yfft)
print np.sqrt(p)*N
print yfft
a= np.abs(yfft)*np.exp(1j*np.angle(yfft))
b= np.sqrt(p)*np.exp(1j*ph1)*N
print a
print b
print np.std(np.abs(a-yfft))
print np.std(np.abs(a-b))
In [769]:
#debug for development of psd2prof
#direct and inverse transform in real case
#x1,y1=make_signal(amp=10.,L=30.,N=10,nwaves=5.7,ystartend=(20,30),noise=0)
#f1,p1,ph1=psd(x1,y1,retall=True)
yfft=np.fft.rfft(y1)
yy=np.fft.irfft(yfft,N)
print y1
print yy
p=p1/2
p[0]=p1[0]*2
p=np.sqrt(p)*np.exp(1j*ph1)*N
yyy=np.fft.irfft(p,N)
print yyy
In [729]:
N*np.sqrt(p1/2)*np.exp(1j*ph1),np.fft.rfft(y1)
Out[729]:
In [829]:
#test for normalization
plt.clf()
x1,y1=make_signal(amp=10.,L=30.,N=10,nwaves=5,ystartend=(0,0),noise=0)
f1,p1,ph1=psd(x1,y1,retall=True)
print "profile 1"
print "surface rms**2",np.std(y1)**2
print "psd sum",p1[1:].sum() #np.trapz(p/len(y),f)
print '#--------------------'
axes=plot_sig_psd(x1,y1,label='5.0',marker='o')
axes[-1].set_xscale('linear')
In [841]:
print
In [360]:
#positive freqs (including 0) are in f[:(N-1)/2+1]
# valid for both odd and even N, even if with this construction
# ir will always be odd
y=np.fft.ifft(p)
x=np.fft.fftfreq(len(y),f[1]-f[0])
#x[(N-1)/2+1:]=-x[:(N-1)/2:-1]+x[(N-1)/2]
x=np.fft.fftshift(x)
x=x-x[0]
In [ ]:
In [ ]: