My oscilloscope peaks do not match my FFT peaks.

Hi all, I am trying to perform a frequency analysis on an electric motor in relation to the normalized frequencies relative to the rotational frequency.

Thanks to previous measurements with the oscilloscope I obtained frequencies in certain peaks as in picture.

Now I am trying to obtain the same with our data acquisition device and an fft code, but the data does not match (picture)

In the description I leave my code with the relevant information.

clear all
close all
% % % For WorkStation dSpace recorders!
file_name = 'exp1_HV_eaux_B0_ASQ_N7_vibrations_NT3a_pos_3000rpm_20kHz_foam';
s=load(file_name);
%%% Select window type for FFT
% a=no window, b=Rectangular, c=Hann, d=Hamming, e=Flattop, f=Blackman-Harris, g=Nuttall, h=Chebyshev
winType = 'b';  
length=length(s.exp1_HV_eaux_B0_ASQ_N7_vibrations_NT3a_pos_3000rpm_20kHz_foam.X(1).Data);
x=s.exp1_HV_eaux_B0_ASQ_N7_vibrations_NT3a_pos_3000rpm_20kHz_foam.X(1).Data(1:length);
%%% Important base recording params
Fs=20000;                 %%% sampling frequency [Hz]
T=1/Fs;
%%% important base FFT params
fft_cycles          = 2            ;
speed_aver_window   = 1; %%% in sec
speed_aver_wind_points = speed_aver_window / T;
multipleORHz        = 1;  %%% 1 - multiply, 2 is Hz
%%% main indexis
index_speed_kistler     = 1;
index_torque_an         = 2;
index_torque_an_filt    = 3;
index_vibr              = 4;
index_vibr_filt         = 5;
index_iq_act            = 7;
index_id_act            = 6;
index_speed_sw          = 8;
index_speed_sew         = 9;
%%% Define staret point of FFT:
time_start_fft  = 22;
point_start_fft = round(time_start_fft / T);
%%%% Indexes for analysis
index_main_speed  = index_speed_sw;
index_main_torque = index_torque_an_filt;
index_main_vibr   = index_vibr;
speed_main  = s.exp1_HV_eaux_B0_ASQ_N7_vibrations_NT3a_pos_3000rpm_20kHz_foam.Y(index_main_speed).Data(1:length);
torque_main = s.exp1_HV_eaux_B0_ASQ_N7_vibrations_NT3a_pos_3000rpm_20kHz_foam.Y(index_main_torque).Data(1:length);
vibr_main   = s.exp1_HV_eaux_B0_ASQ_N7_vibrations_NT3a_pos_3000rpm_20kHz_foam.Y(index_main_vibr).Data(1:length);
%%% Extraction of currents
Id_act_1=s.exp1_HV_eaux_B0_ASQ_N7_vibrations_NT3a_pos_3000rpm_20kHz_foam.Y(index_id_act).Data(1:length);
Iq_act_1=s.exp1_HV_eaux_B0_ASQ_N7_vibrations_NT3a_pos_3000rpm_20kHz_foam.Y(index_iq_act).Data(1:length);
%%% Extraction of speed and torque
SEW_speed=s.exp1_HV_eaux_B0_ASQ_N7_vibrations_NT3a_pos_3000rpm_20kHz_foam.Y(index_speed_sew).Data(1:length);
Kistler_speed=s.exp1_HV_eaux_B0_ASQ_N7_vibrations_NT3a_pos_3000rpm_20kHz_foam.Y(index_speed_kistler).Data(1:length);
Kistler_torque_an=s.exp1_HV_eaux_B0_ASQ_N7_vibrations_NT3a_pos_3000rpm_20kHz_foam.Y(index_torque_an).Data(1:length);
Kistler_torque_an_filt=s.exp1_HV_eaux_B0_ASQ_N7_vibrations_NT3a_pos_3000rpm_20kHz_foam.Y(index_torque_an_filt).Data(1:length);
N_pres=s.exp1_HV_eaux_B0_ASQ_N7_vibrations_NT3a_pos_3000rpm_20kHz_foam.Y(index_speed_sw).Data(1:length);
vibr=s.exp1_HV_eaux_B0_ASQ_N7_vibrations_NT3a_pos_3000rpm_20kHz_foam.Y(index_vibr).Data(1:length);
vibr_filt=s.exp1_HV_eaux_B0_ASQ_N7_vibrations_NT3a_pos_3000rpm_20kHz_foam.Y(index_vibr_filt).Data(1:length);
%%% plot results of the test
set(gcf,'color','white')
ax1=subplot(5,1,1);
plot(x,Kistler_torque_an,'r', x,Kistler_torque_an_filt,'b');
title('Torque Kistler');
xlabel('Time [s]')
ylabel('Torque [Nm]');
grid on
ax2=subplot(5,1,2);
plot(x,N_pres,'r',x,SEW_speed,'b',x,Kistler_speed,'g');
title('Speed');
xlabel('Time [s]')
ylabel('Speed [rpm]');
legend('Sensor','SEW','Kistler');
grid on
ax3=subplot(5,1,3);
plot(x,vibr,'.- r', x,vibr_filt,'b');
title('Vibrations');
xlabel('Time [s]')
ylabel('Acceleration [g]');
grid on
ax4=subplot(5,1,4);
plot(x,Id_act_1,'r');
title('Id act vs ref');
xlabel('Time [s]')
ylabel('Id [A]');
grid on
ax5=subplot(5,1,5);
plot(x,Iq_act_1,'r');
title('Iq act vs ref');
xlabel('Time [s]')
ylabel('Iq [A]');
grid on
linkaxes([ax1,ax2,ax3,ax4,ax5],'x');
%%% Making window to analyze FFT
speed_main_rpm_aver     = abs(mean(speed_main(point_start_fft:(point_start_fft + speed_aver_wind_points))));
speed_rps               = abs(speed_main_rpm_aver / 60);
period                  = 1 / speed_rps;
window_sec              = fft_cycles * period;
window_poins            = round(window_sec * Fs);
%%% Select and build desired window
switch lower(winType)
    case 'a'   % No window (raw data)
        win = ones(window_poins,1);
    case 'b'   % Rectangular
        win = rectwin(window_poins);
    case 'c'   % Hann
        win = hann(window_poins);
    case 'd'   % Hamming
        win = hamming(window_poins);
    case 'e'   % Flattop
        win = flattopwin(window_poins);
    case 'f'   % Blackman-Harris
        win = blackmanharris(window_poins);
    case 'g'   % Nuttall
        win = nuttallwin(window_poins);
    case 'h'   % Chebyshev
        win = chebwin(window_poins, 60); % 60 dB sidelobe suppression
    otherwise
        error('Unknown winType "%s"', winType);
