# Surface roughness analysis from raw data

30 views (last 30 days)

Show older comments

##### 0 Comments

### Answers (1)

Mathieu NOE
on 20 Jun 2022

hello

this is a starter , to do the spectral analysis of your data

then you can use some low / high or bandpass filters to filter your analog signal according to your cut off frequencies (and filter order)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% load signal

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

filename = ('nominalprofile.txt');

data = readmatrix(filename); % 2 channels : Time + signal

time = data(:,1);

dt = mean(diff(time));

signal = data(:,2); % select here one or several channels

Fs = 1/dt; % sampling rate

[samples,channels] = size(signal);

% time vector

t = (0:samples-1)*dt;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% define filters (optionnal)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% low pass filter %

option_LPF = 0; % 0 = without low pass filter , 1 = with low pass filter

fc_lpf = 100; % LPF cut off freq

N_lpf = 5; % filter order

% high pass filter %

option_HPF = 0; % 0 = without high pass filter , 1 = with high pass filter

fc_hpf = 100; % HPF cut off freq

N_hpf = 5; % filter order

% bandpass filter %

option_BPF = 0; % 0 = without bandpass filter , 1 = with bandpass filter

f_low = 50; % lower cut off freq Hz

f_high = 500; % higher cut off freq Hz

N_bpf = 5; % filter order

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% FFT parameters

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

NFFT = 1000; %

OVERLAP = 0.75;

%% apply filters

% low pass filter section %%%%%%

if option_LPF ~= 0

w0_lpf = 2*fc_lpf/Fs;

% digital notch (IIR)

[b,a] = butter(N_lpf,w0_lpf);

% now let's filter the signal

signal = filter(b,a,signal);

end

% high pass filter section %%%%%%

if option_HPF ~= 0

w0_hpf = 2*fc_hpf/Fs;

% digital notch (IIR)

[b,a] = butter(N_hpf,w0_hpf,'high');

% now let's filter the signal

signal = filter(b,a,signal);

end

% band pass filter section %%%%%%

if option_BPF ~= 0

[b,a] = butter(N_bpf,2/Fs*[f_low f_high]);

signal = filter(b,a,signal);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% display 1 : time domain plot

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(1),plot(time,signal);grid on

title(['Time plot / Fs = ' num2str(Fs) ' Hz ']);

xlabel(' length (milimeter) ');ylabel('roughness value (micrometer)');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% display 2 : averaged FFT spectrum

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[freq, sensor_spectrum] = myfft_peak(signal,Fs,NFFT,OVERLAP);

figure(2),semilogy(freq,sensor_spectrum);grid on

df = freq(2)-freq(1); % frequency resolution

title(['Averaged FFT Spectrum / Fs = ' num2str(round(Fs)) ' Hz / Delta f = ' num2str(df,3) ' Hz ']);

xlabel('1/length (1/mm)');ylabel('roughness value (micrometer)');

function [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)

% FFT peak spectrum of signal (example sinus amplitude 1 = 0 dB after fft).

% Linear averaging

% signal - input signal,

% Fs - Sampling frequency (Hz).

% nfft - FFT window size

% Overlap - buffer percentage of overlap % (between 0 and 0.95)

[samples,channels] = size(signal);

% fill signal with zeros if its length is lower than nfft

if samples<nfft

s_tmp = zeros(nfft,channels);

s_tmp((1:samples),:) = signal;

signal = s_tmp;

samples = nfft;

end

% window : hanning

window = hanning(nfft);

window = window(:);

% compute fft with overlap

offset = fix((1-Overlap)*nfft);

spectnum = 1+ fix((samples-nfft)/offset); % Number of windows

% % for info is equivalent to :

% noverlap = Overlap*nfft;

% spectnum = fix((samples-noverlap)/(nfft-noverlap)); % Number of windows

% main loop

fft_spectrum = 0;

for i=1:spectnum

start = (i-1)*offset;

sw = signal((1+start):(start+nfft),:).*(window*ones(1,channels));

fft_spectrum = fft_spectrum + (abs(fft(sw))*4/nfft); % X=fft(x.*hanning(N))*4/N; % hanning only

end

fft_spectrum = fft_spectrum/spectnum; % to do linear averaging scaling

% one sidded fft spectrum % Select first half

if rem(nfft,2) % nfft odd

select = (1:(nfft+1)/2)';

else

select = (1:nfft/2+1)';

end

fft_spectrum = fft_spectrum(select,:);

freq_vector = (select - 1)*Fs/nfft;

end

##### 2 Comments

Mathieu NOE
on 22 Jun 2022

hello

yes this code works fine (with both files). Now the x label is not Herz (1 / second) but the inverse of a length (1/ length)

maybe you feel more confortable with the x axis in log scale so the long "wavelength" defaults are easier to see (long wave = lower values of (1/ length)

% %% demo code for fft

% Fs = 8000/2; % Sampling frequency (number of samples/sampling time)

% T = 1/Fs; % Sampling period

% L = 8000; % Length of signal

% t = (0:L-1)*T; % Time vector

% X = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);

%% application to our case

filename = ('nominalprofile2.txt');

S = readmatrix(filename); % 2 channels : Time + signal

time = S(:,1);

dt = mean(diff(time));

X = S(:,2); % roughness signal

Fs = 1/dt; % sampling rate

[samples] = length(X);

% time vector

t = (0:samples-1)*dt;

figure;

plot(t,X)

title('Signal')

xlabel('length (mm)');

ylabel('Roughness value (um)')

Y = fft(X);

P2 = abs(Y/L);

P1 = P2(1:L/2+1);

P1(2:end-1) = 2*P1(2:end-1);

f = Fs*(0:(L/2))/L;

figure;

semilogx(f,P1)

title('Single-Sided Amplitude Spectrum of X(length)')

xlabel('1/length (1/mm)');

ylabel('Roughness value (um)')

### See Also

### Community Treasure Hunt

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

Start Hunting!