FFT of bearing signal

Hi everyone, i have a problem : when i do the fft of a signal i don't see the first harmonic , instead the fft show the first harmonic multiplied for 3. In fact my fundamental harmonic is 16,66 Hz e when i do the FFT it show 49,80 Hz. Is possible this? I insert the code and the signal on this. Also the sampling frequency is 10kHz.

S=size(N);
L2=S(1,1);
NFFT= 2^nextpow2(L2);
Y=fft(N,NFFT)/L2;
f1=f/2*linspace(0,1,NFFT/2+1);
f2=2*abs(Y(1:NFFT/2+1));

[David Goodmanson]

Hi Eris,

An fft of this signal (without the unnecessary zerofilling done by NFFT) shows the first large peak at 50.048 Hz, and nothing much going on at or near 16.66 Hz. There are strong peaks at 250,350,550,650, which if 50 Hz were the fundamental, would be harmonics 5,7,11,13.

[Eris Prifti]

Hi David,

i have a shaft that rotate at 1000 rpm and i received that file from accelerometer and when i do the FFT show that spectrum. The problem is that, when i calculate the fundamentat like 1000/60= 16,66 ; but in the spectrum is 50Hz and it isn't dhe fundamental harmonic for me.

[Jon]

I'm not seeing anything in either your attached code or data file that indicates what the sampling frequency is. I don't know how you can scale the fft in Hz without knowing the sampling frequency.

Also in your attached code, you reference a variable called N, but the attached data file only has one column of data named x11. You also don't tell us the value of f in your code, is that the sampling frequency? If so what is it's value?

[Eris Prifti]

Sorry,

for fs is 10Khz and N indicate the size of the column named x11. I hope this is sufficiente

S=size(x11);

Unrecognized function or variable 'x11'.

L2=S(1,1);
NFFT= 2^nextpow2(L2);
Y=fft(x11,NFFT)/L2;
f1=fs/2*linspace(0,1,NFFT/2+1);
f2=2*abs(Y(1:NFFT/2+1));

Answer 1

[Jon]

In the attached code, I independently calculated the spectrum and compared to your calculation. It seems correct that your data has a very strong component at 50 Hz and higher harmonics of 50 Hz. Note also in the time domain plot, the 0.02 second period (50 Hz) oscillation is clearly visible. This gives added confidence that the scaling is correct. I also see a smaller peak at approx. 28 Hz. I agree that there is nothing significant at harmonics of the rotating frequency, which you report as 1000rev/min*1 min/60s = 16.667 Hz.

So either there is extremely little imbalance in the shaft, or perhaps the shaft frequency is not actually 1000 rpm. Regarding the 50Hz component, are you in a country where the electrical line frequency is 50 Hz. If so you are likely picking up "hum" from the power supply, or nearby power lines.

% % % S=size(N);
% % % L2=S(1,1);
% % % NFFT= 2^nextpow2(L2);
% % % Y=fft(N,NFFT)/L2;
% % % f1=f/2*linspace(0,1,NFFT/2+1);
% % % f2=2*abs(Y(1:NFFT/2+1));

% Assign parameters
fs = 10e3; % sampling frequency

% Load data from data file, data is in x11
load('vettore.mat')

% Find the length of the signal
L = numel(x11);

% Make sure that length of signal is even, if it isn't throw away last
% value to make it even
if rem(L,2)~=0
    x11 = x11(1:end-1);
    L = L-1;
end

% Compute the fft
Y = fft(x11);

% Compute the two sided spectrum P2
P2 = abs(Y/L);

% Compute the single sided spectrum P1
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);

% Define the frequency domain f and plot the single sided spectrum, P1
f = fs*(0:(L/2))/L;

% Plot results

% Plot the time domain signal
t = (0:L-1)*1/fs;
figure
plot(t,x11)
xlabel("time [s]")
ylabel("x11")
title("time domain signal x11")
xlim([0,0.1]); % zoom in

% Plot spectrum
figure
plot(f,P1)
title("Single-Sided Amplitude Spectrum of X(t)")
xlabel("f [Hz]")
ylabel("|P1(f)|")
xlim([0,300])

% Compare to OP's code
S=size(x11);
L2=S(1,1);
NFFT= 2^nextpow2(L2);
Y=fft(x11,NFFT)/L2;
f1=fs/2*linspace(0,1,NFFT/2+1);
f2=2*abs(Y(1:NFFT/2+1));
figure
plot(f1,f2)
title("OP's Spectrum")
xlabel("frequency [Hz]")
ylabel("Magnitude")
xlim([0,300])

[Jon]

Regarding your question/comment on how to remove the 50Hz line frequency "hum", please seehttps://www.mathworks.com/help/signal/ug/remove-the-60-hz-hum-from-a-signal.html, of course you will have to change the frequency in the example from 60 to 50, but this should be quite obvious.

[Eris Prifti]

But the link tell only how to remove 50 Hz, how can i remove the multiples of the harmonic like 100 Hz, 150Hz.

[Jon]

If you have the dsp toolbox you could use the function iircomb. Here's an example

% Assign parameters
fs = 10e3; % sampling frequency
bw = 0.005; % normalized bandwith of comb notch filter
lineFreq = 50; % line frequency in Hz

% Load data from data file, data is in x11
load('vettore.mat')

% Find the length of the signal
L = numel(x11);