end
% Compute mean values for DC removal
torque_main_aver = mean(torque_main(point_start_fft:(point_start_fft + speed_aver_wind_points)));
vibr_main_aver   = mean(vibr_main(point_start_fft:(point_start_fft + speed_aver_wind_points)));
% Preallocate arrays
time_wind    = zeros(window_poins, 1);
torque_wind  = zeros(window_poins, 1);
vibr_wind    = zeros(window_poins, 1);
speed_wind   = zeros(window_poins, 1);
% Create windowed signals
j = 1;
for i = point_start_fft:(point_start_fft + window_poins - 1)
    time_wind(j)    = x(i);
    torque_wind(j)  = (torque_main(i) - torque_main_aver) * win(j);
    vibr_wind(j)    = (vibr_main(i) - vibr_main_aver) * win(j);
    speed_wind(j)   = speed_main(i);
    j = j + 1;
end
figure
set(gcf,'color','white')
bx1=subplot(3,1,1);
plot(time_wind,torque_wind,'r');
title('Torque window');
xlabel('Time [s]')
ylabel('Torque [Nm]');
grid on
bx2=subplot(3,1,2);
plot(time_wind,vibr_wind,'r');
title('Acceleretion window');
xlabel('Time [s]')
ylabel('Acceleration [g]');
grid on
bx3=subplot(3,1,3);
plot(time_wind,speed_wind,'r');
title('Speed window');
xlabel('Time [s]')
ylabel('Speed [rpm]');
grid on
linkaxes([bx1,bx2,bx3],'x');
%%% FFT of Vibrations     
L=window_poins;
t = (0:L-1)*T;        % Time vector
Y1=fft(vibr_wind,L);
P2 = abs(Y1/L).^2;
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
if multipleORHz == 1
    freq_ref = speed_main_rpm_aver/60;
    freq_plot_lim = 100;
