Simulation of Rough Surfaces and Analysis of Roughness by MATLAB

The simulation of physical phenomena in science and engineering has become an important tool because it allows studying a wide range of real problems. On the other hand, it allows resolving problems that, because of its difficulty, it would be not possible to solve by analytical methods. Moreover, simulation is fast and versatile since it permits to vary parameters of the problem easily, allowing analyzing the effect of the modification of them in the response of the system examined. Simulation requires programming, for which there are many different languages. Each of them has a particular internal structure that distinguishes it from others. Therefore, depending on the problem to be study, it may be advisable to use a specific programming language. In the scientific-technical context MATLAB has been increasingly used by the great advantages that it offers. For example, the instructions are interpreted and not compiled, the user to enter commands interactively. The data processing is flexible. They can be read and stored in two different formats, ASCII and MATLAB format. ASCII has the advantage that the data and results may be used for other programs. However, MATLAB format may be faster. On the other hand, many functions and libraries of MATLAB are written MATLAB language, enabling the user access to the source files. It is possible to execute instructions of the operating system without exiting the program. Moreover, this language is portable in platforms as Windows or Apple, commonly employed by the researcher. From the point of view of numerical calculation, the use of matrices as basic elements makes it efficient and easy to employ, being also possible to perform graphics of curves and surfaces. Finally, the operations can be performed with simple and intuitive expressions similar to those used in science and engineering. MATLAB has been used for many applications in general physics, mathematics, optics, electronics, chemistry, biology, medicine, and artificial intelligence, among others. Now we want to employ MATLAB to simulate an optical procedure to measure surface roughness. Thus, the aim of this paper is the determination of the roughness of a surface from the analysis of the speckle pattern obtained in the far field, when the object is illuminated with a monochromatic beam perpendicularly to its surface.


Introduction
The simulation of physical phenomena in science and engineering has become an important tool because it allows studying a wide range of real problems. On the other hand, it allows resolving problems that, because of its difficulty, it would be not possible to solve by analytical methods. Moreover, simulation is fast and versatile since it permits to vary parameters of the problem easily, allowing analyzing the effect of the modification of them in the response of the system examined. Simulation requires programming, for which there are many different languages. Each of them has a particular internal structure that distinguishes it from others. Therefore, depending on the problem to be study, it may be advisable to use a specific programming language. In the scientific-technical context MATLAB has been increasingly used by the great advantages that it offers. For example, the instructions are interpreted and not compiled, the user to enter commands interactively. The data processing is flexible. They can be read and stored in two different formats, ASCII and MATLAB format. ASCII has the advantage that the data and results may be used for other programs. However, MATLAB format may be faster. On the other hand, many functions and libraries of MATLAB are written MATLAB language, enabling the user access to the source files. It is possible to execute instructions of the operating system without exiting the program. Moreover, this language is portable in platforms as Windows or Apple, commonly employed by the researcher. From the point of view of numerical calculation, the use of matrices as basic elements makes it efficient and easy to employ, being also possible to perform graphics of curves and surfaces. Finally, the operations can be performed with simple and intuitive expressions similar to those used in science and engineering. MATLAB has been used for many applications in general physics, mathematics, optics, electronics, chemistry, biology, medicine, and artificial intelligence, among others. Now we want to employ MATLAB to simulate an optical procedure to measure surface roughness. Thus, the aim of this paper is the determination of the roughness of a surface from the analysis of the speckle pattern obtained in the far field, when the object is illuminated with a monochromatic beam perpendicularly to its surface.

www.intechopen.com
The surface analysis of materials is of great importance, since many technological problems require, previously, the study of the surface state. One of the parameters of any material that changes easily with time is the roughness. Indeed, in many sectors, as civil engineering, architecture, mechanical engineering, etc. materials of different forms and properties are commonly employed, which must meet certain requirements to ensure their use. For this reason, the measurement of some surface parameters, as roughness, must be taking into consideration. There are different methods for determining roughness. One of the most employed is the profilometer (see next section). However, this paper deals with an optical method based on the speckle interferometry which has some advantages. The methodology is fast, accurate, and does not contact the specimen. Above all this we will talk in the next sections

