FFT analysis with window of vibration during milling

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

Answer 1

[Star Strider]

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 = readtable('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1319230/T1_1-p2.xlsx')

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

.

[Szymon Kurpiel]

Thanks a lot! Probably it works with order number=3, with other value there was created a lot of peaks. I made a quick checking - the bigest pick in f=13.27Hz. During experiment I had used a tool with 4 cutting edges and rotational speed 3183rpm. Based on this data, a peak is expected in f_exp=(3183)/(60*4)=13.2625.

[Star Strider]

My pleasure!

The difference between what you expect and what the fft call produced is likely due to the finite nature of the frequency resolution. You can increase the frequency resolution by:

NFFT = 2^(nextpow2(L)+k);

where ‘k’ is any positive integer. That may increase the frequency resolution enough to approximate the value you expect.

Also, the frequency vector will have to be adjusted to account for that. The approach I use is:

f = linspace(0, 1, NFFT/2+1)*(fs/2);

Iv = 1:numel(f);

and then use ‘Iv’ to index into the ‘P’ vector.

Using findpeaks can make this automatic:

[pks,locs] = findpeaks(P(Iv), 'MinPeakProminence',0.1);

fpeaks = f(locs)

That may need to be tweaked for each data column to give the best result for each ‘P’ vector.

.

[Paul]

@Szymon Kurpiel, if you zoom in on the plot, you'll see the peak is at 11.13 Hz, not 13.12 Hz. However,

@Star Strider, shouldn't the equation for fs be

%f = fs*(0:(N/2))/N;
f = fs*(0:(NFFT/2))/NFFT;

Similarly, shouldn't the Pi be indexed to fix(NFFT/2)

[Star Strider]

I wrote a different version of ‘f’ using linspace.

The value of ‘NFFT’ is always going to be even, by definition, so the result of dividing it by 2 will always be an integer.

[Paul]

AFAICT, the original answer does not define f using linspace. In the answer, the fft length is NFFT, which is not equal to N. So f should be determined from NFFT, not N.

[Star Strider]

My recoding of ‘f’ does just that.

[Bruno Luong]

I agree with Paul, it is wrong using N in one-sided FFT truncation and F

P1 = P(1:fix(N/2)+1);
P1(2:end-1) = 2*P1(2:end-1);
f = fs*(0:(N/2))/N;

That's why the frequncy peak found is wrong at 13.2 Hz instead of 8.5 Hz

[Star Strider]

Calculating the frequencies —

T = readtable('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1319230/T1_1-p2.xlsx')

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

.

Answer 2

[Bruno Luong]

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.

T = readtable('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1319230/T1_1-p2.xlsx');
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]

To illustrate that the filtering is free from FFT domain

T = readtable('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1319230/T1_1-p2.xlsx');
fs = 1000; % sampling frequency
L = size(T,1);
T1 = T;

NFFT = 2^nextpow2(L);
S1=T1.VarName1;
X=fft(hann(length(S1)).*S1, NFFT);
cutoff = 40; % or 20 [Hz] unit
ic = round(cutoff*NFFT/fs)+1;
X(ic+1:end-ic+1) = 0;
Sf1 = ifft(X,NFFT);

tfun=@(s) 1/fs*(0:length(s)-1);
plot(tfun(S1),S1,'g',tfun(Sf1),Sf1,'r')
xlim([40 42])
legend('original','fft filtering')

[Paul]

@Bruno Luong

Did you notice that the FFT has spikes at integer multiples of ~2.12-2.13 Hz all the way out to the Nyquist frequency? The "main" peak at low frequency is at 4*2.12 Hz, but there are other peaks of similar magnitude at higher frequency, like at the 140th multiple, in the y-axis.

[Bruno Luong]

Yes I notice the harmonics. The 2.12 Hz is probably mill frequency, and other frequency is something are forced to be resonate.