elseif multipleORHz == 2
    freq_ref = 1;
    freq_plot_lim = 4000;
end
f1 = (Fs*(0:(L/2))/L)/(freq_ref);
figure 
set(gcf,'color','white')
bar(f1,P1) 
title('Single-Sided Amplitude Spectrum of X(t)')
xlabel('f (multiple of mechanical frequency)')
xlim([0 freq_plot_lim])
ylabel('Absolute value of Harmonic VIBRATIONS [g]')

Here also some info about the signal

I would appreciate if in case you see a mistake in my code please let me know, I have no more ideas.

I attached a link where you can download the datafile here: exp1_HV_eaux_B0_ASQ_N7_vibrations_NT3a_pos_3000rpm_20kHz_foam.mat

Thank you in advance

[Mathieu NOE]

hello

it would help if you could also share the data file (exp1_HV_eaux_B0_ASQ_N7_vibrations_NT3a_pos_3000rpm_20kHz_foam)

[Eduardo]

Hi. Thank you for the reply and advice, sadly the data file is to big for the platform to handle, I transformed a description of the data array the data contains in the pdf attached.

I hope this helps.

If you need further things let me know

[Mathieu NOE]

that doesn't help me alot

maybe you can share a link to your big file

[Eduardo]

can you try now? the link should lead you to the file

[Mathieu NOE]

hello

this is a first attempt to double check your code - I built quickly a simplified version just to figure out where could the problem lie

Nb that the two speed data fields does not contain any value (only 0) (see lower subplot of figure 1) so I cannot normalize the freq data to the rpm values (I assumed 3000 RPM as it appears in your post)

• top subplot : raw vibration in [g] as assumed from your code
• bottom subplot : SEW_speed,Kistler_speed (zero !!)

at first glance my fft plot is quite similar to yours in terms of frequencies and amplitudes , but I am not surprised because if the time signal peak amplitude is about +/- 10 g then the sum of the fft peaks are also close to 10 (g)

on the other hand , this fft plot (below) shows peaks which summed amplitudes are way below the +/-10 g amplitude from the time plot (in red in your post)

so there is a few questions :

• where does this plot comes from / data / code ?
• in both cases , did you take into account the sensor sensivity, channel gains , ADC gains , so that a 1 g sine wave analysed from both systems show the same fft outputs ? that does not appear in your code but maybe it's done somewhere else (in your dSpace file ?)

FYI, below is my test code

NB the trick to avoid using very long filenames repeatdly in the code, inspired from How do I rename a structure? - MATLAB Answers - MATLAB Central

Code :

%% For WorkStation dSpace recorders!
file_name = 'exp1_HV_eaux_B0_ASQ_N7_vibrations_NT3a_pos_3000rpm_20kHz_foam';
s=load(file_name);

%% use this trick to avoid using the long filename every time 
StrName=fieldnames(s);
StrName=StrName{1};

%% main Code
% time = s.exp1_HV_eaux_B0_ASQ_N7_vibrations_NT3a_pos_3000rpm_20kHz_foam.X.Data;
time   = s.(StrName).X.Data;

dt = mean(diff(time));
Fs = 1/dt;

%% main indexes
index_speed_kistler     = 1;
index_torque_an         = 2;
index_torque_an_filt    = 3;
index_vibr              = 4;
index_vibr_filt         = 5;
index_iq_act            = 7;
index_id_act            = 6;
index_speed_sw          = 8;
index_speed_sew         = 9;
%%%% Indexes for analysis
index_main_speed  = index_speed_sw;
index_main_torque = index_torque_an_filt;
index_main_vibr   = index_vibr;

