import numpy as np import matplotlib.pyplot as plt from scipy import fftpack def N_for_lc (lambdac, dx): """ Parameters ---------- lambdac (float): long-wavelength roughness cutoff in mm dx (float): sampling interval in mm Returns ------- N_for_lc : number of profile points to be considered in order to get a filtered signal of length one lambdac in x """ q = int(lambdac/dx) N_for_lc = int((lambdac/dx) + (2*q)) return N_for_lc def N_for_ln (lambdac, dx, int_times_lc=5): """ Parameters ---------- lambdac (float): long-wavelength roughness cutoff in mm dx (float): sampling interval in mm int_times_lc : number of times lc is contained within ln. Default: 5 Returns ------- N_for_ln : number of profile points to be considered in order to get a filtered signal of length one ln in x """ ln = int_times_lc*lambdac q = int(lambdac/dx) N_for_ln = int((ln/dx) + (2*q)) return N_for_ln def gauss_low_pass_filter(lambdac=0.8, dx=0.001): """ Gaussian filter S in spatial domain in the interval -lc <= x <= lc (lc is lambdac) Parameters ---------- lambdac : TYPE, float DESCRIPTION. The default is 0.8. dx : TYPE, float DESCRIPTION. The default is 0.001. Returns ------- (x, S): x --> np.array of size (m,) and S --> np.array of size (m,) """ alpha = 0.4697 # Theoretical value # Generate evenly spaced (dx) in interval (-lambdac, lambdac). x = np.arange(-lambdac, lambdac, dx) # Generate the gaussian filter S = 1/(lambdac*alpha)*np.exp(-np.pi*(x/(alpha*lambdac))**2) S = S / S.sum() return (x, S) def plot_gauss_low_pass_filter(lambdac=0.8, dx=0.001): """ Plots the gaussian filter S in spatial domain in the interval -lc <= x <= lc (lc is lambdac) Args: lambdac (float): long-wavelength roughness cutoff in mm dx (float): sampling interval in mm """ # The same as for "gauss_low_pass_filter()": alpha = 0.4697 x = np.arange(-lambdac, lambdac, dx) S = 1/(lambdac*alpha)*np.exp(-np.pi*(x/(alpha*lambdac))**2) # generate the gaussian filter S = S / S.sum() # Plot the FILTER in the interval -lc <= x <= lc: plt.plot(x, S) plt.xlabel('Distance (mm)') plt.ylabel('Weighting function') plt.show() def plot_transmission_gauss(l_profile=8, lambdac=0.8, dx=0.001): """ Plots the transmission characteristics of the Gaussian filter S Args: l_profile (float): length of profile. We assume it to be 8 mm lambdac (float): long-wavelength roughness cutoff in mm dx (float): sampling interval in mm """ # The same as for "gauss_low_pass_filter()": alpha = 0.4697 x = np.arange(-lambdac, lambdac, dx) S = 1/(lambdac*alpha)*np.exp(-np.pi*(x/(alpha*lambdac))**2) # generate the gaussian filter S = S / S.sum() # Plot TRANSMISSION GAUSS: m = S.shape[0] # length of Gaussian filter n = np.int(l_profile/dx) # number of profile point # Center the filter and zero-pad to n long: S = np.pad(S,(n//2 - m//2,n//2 - m//2), 'constant', constant_values=(0,0)) Sf = fftpack.fft(S) # DFT of S j = np.arange(2,n//2+2,1) wave = n*dx/(j-1) plt.semilogx(wave, 100*np.abs(Sf[0:n//2])) plt.xlabel('Wavelength (mm)') plt.ylabel('Amplitude (%)') plt.show() def ampl_trans_sinusoid(gauss_low_pass_filter, wave_sin, l_profile=8): """ Evaluates the amplitude transmission of a sinusoid whose wavelength is within certain limits: min_wavelength <= wave_sin <= max_wavelength, where min_wavelength and max_wavelength are set by the the values (l_profile, dx) Args: gauss_low_pass_filter (x, S): x-->np.array of size (m,) and S--> np.array of size (m,) l_profile (float): length of profile. We assume it to be 8 mm wave_sin (float): wavelength of the sinusoid in mm Returns ------- (float, float): minimum and maximum aplitude transmission """ lambdac = gauss_low_pass_filter[0].max() dx = gauss_low_pass_filter[0][1] - gauss_low_pass_filter[0][0] S = gauss_low_pass_filter[1] n = np.int(l_profile/dx) # number of profile points j = np.arange(2,n//2+2,1) wave = n*dx/(j-1) # wavelength corresponding to (l_profile, dx) # Raise an exception if wave_sin is not within the allowed limits if (wave_sin < wave.min()) or (wave_sin > wave.max()): raise ValueError('Wavelength of sinusoid out of range.') else: x = np.arange(0, l_profile, dx) #print('x: ',len(x)) z = 1*np.sin(2*np.pi*x/wave_sin) # amplitude: 1 micron, wavelength: wave_sin #print('z: ',len(z)) w = np.convolve(z, S) # convolve profile and filter #print('w: ',len(w)) q = np.int(lambdac / dx) # q is half the length of the filter #print('q: ',q) w = w[np.int(q):np.int(w.shape[0] - q - lambdac)] plt.plot(x,z) plt.plot(x,w) plt.plot(x[q:len(x)-q+lambdac], w[q:len(x)-q+lambdac]) plt.show() return w[q:len(x)-q+lambdac].min(), w[q:len(x)-q+lambdac].max() def ampl_trans_profile(gauss_low_pass_filter, x, z): """ Evaluates the amplitude transmission of a profile (x,z) Args: gauss_low_pass_filter (x, S): x-->np.array of size (m,) and S--> np.array of size (m,) x (ndarray): np.array of length (n,) whose max. value is the length of the profile (given in mm) z (ndarray): profile values in mm. Returns ------- (x, z_low, z_roughness): x, z_low, z_roughness-->np.arrays of length (n-2q+1,) each. It only represents the relevant portion of the profile after removing the points from both ends. """ lambdac = gauss_low_pass_filter[0].max() dx = gauss_low_pass_filter[0][1] - gauss_low_pass_filter[0][0] #dx = x[1] - x[0] S = gauss_low_pass_filter[1] n = len(z) # number of profile points j = np.arange(2,n//2+2,1) wave = n*dx/(j-1) # wavelength array corresponding to (n, dx) w = np.convolve(z, S) # convolve profile and filter #print('w: ',len(w)) q = np.int(lambdac / dx) # q is half the length of the filter #print('q: ',q) w = w[np.int(q):np.int(w.shape[0] - q - 1)] #print(len(w) # ---> n # Get the roughness from the profile after the use of filter: z_roughness = z[q:len(x)-q+1] - w[q:len(x)-q+1] #print(len(z_roughness)) #plt.show() # Plot and show PRIMARY PROFILE and WAVINESS: # plt.figure(dpi = 600) # plt.plot(x,z) #plt.plot(x,w) # plt.plot(x[q:len(x)-q+1], w[q:len(x)-q+1]) #we only extract the relevant portion # plt.xlabel('Distancia del Perfil (mm)') # plt.ylabel('Perfil Primario (mm)') # plt.show() # Plot and show ROUGHNESS: # plt.figure(dpi = 600) # plt.plot(x[q:len(x)-q+1], z[q:len(x)-q+1] - w[q:len(x)-q+1]) # plt.xlabel('Distancia del Perfil (Parte filtrada) (mm)') # plt.ylabel('Rugosidad (mm)') # plt.show() return (x[q:len(x)-q+1], w[q:len(x)-q+1], z_roughness) #remember that the side # gauss_low_pass = gauss_low_pass_filter(lambdac=0.8, dx=0.001) # dx = 0.001 # l = 8 # n = l/dx # x = np.arange(0, l, dx) # z = np.sin(2*np.pi*x/2) # x_subset, z_filtered = ampl_trans_profile(gauss_low_pass, x, z) ############################################################################# # AMPLITUDE PARAMETERS ############################################################################# def amplitude_params_Rp_Rv_Rt(z): """ The function accepts a profile (x, z) and returns the amplitude parameters computed on the profile and evaluated over the evaluation length Parameters ---------- #x (ndarray): np.array of length (n,) whose max. value is the length of the profile (given in mm) # DESCRIPTION. z (ndarray): profile values in mm. DESCRIPTION. Returns ------- (Ra, Rq, Rp, Rv, Rt): (float, float, float, float, float) Some of the most representative amplitude parameters (American Society of Mechanical Engineers 2002) """ # n = z.shape[0] # length of profile z z = z - z.mean() # shift profile to zero mean #Ra = np.abs(z).sum() / n # compute Ra #Rq = np.sqrt(np.sum(z*z)/n) # compute Rq Rp = np.max(z) # compute Rp Rv = np.abs(np.min(z)) # compute Rv Rt = Rp + Rv # compute Rt return (Rp, Rv, Rt) def amplitude_params_Ra_Rq (z): """ The function accepts a profile (x, z) and returns the amplitude parameters computed on the profile and evaluated over the evaluation length Parameters ---------- #x (ndarray): np.array of length (n,) whose max. value is the length of the profile (given in mm) # DESCRIPTION. z (ndarray): profile values in mm. DESCRIPTION. Returns ------- (Ra, Rq): (float, float, float, float, float) Some of the most representative amplitude parameters (American Society of Mechanical Engineers 2002) """ n = z.shape[0] # length of profile z z = z - z.mean() # shift profile to zero mean Ra = np.abs(z).sum() / n # compute Ra Rq = np.sqrt(np.sum(z*z)/n) # compute Rq return (Ra, Rq) ############################################################################# # SPACING PARAMETERS ############################################################################# ############################################################################# # SURFACE FINISH PARAMETERS ############################################################################# def compute_ACF(z_roughness, dx): """ The function accepts a roughness profile and the sampling interval and returns the autocorrelation function ACF (autocovariance normalized by the square of the root mean square roughness) Parameters ---------- #dx (float): sampling interval in mm z_roughness (ndarray): roughness profile value in mm. The waviness has been previously extracted the raw profile. Returns ------- (ACV, ACF, shift): ACV --> np.array of size (n,), ACF --> np.array of size (n,), shift --> np.array of size (n,) """ n = z_roughness.shape[0] # length of profile z Rq = np.sqrt(np.sum(z_roughness*z_roughness)/n) # compute Rq ACV = np.zeros(n) for k in range(n): temp = 0 for j in range(0, n-k): # compute the sum of products temp = temp + z_roughness[j] * z_roughness[j+k] ACV[k] = temp/n # divide by number of points shift = np.linspace(0,n-1,n)*dx # generate the x axis array of shift distances ACF = ACV / (Rq)**2 return (ACV, ACF, shift) def compute_CCF(z1_roughness, z2_roughness, dx): """ The function accepts two roughness profiles and the sampling interval and returns the cross correlation function CCF (CCV(k) / Rq1*Rq2 ) Parameters ---------- #dx (float): sampling interval in mm z_roughness (ndarray): roughness profile value in mm. The waviness has been previously extracted the raw profile. Returns ------- (CCV, CCF, shift): CCV --> np.array of size (n,), CCF --> np.array of size (n,), shift --> np.array of size (n,) """ n = z1_roughness.shape[0] # length of profile z Rq1 = np.sqrt(np.sum(z1_roughness*z1_roughness)/n) # compute Rq1 Rq2 = np.sqrt(np.sum(z2_roughness*z2_roughness)/n) # compute Rq1 CCV = np.zeros(n) for k in range(n): temp = 0 for j in range(0,n-k): # compute the sum of products temp = temp + z1_roughness[j]*z2_roughness[j+k] CCV[k] = temp/n # divide by number of points shift = np.linspace(0,n-1,n)*dx # generate the x axis array of shift distances CCF = CCV / ((Rq1)*(Rq2)) return (CCV, CCF, shift) # dx = 0.001 # l = 8 # n = l/dx # x = np.arange(0, l, dx) # z = np.sin(2*np.pi*x/2) # (ACV, ACF, shift) = compute_ACF(z, dx) # plt.plot(shift, ACF)