# FFT analysis with window of vibration during milling

3 views (last 30 days)

Show older comments

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

##### 0 Comments

### Answers (2)

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 = 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;

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

.

##### 8 Comments

Bruno Luong
on 10 Mar 2023

Edited: Bruno Luong
on 10 Mar 2023

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
on 10 Mar 2023

Calculating the frequencies —

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;

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)

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)

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)

.

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.

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])

##### 3 Comments

Paul
on 11 Mar 2023

Bruno Luong
on 11 Mar 2023

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!