% vibr_main   = s.exp1_HV_eaux_B0_ASQ_N7_vibrations_NT3a_pos_3000rpm_20kHz_foam.Y(index_main_vibr).Data;
vibr_main   = s.(StrName).Y(index_main_vibr).Data;

%% Extraction of speed and torque
% SEW_speed=s.exp1_HV_eaux_B0_ASQ_N7_vibrations_NT3a_pos_3000rpm_20kHz_foam.Y(index_speed_sew).Data;
SEW_speed=s.(StrName).Y(index_speed_sew).Data;
% Kistler_speed=s.exp1_HV_eaux_B0_ASQ_N7_vibrations_NT3a_pos_3000rpm_20kHz_foam.Y(index_speed_kistler).Data;
Kistler_speed=s.(StrName).Y(index_speed_kistler).Data;

% extract time portion between t = 22 and 25 s
time_start_fft  = 22;
duration = 3;
ind = (time>=time_start_fft) & (time<=time_start_fft+duration);
time = time(ind);
vibr_main = vibr_main(ind);
SEW_speed = SEW_speed(ind);
Kistler_speed = Kistler_speed(ind);

speed_RPM_estimated = 3000;  % /!\ to be replaced with actual data (when available)

figure(1)
subplot(2,1,1),plot(time,vibr_main);
subplot(2,1,2),plot(time,SEW_speed,time,Kistler_speed); % NO SPEED DATA !!

%% FFT plot 
nfft = 1024*4;    % 
Overlap = 0.75; % (75% ) overlap 
df = Fs/nfft;

% vibr_main = sin(2*pi*10*df*time); % a simple test with a unitary amplitude sine wave

[freq,fft_spectrum] = myfft_peak(vibr_main, Fs, nfft, Overlap);

figure(2)
freq_normalized = freq*60/speed_RPM_estimated;
plot(freq_normalized,fft_spectrum)
title('FFT Spectrum')
ylabel('Amplitude')
xlabel('Frequency [normalized]')
xlim([0 25]);

figure(3) % check signal stationnarity with spectrogram : OK
nfft = 1024;    % 
Overlap = 0.75; % (75% ) overlap 
spectrogram(vibr_main,hanning(nfft),round(Overlap*nfft),nfft,Fs,'yaxis','minthreshold',-80)  % WINDOW,NOVERLAP,NFFT,Fs
colormap('jet')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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);

% transpose signal dimensions if channels >> samples
if channels > samples
    signal = signal';
    [samples,channels] = size(signal);
end

% 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

[Mathieu NOE]

I replaced the index index_speed_sw instaed of index_speed_sew

and then , bingo , I got some speed data ,the averaged speed is then

speed_RPM_estimated = 3.000577608046834e+03

which is very close to the supposed RPM (3000) , but basically that does not change anything to the plots and my comments above

updated code :

%% For WorkStation dSpace recorders!
file_name = 'exp1_HV_eaux_B0_ASQ_N7_vibrations_NT3a_pos_3000rpm_20kHz_foam';
s=load(file_name);

%% use this trick to avoid using the long filename every time 
StrName=fieldnames(s);
StrName=StrName{1};

%% main Code
% time = s.exp1_HV_eaux_B0_ASQ_N7_vibrations_NT3a_pos_3000rpm_20kHz_foam.X.Data;
time   = s.(StrName).X.Data;

dt = mean(diff(time));
Fs = 1/dt;

%% main indexes
index_speed_kistler     = 1;
index_torque_an         = 2;
index_torque_an_filt    = 3;
index_vibr              = 4;
index_vibr_filt         = 5;
index_iq_act            = 7;
index_id_act            = 6;
index_speed_sw          = 8;
index_speed_sew         = 9;
%%%% Indexes for analysis
index_main_speed  = index_speed_sw;
index_main_torque = index_torque_an_filt;
index_main_vibr   = index_vibr;

