Surface roughness analysis from raw data

115 views (last 30 days)
Hello! I need some help with this: I have roughness values of a surface, and now I want to filter the values in order to differentiate between roughness and waviness. I read I should use FFT and short and long wave filtering but I am a bit lost with it. Find attached the data (first column is length (milimeter) and second roughness value (micrometer)). I guess the first column could be considered as time... Any help is welcome! Thanks in advance.

Answers (1)

Mathieu NOE
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
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)')

Sign in to comment.

Categories

Find more on Vibration Analysis in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!