Texture wavelength spectra



#from scipy import *

from numpy import *

from pylab import *



print '==============================='

print 'ULIP texture wavelength spectra'

print 'AdF Smit - UTA - APR 08'

print '==============================='

print ''

# 1/3rd octave bands

f13 = array([1.000, 1.250, 1.600, 2.000, 2.500, 3.150, \

4.000, 5.000, 6.300, 8.000, 10.00, 12.50, \

16.00, 20.00, 25.00, 31.50, 40.00, 50.00, \

63.00, 80.00, 100.0, 125.0, 160.0, 200.0, \

250.0, 315.0, 400.0, 500.0, 630.0, 800.0, 1000.0])

print 'Determine lower and upper limits of each 1/3 octave band ...'

b13=zeros((2,31))

b13[0,]=f13/2.**(1./6.)

b13[1,]=f13*2.**(1./6.)

print 'Read in TOT file ...'

file='DATA01.TOT'

y=load(file)

m=len(y)

print 'Data points in file:',m

y=y/1000 # convert to m

base=50 #base length for moving average

print 'Apply moving average with base length',base,'...'

print 'This is a low-pass filter to extract the macrotexture profile'

print 'from the megatexture or roughness profile'

y1=movavg(y,base)

n=len(y1)

y=y[base/2:m-base/2]

y1=y1[1:n]

y2=y1-y

print 'Run FFT ...'

y3=fft(y2)

print 'Calculate power amplitudes and wavelengths ...'

o=len(y3)

power=abs(y3[0:(o/2.)])**2.

sr=2 # cycles/mm sampling rate

p=o*1.

ftx(i)=ftx(i-1)+%SR/(2*n)

x=arange(sr/p,sr/(2.*p),sr/p)

freq=1./x

speed=2./1000. #

wave=speed/freq

plot(freq,power)

show()