% vibr_main   = s.exp1_HV_eaux_B0_ASQ_N7_vibrations_NT3a_pos_3000rpm_20kHz_foam.Y(index_main_vibr).Data;
vibr_main   = s.(StrName).Y(index_main_vibr).Data;

%% Extraction of speed and torque
% SEW_speed=s.exp1_HV_eaux_B0_ASQ_N7_vibrations_NT3a_pos_3000rpm_20kHz_foam.Y(index_speed_sew).Data;
SEW_speed=s.(StrName).Y(index_speed_sw).Data;
% Kistler_speed=s.exp1_HV_eaux_B0_ASQ_N7_vibrations_NT3a_pos_3000rpm_20kHz_foam.Y(index_speed_kistler).Data;
Kistler_speed=s.(StrName).Y(index_speed_kistler).Data;

% extract time portion between t = 22 and 25 s
time_start_fft  = 22;
duration = 3;
ind = (time>=time_start_fft) & (time<=time_start_fft+duration);
time = time(ind);
vibr_main = vibr_main(ind);
SEW_speed = SEW_speed(ind);
Kistler_speed = Kistler_speed(ind);

speed_RPM_estimated = mean(SEW_speed);  

figure(1)
subplot(2,1,1),plot(time,vibr_main);
subplot(2,1,2),plot(time,SEW_speed,time,Kistler_speed);

%% FFT plot 
nfft = 1024*4;    % 
Overlap = 0.75; % (75% ) overlap 
df = Fs/nfft;

% vibr_main = sin(2*pi*10*df*time); % a simple test with a unitary amplitude sine wave

[freq,fft_spectrum] = myfft_peak(vibr_main, Fs, nfft, Overlap);

figure(2)
freq_normalized = freq*60/speed_RPM_estimated;
plot(freq_normalized,fft_spectrum)
title('FFT Spectrum')
ylabel('Amplitude')
xlabel('Frequency [normalized]')
xlim([0 25]);

figure(3) % check signal stationnarity with spectrogram : OK
nfft = 1024;    % 
Overlap = 0.75; % (75% ) overlap 
spectrogram(vibr_main,hanning(nfft),round(Overlap*nfft),nfft,Fs,'yaxis','minthreshold',-80)  % WINDOW,NOVERLAP,NFFT,Fs
colormap('jet')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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);

% transpose signal dimensions if channels >> samples
if channels > samples
    signal = signal';
    [samples,channels] = size(signal);
end

% 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

[Eduardo]

Hello Mathieu.

Thank you so much for your reply.

I will try to give as much detail as possible.

Im currently writting a paper whereI investigate the vibration behavior of an PMSM when mounted in a testbench.

Which brings us to point number 1: the reason why i have 3 values of speed is because i have one speed coming from a loadmachine that is connected to my system applying speed, then we have the speed of my motor thats being tested. and i have a software convertion for the speed. In this specific test i Isolated my PMSM and its the only thing active in the system, that is why only 1 out of 3 speeds has value non zero.

Now the second point. you mention that the sume of each individual peak has to sume to the value of the raw signal. is there a reason why is this so where i can verify it?

The 2nd fft plot that im comparing to is one made a year ago in a previous study using the exact same motor and testbench, but a 3 axial sensor that we borrowed instead of our own 1 axial one.

Regarding your question about the sensitivity of the sensors i dont have right now the info about the one in the previous study but i can take a look into it.

I also thank you for trying your own code but since you get the same data then im starting to question if the kast study made was done correctly. I will double check them but im glad that at least im not making a syntax error in my code

I appreciete the help and I will post as soon as i get better news!!

[Eduardo]

could you maybe share an image of the fft you are getting?

[Mathieu NOE]

hello again Eduardo

here is my fft plot (code is the same as in my last answer above)

according to what I know from PMSM motors vibration , we should see a dominant frequency which is the pole passing frequency (with some harmonics but with much lower amplitude)

so all those peaks that we see in this fft plot (above) are not a pole passing frequency + harmonics - this must come from either your load or degraded bearings or whatever....NB I don't know how many poles you have in your motor