% Make sure that length of signal is even, if it isn't throw away last
% value to make it even
if rem(L,2)~=0
    x11 = x11(1:end-1);
    L = L-1;
end

% Compute the fft
Y = fft(x11);

% Compute the two sided spectrum P2
P2 = abs(Y/L);

% Compute the single sided spectrum P1
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);

% Define the frequency domain f and plot the single sided spectrum, P1
f = fs*(0:(L/2))/L;

% Plot results

% Plot the time domain signal
t = (0:L-1)*1/fs;
figure
plot(t,x11)
xlim([0,0.1]); % zoom in

% Plot spectrum
figure
plot(f,P1)
title("Single-Sided Amplitude Spectrum of X(t)")
xlabel("f [Hz]")
ylabel("|P1(f)|")
xlim([0,300])

% Implement comb filter to filter 50Hz hum and harmonics
n = ceil((fs/2)/lineFreq)*2; % number of harmonics
[b,a] = iircomb(n,bw,'notch');


h = freqz(b,a,L/2+1,fs);
figure
plot(f,20*log10(abs(h)))
xlabel("f [Hz]")
ylabel("Magnitude [dB]")
title('Comb Filter Magnitude Response ')

% Filter the signal
x11f = filter(b,a,x11);

% FFT the filtered signal
% Compute the fft
Yf = fft(x11f);

% Compute the two sided spectrum P2
P2f = abs(Yf/L);

% Compute the single sided spectrum P1
P1f = P2f(1:L/2+1);
P1f(2:end-1) = 2*P1f(2:end-1);

% Define the frequency domain f and plot the single sided spectrum, P1
f = fs*(0:(L/2))/L;
% Plot spectrum
% % % figure
% % % plot(f,P1f)
% % % title("Single-Sided Amplitude Spectrum of notch filtered X(t)")
% % % xlabel("f [Hz]")
% % % ylabel("|P1(f)|");

% Plot it again, zoomed in
figure
plot(f,P1f)
title("Single-Sided Amplitude Spectrum of notch filtered X(t)")
xlabel("f [Hz]")
ylabel("|P1(f)|");

xlim([0,300])

[Jon]

If you are also running the experiments and/or have some involvement with the instrumentation and data acquisition, it would be far better to resolve your grounding/shielding issues to avoid the 50Hz hum and its harmonics in the first place.

[Eris Prifti]

Thank you Jon,

i don't know if my acquisition system is screened. Also i have anothwer question: is possbile to select a range of frequency like 90-93 and the harmonics of this range?

[Jon]

You can modify the line where the notch frequencies of the comb filter are calculated, for example,

 n = floor((fs/2)/91.5)*2
 [b,a] = iircomb(n,bw,'notch')

[Eris Prifti]

But if i apply this to the fft like the code below , it delete the information that i want to display?

% Filter the signal
x11f = filter(b,a,x11);
% FFT the filtered signal
% Compute the fft
Yf = fft(x11f);

[Jon]

I'm sorry, I don't understand what you are asking. Please explain what issue you are having now

[Eris Prifti]

the vector.mat is a measurement taken by an accelerometer, i don't know if the accelerometer is shielded. So I would like to remove all the disturbing frequencies like 50Hz of the current. The signal was taken from a measurement of a rolling bearing and I want to see only the failure frequency of the bearing which is 93.33Hz. Is this possible?

[Jon]

OK, I misunderstood your earlier question. So you only want to look at frequencies around 93.33 Hz and possible the harmonics of that.

The iircomb function offers a 'peak' option, which presumably just keeps frequencies at the selected values, rather than eliminates values at the selected frequency. You can read the documentation (type doc iircomb) and try that.

I think you may have practical problems though, since your frequency of interest, 93.33, Hz is so close to your noise, harmonic at 100 Hz. You would need a very narrow band to cut the 100 Hz and preserve the 93.33. You can try, but there may be limits to what you can do with digital signal processing when the frequency of your signals of interest are so close to your noise.

As I said earlier, you really need to determine why your accelerometer signals are so contaminated with 50 Hz line hum, need to look at shielding and grounding (beyond the scope of MATLAB answers).

[Jon]

Hi, if this answered your question, can you please accept the answer, so that others will know that an answer is available.

[Jon]

OK good luck with your project, if you have more questions on this topic you can always add to the thread.

[Eris Prifti]

Jon i have a question,

is possible to do a multiband pass filter to select the range of frequencies like 80-100 , 160-180?

[Jon]

Maybe these contributions from @Star Strider give you some ideas:

https://www.mathworks.com/matlabcentral/answers/704287-how-to-design-a-fir-filter-with-two-pass-bands?s_tid=prof_contriblnk

https://www.mathworks.com/matlabcentral/answers/503999-filtering-for-multiple-band-of-frequncies

[Eris Prifti]

Hi jon,

i have another question : is possible to impement a Homomorphic filtering demodulation for signal like bearing in matlab and can you write the code if possible pls?

[Jon]

Sorry, I'm not familiar with that technique. In general, I can offer you advice on particular problems you are having, but don't have time to write actual applications for people. If you start to work on your Homomorphic demodulation, and have specific MATLAB problems please open up a new question so others can see it and help you.

Answer 2

[Eris Prifti]

Hi david,

the 50 Hz is the hum frequency of my electrical line frequency . Is possible to filter this harmonic and its harmonics?

[Jon]

To keep the flow of this thread, it would be better to please put this as a comment to an anwer, rather than as an answer to your own question.