Discretization of the problem. Roughness
A classic device for measuring surface roughness is the mechanical profilometer which is formed by a tiny stylus (with a small ball), and a displacement sensor. The typ moves along a straight line parallel to the surface plane and records the displacements in the perpendicular direction, tracing out the outline of the surface. If the ball has a diameter bd it can not be inserted between two grooves whose distance is less than bd, being only possible to detect the topographic level with a distance between grooves greater than bd.
z=h(x,y) X bd Fig. 1. Classical device for measuring roughness. Observe that when the diameter of the end needle bd is greater than the groove, the transducer can not reproduce the high frequencies of the surface outline.
In this article we are interested in using a speckle technique to measure the roughness of a surface. From a didactic viewpoint, the explained idea of the profilometer may be employed to understand the sampling, when a rough surface is simulated by MATLAB. To start let us suppose a one-dimensional rough surface, and then we extend the results to the case of two variables. The height of the rough surface can also be measured by sampling. With this aim let us consider a curve z=h(x) as shown in Fig. (2), aligned on the OX axis. For sampling the function h(x) the X co-ordinate axis is divided into intervals of length u measured with respect the origin O resulting in a system of aligned points of co-ordinates 0, u, 2u, 3u,… The distance u between two neighbor points, i.e. the sample interval, is called the sampling period and its value is chosen depending of the function to be investigated (in our case the form of the surface). The distance u between two points may be likened to the ball diameter bd of the profilometer. If N samples are taken, they form a string of N integers for which a value z p =h(xp) is given. This set of numbers is collecting in a matrix IF of dimension N × 1. The range of variation of the index p is 1, 2, 3,..., N, and it represents the element p of the string. Therefore this index p is related with the sampling period as follows: x = 0, u, 2u, 3u, ..., (p-1)u. Two neighboring elements of the IF matrix contain the values of the surface heights of the grooves of two points on the reference plane separated u meters. As it will see, when studying the phenomenon of diffraction in the far field, the Fourier transformation must be applied. Therefore, we need to study also the sampling in the frequency domain. When calculating the finite discrete Fourier transform (DFT) of the IF matrix of N elements, a new set of N numbers is obtained which is grouped in another matrix FO of dimension N × 1. Due to the Fourier transform is performed from the discrete values of IF, the result is also discrete. As a result the distance between two points of the transformed numbers in frequency domain is also quantified. Denoting by ν (1,2,…N) the index for the matrix FO, the row index represents the harmonic components whose frequencies are α= 0, 1 / (Nu), 2 / (Nu), ..., (ν-1) / (Nu). The sampling frequency is defined as f s ≡ 1/u, measured in m -1 , and represents the number of measures per unit length. By using this expression, spatial frequency components may also be written in the form α = 0, f s / N, 2 f s / N, ..., (ν-1) f s /N . In general, the matrix element ν of FO represents the harmonic ν α in the space of (1) Thus there is a correspondence between the index ν = 1, 2, 3, ..., N, and the spatial frequency α by means of the factor (f s /N). Obviously, the sampling process implies that some information about the sampled function is loosed, because no value between two neighboring points is known. However, depending on the physical problem studied, using the Shannon theorem, the sample interval can be modified so that the data be enough for numerical calculus.
Taking into account that the expression obtained for representing a point on the OX-axis has the form (p-1)u, it seems to be appropriate to change to non-dimensional variables. To do so we divide x for u resulting (x/u)=(p-1), p being 1,2,…,N. This new variable represents the distance from the origin O to an arbitrary point on OX (Fig.3  By setting α(Nu)= α/(1/Nu) =(ν-1) the spatial frequency α is converted in a dimensionless number representing the basic unities of measurement in this space. Furthermore, α/(1/Nu) is directly the frequency measured in unities of (1/Nu) corresponding to the element ν of the matrix FO (Fig.4). When generating a rough surface the components of IF are real numbers (Fig.3). However, as we will see, the diffraction of a light beam by a surface can introduce phase factors resulting in complex numbers in the elements of matrix IF. In any case the resulting N × 1 string of IF and its fast Fourier transform (FFT), are calculated without difficulty with MATLAB. One advantage of the aforementioned procedure is that the sampling distance u between two points of IF is not directly involved in the numerical calculation, and then it may be considered as a parameter.  (Fig. 3). Two points of coordinates (x i , y j ) and (x k , y l ) respectively are separated in the matrix IF a distance ((i-j) 2 + (k-l) 2 ) 1/2 u, and in dimensionless co-ordinates ((i-j) 2 +(k-l) 2 ) 1/2 . In relation to FO similar expressions may be obtained, but in frequency space. So the spatial frequencies between two points whose coordinates are ( ) is ((h-k) 2 + (l-m) 2 ) 1/2 (1/Nu) in m -1 , and without dimensions ((h-k) 2 + (l-m) 2 ) 1/2 .

Fraunhofer diffraction with MATLAB
In this section we are interested in the phenomenon of diffraction of light, given the importance to understand the speckle patterns. With this objective let us use the experimental lay-out depicted in Fig.(5). A collimated monochromatic laser beam LB of wavelength λ is directed to a beam splitter BS, which projects the light perpendicularly on a diffracting rough sample S located on the OXY plane. The surface has, in principle, a variable reflectance R(x,y). This means the different scatters that form the surface may have distinct reflection properties. The shape of the radiation beam used determines the geometry and the intensity inside of the illuminated area. If we suppose a beam of homogeneous intensity, its geometry can be expressed easily (in view of the simulation) as an opaque mask M placed on the surface, which has the function to define the illuminated area (Fig. 6). Taking into consideration the most cases studied in optics, we will choose a circular mask of diameter D. An observation screen is placed parallel to the diffracting surface at a distance of z from its plane. The points on the observation plane are specified by means of its x', y' coordinates, with respect to an O'X'Y' coordinate system (on the CCD camera). Supposing that the scalar diffraction theory applies, the Fresnel-Kirchhoff integral and the theories of Rayleigh-Sommerfeld can be used. However, the calculation of the diffraction pattern through these theories is not always easy to carry out. Sometimes the procedure may be simplified under certain conditions of the problem. So if the linear dimensions of the aperture (mask in our case) is much greater than the wavelength, i.e. D>>λ, and the distance z between the surface and the observation plane is great enough, the paraxial theories apply. In this case, the mathematical expression for the diffracted field depends on the specific www.intechopen.com dependence between D and z. When expanding the phase term in the Fresnell-Kirchhoff integral is not possible to neglect the quadratic terms that appear we speaks of Fresnel diffraction. On the contrary, if these terms can not be tacked into consideration we have Fraunhofer diffraction. These approximations are the most important cases in the field of the classical optics. A possible quantitative criterion to be employed in order to use the Fraunhofer approximation, or that of Fresnel, is based on taking a circle of diameter D, which only includes the regions of interest (in the present case the hole of the mask). Let r be the distance from a point on the diffracting surface to the observation point. Let ρ be the distance from the centre of the circle to a point inside its circle. If 2πr/λ varies linearly with ρ, the diffraction is called Fraunhofer diffraction; if the variation has non-linear terms of magnitude comparable with π/2, the diffraction is said Fresnel diffraction. Therefore, for Fraunhofer diffraction we obtain z>>D 2 /(4λ). In short, the diffracting area must be greater than λ and the observation of the intensity pattern must be carried out from a large distance with respect to the scatter surface. In other circumstances, i.e. if the distance z does not fulfil the conditions needed, non-paraxial terms of the phase must be included in the integrand of the Fresnell-Kirchhoff integral (higher expansion coefficients). Fraunhofer diffraction is related with the Fourier transform which takes an angular spectrum of the reflectance (or transmittance) to be considered. From a physical point of view it is equivalent to observe the phenomenon in the far field (another possibility is to employ a lens and locate the observation plane on its back focal plane). This angular spectrum means that the Fraunhofer diffraction gives the behaviour of the field amplitude for the directions in space. If we use two variables, the amplitude of the diffracted field done through the Fourier transform depends on α and β, which are related with the directions ( ) (2) As we will see in the following section, the proposed method for measuring roughness is developed under the supposition that the conditions of the Fraunhofer diffraction apply. Therefore, this case must be translated to the context of MATLAB. With this aim, the basic results of the preceding section should used. The elements fo ij of the matrix FO belonging to a row or column represent the complex amplitude of two harmonics separated 1/(Nu). Therefore the first angular direction is , which corresponds to the frequency 0 = α and the direction for the least coefficient of FO is θ x = cos -1 [(N-1)λf s /N] corresponding to the higher frequency α= (N-1) f s /N. In the case of nondimensional variables we can use for the two axes Nu cos θ x /λ, and Nu cos y θ /λ, respectively. If the diffraction pattern is observed on a plane screen a distance z from the diffusing surface, the spatial frequencies may be related with points on that plane.
For small angles θ it can be written: www.intechopen.com and cos ' Due to the properties of the Fourier transform, the FFT of the reflectance will contain N/2 of positive frequencies, and N/2 negative, whose zero spatial frequency occurs at ν =1. In the FFT, the independent variable is the frequency, and in the representation with positive and negative frequencies its maximum value will be f s /2. Based on a reflectance matrix of NxN elements located at the XOY axes associated, practically centred in the middle, we calculate FO by means of the FFT, obtaining another matrix from the centre of which the amplitude of the null frequency harmonic component is indicated.
The intensity registered over a direction ( )

Speckle pattern generation
When a laser beam illuminates a rough surface at scale of the wavelength, the diffraction pattern consists of a random distribution of intensity called speckle. The apparition of speckle may be understood by the fact that the coherent waves falling on the rough surface travel a different optical path from the diffusing surface to the observation point. When the object is rough, the reflectance is a random function on the aperture, and then the corresponding optical paths for the different scatters vary rapidly. As a result, the intensity on the observation screen (or space) also varies very quickly from one point to another of its surroundings, giving brilliant and dark spots irregular in shape. A model of diffusing and non-absorbent surface is proposed, in which the height of the scatters with respect to a reference plane are supposed as a random variable, and with a gaussian probability density function. A surface of these characteristics is, for example, a metal which is not well polished. We suppose that the rough surface is illuminated by a collimated light beam perpendicularly to its plane resulting in a speckle pattern which is calculated by means of the FFT (Fig. 7).
Due to that optical path δ followed by the different points of the wavefront is not the same, consequently, neither is the phase 2πδλ. As we have to count the return path, the path length and height h(x,y) of the surface referred to the plane z=0 are related by the expression 4( , ) hxy π λ . Thus, the reflectance will be proportional to the exponential of this phase factor, adopting the form where R 0 (x,y) is the reflection coefficient of the surface, and ( ) 4( , ) ih x y π λ is the phase. In the simulation presented in this paper we choose R 0 (x,y)=1. By measuring the random height h(x,y) of the sampled points (Fig.8), it would be possible to construct the reflectance matrix of N × N elements. Following the nomenclature of the www.intechopen.com  www.intechopen.com preceding section we call this array the IF matrix (Fig.9). The elements of IF contain the complex reflectance R(x,y) corresponding to each point of the surface, which are separated from their neighbors a distance equal to the sampling period u. The area of the delimiting mask will be represented by points outside a circle with zero reflectance. Fig. 9. View of the grid chosen on the OXY reference plane for N=64. The colours represent the surface heights at each point (pixel).
Following the same way as in preceding paragraphs, we employ the ratio h/λ as a nondimensional variable, which will be very useful when changing the wavelength. In the model this variable is equal to a constant multiplied by a random number, which will provide information on the roughness in the simulation. We will call in the program this constant RU and it represents a roughness modulating factor. Random numbers with Gaussian distribution are generated in MATLAB by the command randn. The mathematical expression for reflectance is and the phase () www.intechopen.com Thus, an element of RUrandn is a number equal to an optical path measured in wavelengths. For example, an RU=1 and a randn=2 give rise max(RUrandn)=2, which indicates a maximum path difference of 2λ, that is to say, a groove on the reflecting surface with depth equal to λ. However with the same randn, but with the modulating factor equal to 0.1, the roughness would be a tenth part. Hence the RU factor represents the roughness measured in wavelengths.
To account the transversal geometry of the incident laser beam on the surface, the rough surface is delimited by means of a round mask of diameter D (geometry could be different; see section 7.2). The diameter D must be greater than the wavelength and the sampling period u. On the other hand it is supposed that the number of sampled points inside the diameter D is large enough, in order be sure that the statistics applies.
Once that the characteristics of the surface and beam are defined, the diffraction pattern is obtained by means of the FFT of the reflectance matrix IF. The registered intensity of the diffracted light by the rough surface is proportional to the square modulus of the diffracted amplitude, e.g.

Definitions of roughness
In this section we try to adapt some definitions of roughness to our specific problem. We start the quantitative definition of the average roughness R a from the mean surface level, as the average absolute value of the height, for all the points along a straight line (remember the profilometer). Then in a circular matrix of diameter D inside the IF, corresponds 2 4 BD π elements. Therefore, the roughness of the sample may be expressed by the following formulae where the sum is extended to the sampled points within the circle of diameter BD. As previously, if we transform this Eq.(8) to non-dimensional variables, we get The number of elements G within the beam of diameter D (BD) is less than the N × N elements of IF. Say L the length of the square side where the surface is defined. In any case BD L Similarly, the roughness R q (root mean square) could be expressed as function of BD. In fact, considering the usual definition of this parameter, the following formulae may be written and its non-dimensional value

Programming using MATLAB
We will see that the simulated specklegram corresponding to the diffraction of a monochromatic radiation by a rough surface is altered by the roughness of the object within a certain range, which depends on the wavelength of the beam used. Therefore, by analyzing some characteristics of the intensity pattern it would be possible to measure roughness.
To understand the idea let us suppose a flat surface, well polished, delimited by an aperture (mask). If a beam strikes on the surface, the delimiting aperture diffracts it resulting in an intensity pattern that depends on the geometry of the obstacle. Now if the surface is scratched, the intensity registered changes, although the aperture maintains its geometry. In both cases the autoconvolution of the intensity is different, which means that the roughness produced on the surface is the cause of the change. Therefore, the convolution of the diffraction pattern could be indicative of the degree of surface polish.
To test the hypothesis, first we constructed a computer model of a rough surface, and second we simulate the diffraction of a collimated monochromatic beam by this surface. The resulting random intensity, that is, the speckle, is stored in a matrix (FO) and its autoconvolution (CO) is performed. Once all data of CO are obtained, the functional relationship of the maximum value of the autoconvolution and its relation with the roughness is analyzed. The program consists of the following steps: 1. Begin by setting the number of samples N along each axis. 2. The matrix IF is constructed by using the command RAN = randn (N). 3. The diameter of the laser beam BD is specified, measured in number of array elements. 4. A value to the RU is assigned. 5. The BS array is constructed. The mask is 0 outside the circle and 1 inside. 6. The matrix RURAN = RU * RAN is introduced, representing the surface heights for each pixel on the area NxN. 7. The matrix hs is defined as hs = RURAN.*BS. It represents the height of the points inside the circle (mask M). 8. The reflectance matrix is obtained. Its expression is ts = exp (4πi RURAN). 9. The array FO is calculated, which is the FFT of ts.
10. The intensity of the diffraction pattern is determined , 2 FIDI FO = .
www.intechopen.com 11. The autoconvolution CO of FIDI, and its maximum COV is computed. 12. In order to manipulate the data more easily, the logarithm of COV is given (log(COV)). The detailed program may be found in appendix A On the other hand, the first minimum given in the theory of diffraction by a circular hole is Both results agree and differ in a small amount. The difference can be attributed to the small number of values chosen. The second row refers to the same mirror, but not completely polished, and with a coefficient RU = 0.1. The profile shows small heights and valleys. The Airy disk is a little blurred, and not as clear as in the previous case. In the third row RU = 0.2 the central disk appears deformed and a speckleled. In the fourth and fifth rows the figure is quit different with respect to the first one, and the speckles are on all the pattern. In the last row only speckle may be seen, and no traces of the Airy disk are present. When the roughness is RU = 0.5 (bottom row) yields a rough surface with high grooves. The intensity is formed by irregular random spots being unknown directly the form of the mask, e.g. the symmetrical intensity circle of the Airy function. In summary, diffraction by a specular surface delimited by an aperture produces an intensity pattern concentrated around the direction of the reflected beam, but if the roughness is increasing, the light is diffracted producing speckle which structure is random. The third figure of each row (third column) corresponds to autoconvolución (CO), which has a maximum at the center (COV). In effect, the values for the logarithm of COV are, respectively: logCOV(RU=0.0) = 7.78, logCOV(RU=0.1) = 7.10, logCOV(RU=0.2) = 6.86, logCOV(RU=0.3)=6.90, logCOV(RU=0.4)= 6.88, and logCOV(RU=0.5)=6.91. In this calculus the logarithm of the autoconvolution hass been used because the maximum value of CO is very large. Employing log(COV), the data are easier to manipulate. These results show that the values of their maxima are not the same. On the contrary, the maximum value for each one depends on the surface roughness. For this reason it seems suitable to employ the maximum value of the autoconvolution of the speckle pattern, as a possible procedure for measuring the roughness of a surface, if the roughness is smaller than the wavelength of light used in the experiment. At the same time, COV depends on the diameter of the beam used (BD) also. To see the effect in the autoconvolution when the wide of the laser is changed, we computed logCOV with N and D for two different number of data and diameters. For example if N= 64, and BD = 32, it yields logCOV = 12. 0, 11.7, 10.7, 9.8, 9.7, 9.7, whereas  From these results may be inferred that if the area of the illuminated surface is known, measuring experimentally the autoconvolution of the speckle pattern, it would be possible measuring the roughness of this surface. To conduct laboratory experiments would be necessary to build larger tables with more values, for different incident beam intensities. By examining the calculations it also follows that for values of RU close to zero, the difference between logCOV, corresponding to a BD, and a diameter corresponding to half value, is approximately constant and equal to 2, i.e.

Circular beam
This property will be important for ulterior calculations.

Square mask
In the preceding developing calculations, a circular geometry for the beam was supposed. However, other possibilities may occur. For instance, when a laser ray is directed onto a sample under an angle of incidence θ , the effective area intersected by the beam has a quasielliptical form. Although an elliptical mask is easy to simulate with MATLAB, this paragraph deals with the study of the effect of employing a square aperture. This approximation simplifies the program, since there is no need the beam diameter datum. Moreover, from the point of view of the results, it has little influence in the final values when comparing these values with those obtained for an elliptical mask. The simulation gives the results of log(COV) for RU = 0.0, 0.1, ... 0.5, and N = 16, 32, 64, 128, 256, 512, that appear in the following  Table 2. Values of the autoconvolution log(COV) for different data and roughness parameter RU. The results in green do not give information. Figure 11 represents the values of the attached table II. These curves show the dependence of log(COV) with the roughness for different values of N, provide that the roughness is less than 0.4. Therefore, the trend is maintained even if the aperture is different. From the figure it follows that, except for values marked in green on the table, the dependence of logCOV with roughness is approximately parabolic, and can be approximated by the equation As in section 7.1., the difference of the log(COV) for consecutive values of RU, follows certain regularity. In fact, if the values of log(COV) for RU=0 are examined (see figure 11), we observe that for adjacent values of this variable, the differences between two consecutive points (corresponding to double N) are: 2.41, 2.41,2.41, 2.41, 2.40. Taking into consideration these differences, the following mathematical relationship is verified: www.intechopen.com The advantage of this formula is that it allows calculating the value of the roughness for each N and D.
To It may be seen that the differences between these values calculated with the formula (20) and those displayed in table II are equal to or less than 0.02, except for the case N = 32, RU = 0.3, which is 0.15. But as in this case the value of the table does not correspond to the difference of logarithms (marked in green), it follows that the equation obtained is suitable for the specified intervals.    A well-known statement says that the PID controller is the â€oebread and butterâ€ of the control engineer. This is indeed true, from a scientific standpoint. However, nowadays, in the era of computer science, when the paper and pencil have been replaced by the keyboard and the display of computers, one may equally say that MATLAB is the â€oebreadâ€ in the above statement. MATLAB has became a de facto tool for the modern system engineer. This book is written for both engineering students, as well as for practicing engineers. The wide range of applications in which MATLAB is the working framework, shows that it is a powerful, comprehensive and easy-to-use environment for performing technical computations. The book includes various excellent applications in which MATLAB is employed: from pure algebraic computations to data acquisition in real-life experiments, from control strategies to image processing algorithms, from graphical user interface design for educational purposes to Simulink embedded systems.

How to reference
In order to correctly reference this scholarly work, feel free to copy and paste the following: F. Gascoń and F. Salazar (2011). Simulation of Rough Surfaces and Analysis of Roughness by MATLAB, MATLAB -A Ubiquitous Tool for the Practical Engineer, Prof. Clara Ionescu (Ed.), ISBN: 978-953-307-907-3, InTech, Available from: http://www.intechopen.com/books/matlab-a-ubiquitous-tool-for-the-practicalengineer/simulation-of-rough-surfaces-and-analysis-of-roughness-by-matlab