version 1.6.1.0 (2.98 KB) by
Mona Mahboob Kanafi

Calculates 1D surface roughness power spectrum for a surface topography/or multiple line profiles

In surface roughness analysis, one of the powerful tools for roughness characterization is surface roughness power spectrum (PSD). If you have multiple line scans, for examples obtained by stylus profilers, or, if you have a 3D surface topography in which you are interested to calculate power spectrum in x direction and/or in y direction (e.g. to check anisotropy of your surface profile), you need to calculate surface roughness power spectrum for each line profile and then average calculated power spectrums together.

With this function you can calculate 1D PSD of a surface topography or multiple line scans stored in a matrix.

In case of 1D line profiles, store all your lines in one matrix, i.e. if you have ‘n’ measurements where each measurement consists of ‘m’ points, you must first store all your measurements in a matrix with size n*m, called here z, and just use below line to get your calculated 1D PSD for you:

[q , C] = psd_1D(z , PixelWidth , 'x') ;

Be aware that, you need to know you PixelWidth (spatial resolution) which is obtained easily by dividing you measurement length (L) to the number of pixels/points in length (n), i.e. :

PixelWidth = L / n; % in SI units

In case you have a surface topography matrix of size n*m, called here z, you can calculate 1D PSD considering profiles parallel to x axis or y axis. Just use below lines:

[qx , Cx] = psd_1D(z , PixelWidth , 'x') ; % profiles will be z(1,:), z(2,:), ... and so on.

[qy , Cy] = psd_1D(z , PixelWidth , 'y') ; % profiles will be z(:,1), z(:,2), ... and so on.

In order to plot the output:

loglog(q,C)

*** The "Warning: Negative data ignored " is normal and is due to shifting the frequency components to have symmetry with respect to zero frequency.

For 2-dimensional surface roughness power spectrum on surface topographies, refer to:

Radially Averaged Surface Roughness Power Spectrum

(PSD) http://se.mathworks.com/matlabcentral/fileexchange/54297-radially-averaged-surface-roughness-power-spectrum--psd-

Radially averaged surface roughness power spectrum (PSD) only on top or bottom part of a topography

http://se.mathworks.com/matlabcentral/fileexchange/60771-radially-averaged-surface-roughness-power-spectrum--psd--only-on-top-or-bottom-part-of-a-topography

To generate an artificial surface topography to try the code, refer to:

Surface generator: artificial randomly rough surfaces

http://se.mathworks.com/matlabcentral/fileexchange/60817-surface-generator--artificial-randomly-rough-surfaces

