Surface roughness analysis from raw data
79 views (last 30 days)
Show older comments
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.
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!