Now the second point. you mention that the sume of each individual peak has to sume to the value of the raw signal. is there a reason why is this so where i can verify it?

the only and true relationship is the Parseval theorem (Parseval's theorem - Wikipedia) that states that the signal power (or rms, power = rms squared) is the same if you compute the power (or rms) in the time domain or in the frequency domain

so here , take the dominant peaks from your fft graph, compute the rms sum of them and you should get a similar number with the rms of the time signal

this is a good way to make a quick coherence check between the time domain plot and the fft plot

obviously here, with a signal that has a range of +/- 10g in the time domain (rms estimated around 5 / 6 g) , the very first fft plot in your post does not have peaks large enough so that their rms sum matches (the rms sum of the fft graph is only 1.3 g ( sqrt(0.15^2+0.92^2+0.42^2+0.22^2+0.35^2 + 0.73^2) )

NB : if you have a lot of broadband noise and peaks that do not emerge a lot from the noise, this comparison is less obvious to do, as picking only the fft peaks does not represent the entire power of the signal (underestimation))

on the other side the frequency content seems more like a motor pole passing frequency (p = 10 ?) + some harmonics

[Eduardo]

My motor has a structur 10p12s

[Eduardo]

do you have any idea why could be that there are big sidebands at 10 (i assume 8 and 11 if its same as my results) and then around 5 also

[Mathieu NOE]

I looked again at the fft plot

normally there should be no vibration at order 5 and its vicinity

there should be one peak exactly at order 10 and no sidebands

what puzzles me is that the sidebands are not occuring at any integer value (orders) but there is still a logic about what distance they have vs the supposed frequency like :

• "order 5" : sidebands at 4.3 and 5.7 (delta = 0.7 in both directions)
• "order 10" : sidebands at 8.6 and 11.4 (delta = 1.4 in both directions)

do you have a good analog antialiasing filter before data enter the ADC ? I ask this because I doubt that dSpace systems have one by default

[Eduardo]

Hi.

We do have a filter as far as I know. there is a responsable guy that does that but he comes back on monday.

In any case what is it particulary interesting about that filter exactly? for me to able to build some knowledge

[Mathieu NOE]

anti aliasing filters are analog low pass filters that should be used to avoid signals above the Nyquist frequency being replicated in the base band , there is ton of litterature on the internet on this subject

Anti-Aliasing Filters and Their Usage Explained - NI

Anti-aliasing filter - Wikipedia

adcs-aliasing-and-anti-aliasing-filters-presentation-quiz.pdf

of course , if you are 100% sure that there is no signal with significant energy in the very high frequencies , there is probably no issue with aliasing

there is maybe another problem with the setup so I was about to suggest that next time you run the test, you should record / analyse the signals with two separate acquisition systems to compare the results

as usual when in doubt there are the basic checks to do :

• sensor , cable ok ? (not damaged)
• sensor properly secured to the structure (not loosy ?)
• signal range ok ? no overload (too high gain before adc) ?
• acquisition system calibrated (with known test signals ) ?

[Eduardo]

Hi,

I see what you mean.

As far as i have read about this research there shouldnt be frequencies above the nyquist frequecy expected.

On top of that im using a sampling frequency of 40kHz which is already a bit less than double my max expected frequency

As far as another recorder as i mention in the titel of the post. I did an ozciloscope meassurment which showed me the same signal as the one I am trying to get

About the sensor setup i made sure cable is magnetic isolated and also there is a string that prevents the cable to pull down the sensor

[Mathieu NOE]

hello again

just for my own fun I analysed three portions of the record , that correspond to 600 / 1200 & 3000 RPM

% extract time portion 
time_start_fft  = 22;  % 4 / 11 / 22  (RPM = 600 / 1200 / 3000)
duration = 3;

what is really puzzling here is that there are always two peaks symetricaly positionned vs the theoretical dominant order (10) , and this distance decreases as RPM increases ... there must be an explanation somewhere...

600 RPM

1200 RPM

3000 RPM

[Eduardo]

Hi Mathieu,

Thanks for the observation.

I wanted to try with the codes you sent above but the spectrum for some reason is not matching, are you using a different script?

i get this for example in the 1200 rpm

[Mathieu NOE]

hello again

here's the code - condition 1200 RPM :

%% For WorkStation dSpace recorders!
file_name = 'exp1_HV_eaux_B0_ASQ_N7_vibrations_NT3a_pos_3000rpm_20kHz_foam';
s=load(file_name);

%% use this trick to avoid using the long filename every time 
StrName=fieldnames(s);
StrName=StrName{1};

%% main Code
% time = s.exp1_HV_eaux_B0_ASQ_N7_vibrations_NT3a_pos_3000rpm_20kHz_foam.X.Data;
time   = s.(StrName).X.Data;

dt = mean(diff(time));
Fs = 1/dt;

%% main indexes
index_speed_kistler     = 1;
index_torque_an         = 2;
index_torque_an_filt    = 3;
index_vibr              = 4;
index_vibr_filt         = 5;
index_iq_act            = 7;
index_id_act            = 6;
index_speed_sw          = 8;
index_speed_sew         = 9;
%%%% Indexes for analysis
index_main_speed  = index_speed_sw;
index_main_torque = index_torque_an_filt;
index_main_vibr   = index_vibr;

% vibr_main   = s.exp1_HV_eaux_B0_ASQ_N7_vibrations_NT3a_pos_3000rpm_20kHz_foam.Y(index_main_vibr).Data;
vibr_main   = s.(StrName).Y(index_main_vibr).Data;

%% Extraction of speed and torque
% SEW_speed=s.exp1_HV_eaux_B0_ASQ_N7_vibrations_NT3a_pos_3000rpm_20kHz_foam.Y(index_speed_sew).Data;
SEW_speed=s.(StrName).Y(index_speed_sw).Data;
% Kistler_speed=s.exp1_HV_eaux_B0_ASQ_N7_vibrations_NT3a_pos_3000rpm_20kHz_foam.Y(index_speed_kistler).Data;
Kistler_speed=s.(StrName).Y(index_speed_kistler).Data;

% extract time portion 
time_start_fft  = 11;  % 4 / 11 / 22  (RPM = 600 / 1200 / 3000)
duration = 3;
ind = (time>=time_start_fft) & (time<=time_start_fft+duration);
time = time(ind);
vibr_main = vibr_main(ind);
SEW_speed = SEW_speed(ind);
Kistler_speed = Kistler_speed(ind);

speed_RPM_estimated = mean(SEW_speed);  

figure(1)
subplot(2,1,1),plot(time,vibr_main);
subplot(2,1,2),plot(time,SEW_speed,time,Kistler_speed);

%% FFT plot 
nfft = 1024*8;    % 
Overlap = 0.75; % (75% ) overlap 
df = Fs/nfft;

% vibr_main = sin(2*pi*10*df*time); % a simple test with a unitary amplitude sine wave

[freq,fft_spectrum] = myfft_peak(vibr_main, Fs, nfft, Overlap);

figure(2)
freq_normalized = freq*60/speed_RPM_estimated;
plot(freq_normalized,fft_spectrum)
title('FFT Spectrum')
ylabel('Amplitude')
xlabel('Frequency [normalized]')
xlim([0 25]);

figure(3) % check signal stationnarity with spectrogram : OK
nfft = 1024;    % 
Overlap = 0.75; % (75% ) overlap 
spectrogram(vibr_main,hanning(nfft),round(Overlap*nfft),nfft,Fs,'yaxis','minthreshold',-80)  % WINDOW,NOVERLAP,NFFT,Fs
colormap('jet')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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);

% transpose signal dimensions if channels >> samples
if channels > samples
    signal = signal';
    [samples,channels] = size(signal);
end

% 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

[Eduardo]

Hi its me again, is there a way around to define the Fs with a fixed value instead of the formula used there?

[Mathieu NOE]

hello

sure

yiu can alway comment / remove these lines

dt = mean(diff(time));
Fs = 1/dt;

and replace them with Fs = whatever value (as you did in your code)