Mona Mahboob Kanafi (2021). 1-dimensional surface roughness power spectrum of a profile or topography (https://www.mathworks.com/matlabcentral/fileexchange/54315-1-dimensional-surface-roughness-power-spectrum-of-a-profile-or-topography), MATLAB Central File Exchange. Retrieved .

Created with
R2016a

Compatible with any release

**Inspired by:**
Radially averaged surface roughness/topography power spectrum (PSD), Radially averaged surface roughness power spectrum (PSD) only on top or bottom part of a topography, Surface generator: artificial randomly rough surfaces

**Inspired:**
Radially averaged surface roughness/topography power spectrum (PSD), Surface generator: artificial randomly rough surfaces, Radially averaged surface roughness power spectrum (PSD) only on top or bottom part of a topography

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!Create scripts with code, output, and formatted text in a single executable document.

Khaustov AlexandrElijah UcheHi Mona,

Please I would like to be able to obtain the underlying equations and maths on which this rough surface model was built on so I could fully appreciate and understand how it works.

Thanks

Sali SEHi Mona,

thanks for the code, it helps a lot. however, the results I get using this code is different than what I get from a software designed for this! is your code doing the PSD in uniform window?

Thanks,

Sali

Sali SEHi Mona,

thanks for the code. I have a problem with this which is a bit strange. I tried a couple of different surfaces but it always give me the same results! do you know what is wrong?

Thanks,

Z

fabingJairan Nafar DastgerdiThanks for the code.

Tae-Yong DoDennis ErnensTo check the correctness of the function output I generated pure sinusoidal surfaces. Similar to Wade (entry below) I discovered that there is a 2*pi scaling error for the frequency or wavelength. Hence dividing by 2*pi places the power peak at the correct wavelength location for the generated sine wave. For more details the paper of Persson et al. is cited in the code: Persson BNJ, Albohr O, Tartaglino U, Volokitin AI, Tosatti E. On the nature of surface roughness with application to contact mechanics, sealing, rubber friction and adhesion. Journal of Physics Condensed Matter : An Institute of Physics Journal 2005;17:R1–62. doi:10.1088/0953-8984/17/1/R01.

After reading the paper it is still not clear to me why q is scaled by 2*pi and I am probably missing something obvious. It might have to do with how Persson defines his upper and lower bounds for q as 2*pi/a < q < 2*pi/L. These bounds correspond to what can be reliably represented on the domain size (remember Nyquist). See figure 6 for a graphical representation. Thus the bounds will always be shifted by 2*pi and if one wants to compare to those you need to scale the whole q-vector.

Taking the example in section 3.4 for the asphalt surface, q0 = 1000 m^-1 when scaling with 2*pi. This corresponds to the largest stones of 1 cm. Which I would have expected to give a q0 = 100 m^1. So, I think the conclusion on the scaling is right. For now I will just divide by 2*pi to have a direct connect to the real specimen, instead of needing to keep in mind the 2*pi. However, please keep in mind the limitations that the domain size puts on the spectral content that can be represented.

wadeHi Mona, can you show me the formula of PSD in your matlab code? When I use the code, I find that the x and y values in my results are always larger than AFM results by several times. I find the definition of frequency. In your definition, f=2*pi/lambda while frequency is defined by f=1/lambda in AFM software. When I divide x values (f), obtained by your code, by 2*pi, the x values of two results agrees well. I just change the x definition provided that PSD definition and values are same in the 2 methods. Now I would like to see the formula of PSD in your matlab code to check whether there are same in 2 methods.

In addition, my result is very small (1e-32) if the unit is SI. I would like to confirm whether the common computer can calculate the small values. Thanks a lot!

Juan PeñaMany thanks sharing for the code.

Just a couple of questions:

is the q expressed in (rads/m)?

is the vertical axis in G(n) or in G ( omega) as expressed in ISO 8608.

just asking because I saw in this thread* that X axis is given in (1/m)

*https://es.mathworks.com/matlabcentral/answers/262428-how-can-i-get-this-pattern#answer_205019

many thanks

-----------------------------------

To teron:

-You dont have to worry by negative data ignored, if you plot you r data you will see that you have same values at both sides of the Y axis.

-I guess you have to do a analysis for each line at a time.

Juan PeñaTo teron:

-You dont have to worry by negative data ignored, if you plot you r data you will see that you have same values at both sides of the Y axis.

-I guess you have to do a analysis for each line at a time.

Teron NguyenThe code is working, but there is also warning about negative data ignored.

I have a z matrix of left and right profiles (size = 2xm), will the final result be the average PSD of 2 profiles?

Thank you.

Teron NguyenDennis Ernenssandra ariasThe code is not working. Some reason I got the following warning: Negative data ignored

Thurian Le DûMona Mahboob KanafiHi Matheu,

You need to provide 3 inputs for the code to work. First is your height profile (z), second is the lateral distance between each of your data points. And the third one, is just to inform whether you have stored your 1D z profiles column-wise or row-wise in your z matrix. If for instance, you stored them row-wise in your z matrix, just use the function like this:

[qx , Cx] = psd_1D(z , PixelWidth , 'x') ;

Let me know if this isn't what you're looking for.

matheu BroomHi

Sorry cant seem to get this to work. How do you define the dim? I have tried everything

Matheu

Mona Mahboob KanafiHi Xaris,

If you are looking for road surface roughness features in PSD like Hurst exponent,... you need to convert your time signal into displacement signal. I assume while getting your signal, you have recorded the velocity as well. So, you can get the surface roughness profile with it and then use the code as instructed.

If you still are interested in getting some frequency features from your time signal, you can use the code, but you need to change some coefficients in Cq. Let me know if the the latter is what you are looking for.

BobCan I use this code for PSD Road Roughness?

I mean I have Displacement vs Time which is a Road Profile and I want to take the PSD

Eric_schThanks! Works very well.

Eric_sch