File Exchange

image thumbnail

1-dimensional surface roughness power spectrum of a profile or topography

version (2.98 KB) by Mona Mahboob Kanafi
Calculates 1D surface roughness power spectrum for a surface topography/or multiple line profiles


Updated 19 Dec 2016

View Version History

View License

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:


*** 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

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

Cite As

Mona Mahboob Kanafi (2021). 1-dimensional surface roughness power spectrum of a profile or topography (, MATLAB Central File Exchange. Retrieved .

Comments and Ratings (22)

Khaustov Alexandr

Elijah Uche

Hi 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.

Sali SE

Hi 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?


Sali SE

Hi 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 for the code.

Tae-Yong Do

Dennis Ernens

To 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.


Hi 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ña

Many 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)
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ña

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.

Teron Nguyen

The 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 Nguyen

Dennis Ernens

sandra arias

The code is not working. Some reason I got the following warning: Negative data ignored

Thurian Le Dû

Mona Mahboob Kanafi

Hi 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 Broom


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


Mona Mahboob Kanafi

Hi 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.


Can 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


Thanks! Works very well.


MATLAB Release Compatibility
Created with R2016a
Compatible with any release
Platform Compatibility
Windows macOS Linux

Community Treasure Hunt

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

Start Hunting!