# FFT analysis with window of vibration during milling

9 views (last 30 days)
Szymon Kurpiel on 9 Mar 2023
Commented: Bruno Luong on 11 Mar 2023
Hello,
During milling experiment I had measured an acceleration on x, y, z axis. My problem is that after generating FFT graph there is a fair amount of noise. Even if I use a window for filtering the graph shape remains unchanged - as if the windows (hanning, hamming, blackman etc.) were wrong implemented. I have tried correcting the syntax in various ways.
I am attaching a data file with the oscillation waveform. In the columns are the accelerations on the axes - x, y, z respectively. I import the data into matlab in the form of matrices and perform the analysis on them. Could somebody help with finding some errors?
fs = 1000; % sampling frequency
N=99999; % Block size
subplot(311)
X=fft(hann(length(VarName1)).*VarName1);
P = abs(X/N);
P1 = P(1:N/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = fs*(0:(N/2))/N;
plot(f,P1, 'r')
xlabel("Częstotliwość (Hz)")
ylabel({'Przyspieszenie', 'oś x [mm/s^2]'})
title('Jednostronne widmo amplitudowe') , shg
grid on
grid minor
subplot(312)
Y=fft(hann(length(VarName2)).*VarName2);
P = abs(Y/N);
P2 = P(1:N/2+1);
P2(2:end-1) = 2*P2(2:end-1);
f = fs*(0:(N/2))/N;
plot(f,P2, 'g')
xlabel("Częstotliwość (Hz)")
ylabel({'Przyspieszenie', 'oś y [mm/s^2]'})
grid on
grid minor
subplot(313)
Z=fft(hann(length(VarName3)).*VarName3);
P = abs(Z/N);
P3 = P(1:N/2+1);
P3(2:end-1) = 2*P3(2:end-1);
f = fs*(0:(N/2))/N;
plot(f,P3, 'b')
xlabel("Częstotliwość (Hz)")
ylabel({'Przyspieszenie', 'oś z [mm/s^2]'})
grid on
grid minor

Star Strider on 9 Mar 2023
The signals have broadband noise, so any sort of frequency-selective filter will not adequately eliminate it. There are two ways to deal with broadband noise, one is with wavelet denoising, and the other (that I usually use) is the Savitzky-Golay filter, sgolayfilt. You have the Signal Processing Toolbox (since you use the hann function) so you also have sgollayfilt. I usually use a 3-order filter, and then adjust the frame length (that has to be an odd number here) until I get the result I want. I cannot determine what that is, so I defer to you to experiment with it. (Analysing the sgolayfilt function output after the filtering operation reveals it to be a sort of comb filter, however unlike FIR filters, it is highly adaptable, and signal length is not a problem for it.)
A lowpass filter with a cutoff of about 25 Hz could also work, however it would not eliminate the noise in the signal below that frequency.
T = 99999×3 table
VarName1 VarName2 VarName3 _________ ________ ________ 1.3257 8.6069 4.8746 0.21254 1.9944 -17.529 -0.034467 1.6997 0.94046 3.151 2.2187 10.767 1.6229 3.1303 -8.7216 2.3658 -1.7097 1.142 0.056947 1.5759 8.7671 2.515 0.054625 -4.9641 2.3706 4.8615 -0.94778 -0.36672 -6.9198 -6.2516 1.5539 -0.11157 0.45935 1.6676 0.6172 8.0744 1.0386 -4.0158 3.6676 -1.3781 -4.2896 -0.51666 1.9165 2.3299 1.2179 0.069988 -3.2052 -5.1812
fs = 1000; % sampling frequency
L = size(T,1);
N = L;
T1 = T;
order = 3;
framelen = 101;
for k = 1:size(T,2)
T1{:,k} = sgolayfilt(T{:,k}, order, framelen);
end
t = linspace(0, L-1, L)/fs;
figure
subplot(2,1,1)
plot(t, T{:,:})
xlim([0 1])
ylim([-1 2])
grid
title('Original Signal')
subplot(2,1,2)
plot(t, T1{:,:})
xlim([0 1])
grid
title('Filtered Signal')
NFFT = 2^nextpow2(L);
figure
subplot(311)
X=fft(hann(length(T1.VarName1)).*T1.VarName1, NFFT);
P = abs(X/N);
P1 = P(1:fix(N/2)+1);
P1(2:end-1) = 2*P1(2:end-1);
f = fs*(0:(N/2))/N;
plot(f,P1, 'r')
xlabel("Częstotliwość (Hz)")
ylabel({'Przyspieszenie', 'oś x [mm/s^2]'})
title('Jednostronne widmo amplitudowe') , shg
grid on
grid minor
subplot(312)
Y=fft(hann(length(T1.VarName2)).*T1.VarName2, NFFT);
P = abs(Y/N);
P2 = P(1:fix(N/2)+1);
P2(2:end-1) = 2*P2(2:end-1);
% f = fs*(0:(N/2))/N;
plot(f,P2, 'g')
xlabel("Częstotliwość (Hz)")
ylabel({'Przyspieszenie', 'oś y [mm/s^2]'})
grid on
grid minor
subplot(313)
Z=fft(hann(length(T1.VarName3)).*T1.VarName3, NFFT);
P = abs(Z/N);
P3 = P(1:fix(N/2)+1);
P3(2:end-1) = 2*P3(2:end-1);
% f = fs*(0:(N/2))/N;
plot(f,P3, 'b')
xlabel("Częstotliwość (Hz)")
ylabel({'Przyspieszenie', 'oś z [mm/s^2]'})
grid on
grid minor
.
Star Strider on 10 Mar 2023
Calculating the frequencies —
T = 99999×3 table
VarName1 VarName2 VarName3 _________ ________ ________ 1.3257 8.6069 4.8746 0.21254 1.9944 -17.529 -0.034467 1.6997 0.94046 3.151 2.2187 10.767 1.6229 3.1303 -8.7216 2.3658 -1.7097 1.142 0.056947 1.5759 8.7671 2.515 0.054625 -4.9641 2.3706 4.8615 -0.94778 -0.36672 -6.9198 -6.2516 1.5539 -0.11157 0.45935 1.6676 0.6172 8.0744 1.0386 -4.0158 3.6676 -1.3781 -4.2896 -0.51666 1.9165 2.3299 1.2179 0.069988 -3.2052 -5.1812
fs = 1000; % sampling frequency
L = size(T,1);
N = L;
T1 = T;
order = 3;
framelen = 101;
for k = 1:size(T,2)
T1{:,k} = sgolayfilt(T{:,k}, order, framelen);
end
t = linspace(0, L-1, L)/fs;
figure
subplot(2,1,1)
plot(t, T{:,:})
xlim([0 1])
ylim([-1 2])
grid
title('Original Signal')
subplot(2,1,2)
plot(t, T1{:,:})
xlim([0 1])
grid
title('Filtered Signal')
NFFT = 2^nextpow2(L);
f = linspace(0, 1, NFFT/2+1)*(fs/2);
Iv = 1:numel(f);
figure
subplot(311)
X=fft(hann(length(T1.VarName1)).*T1.VarName1, NFFT);
P = abs(X/N);
% P1 = P(1:fix(NFFT/2)+1);
P1 = P(Iv);
% P1(2:end-1) = 2*P1(2:end-1);
% f = fs*(0:(N/2))/N;
[pks,locs] = findpeaks(P1, 'MinPeakProminence',0.1);
fpeaks = f(locs)
fpeaks = 8.4915
plot(f,P1, 'r')
xlabel("Częstotliwość (Hz)")
ylabel({'Przyspieszenie', 'oś x [mm/s^2]'})
title('Jednostronne widmo amplitudowe') , shg
grid on
grid minor
subplot(312)
Y=fft(hann(length(T1.VarName2)).*T1.VarName2, NFFT);
P = abs(Y/N);
P2 = P(1:fix(N/2)+1);
% P2(2:end-1) = 2*P2(2:end-1);
P2 = P(Iv);
% f = fs*(0:(N/2))/N;
plot(f,P2, 'g')
xlabel("Częstotliwość (Hz)")
ylabel({'Przyspieszenie', 'oś y [mm/s^2]'})
grid on
grid minor
[pks,locs] = findpeaks(P(Iv), 'MinPeakProminence',0.1);
fpeaks = f(locs)
fpeaks = 8.4839
subplot(313)
Z=fft(hann(length(T1.VarName3)).*T1.VarName3, NFFT);
P = abs(Z/N);
P3 = P(1:fix(N/2)+1);
% P3(2:end-1) = 2*P3(2:end-1);
% f = fs*(0:(N/2))/N;
P3 = P(Iv);
plot(f,P3, 'b')
xlabel("Częstotliwość (Hz)")
ylabel({'Przyspieszenie', 'oś z [mm/s^2]'})
grid on
grid minor
[pks,locs] = findpeaks(P(Iv), 'MinPeakProminence',0.05);
fpeaks = f(locs)
fpeaks = 8.4839
.

Bruno Luong on 10 Mar 2023
Edited: Bruno Luong on 10 Mar 2023
I don't see the need of calling sgolayfilt. It is just a linear filter so it is equivalent to scale down the signal in frequency domain by a filter frequency response.
As OP is interested in signal about 13 Hz, just put empirical amplitude above 20 Hz to 0, it does reveal the amplitude peaks in the low frequency bandwidth as well without the complexity of calling sgolayfilt.
And no the main pic is not 13 Hz but about 8.5 Hz. If you look at @Star Strider filtered signal during 1s, there are more than 8 periods in 1second, not 13.
fs = 1000; % sampling frequency
L = size(T,1);
N = L;
T1 = T;
t = linspace(0, L-1, L)/fs;
NFFT = 16*2^nextpow2(L);
figure
subplot(311)
X=fft(hann(length(T1.VarName1)).*T1.VarName1, NFFT);
P = abs(X/N);
P1 = P(1:fix(NFFT/2)+1);
P1(2:end-1) = 2*P1(2:end-1);
f = fs*(0:(NFFT/2))/NFFT;
P1(f>20)=0; % filter in frequency domain
plot(f,P1, 'r')
xlabel("Częstotliwość (Hz)")
ylabel({'Przyspieszenie', 'oś x [mm/s^2]'})
title('Jednostronne widmo amplitudowe') , shg
grid on
grid minor
subplot(312)
Y=fft(hann(length(T1.VarName2)).*T1.VarName2, NFFT);
P = abs(Y/N);
P2 = P(1:fix(NFFT/2)+1);
P2(2:end-1) = 2*P2(2:end-1);
P2(f>20)=0; % filter in frequency domain
% f = fs*(0:(N/2))/N;
plot(f,P2, 'g')
xlabel("Częstotliwość (Hz)")
ylabel({'Przyspieszenie', 'oś y [mm/s^2]'})
grid on
grid minor
subplot(313)
Z=fft(hann(length(T1.VarName3)).*T1.VarName3, NFFT);
P = abs(Z/N);
P3 = P(1:fix(NFFT/2)+1);
P3(2:end-1) = 2*P3(2:end-1);
P3(f>20)=0; % filter in frequency domain
% f = fs*(0:(N/2))/N;
plot(f,P3, 'b')
xlabel("Częstotliwość (Hz)")
ylabel({'Przyspieszenie', 'oś z [mm/s^2]'})
grid on
grid minor
copyobj(get(gcf,'children'), figure)
set(get(gcf,'children'),'xlim',[7 10])
Bruno Luong on 11 Mar 2023
Yes I notice the harmonics. The 2.12 Hz is probably mill frequency, and other frequency is something are forced to be resonate.