You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
How to set Spectrogram parameters
10 views (last 30 days)
Show older comments
Hello Everyone,
I would like to print the spectrogram of a curve, in order to observe variations. But I don't really understand how to set the parameters of the function "spectrogram" even with the Matlab help.
For now I have this :
But I can't observe any variation. Is it possible to obtain something more like that :
There is my code and the attached file to get all the datas :
clear all;
clc ;
close all;
Valeurs = xlsread('C:\Users\tlam\Desktop\Run 80.xlsx',"Data1");
Temps = Valeurs(:,1);
Accelero = Valeurs(:,2);
Moteur_rpm = Valeurs(:,3);
BV_rpm = Valeurs(:,4);
Delta = Moteur_rpm - BV_rpm;
M_temps = max(Temps);
Moteur_Hz = Moteur_rpm * 0.016667;
BV_Hz = BV_rpm * 0.016667 ;
Delta_Hz = Delta * 0.016667;
%%
F = figure( 'Position', [50 50 1600 900]) ; % Visible 'off' pour ne pas afficher
t = tiledlayout(3,3,'TileSpacing','Compact','Padding','Compact');
title(t,"Spectro FULL ")
% Mot_rpm(t) et BV_rpm(t) Run 1
nexttile([1 2])
plot(Temps,Moteur_rpm,'red',Temps,BV_rpm,'green')
xlim ([0 M_temps])
ylabel('Moteur & BV [rpm]')
legend ('Moteur','Boite de Vitesse')
grid on
Fs = 50000; %Hz
% Spectrogramme Run 1
nexttile([2 1])
spectrogram(Accelero,100,99,100,Fs,'yaxis','power')
hold on
plot(Temps,Moteur_Hz,'r',Temps,BV_Hz,'g',Temps,Delta_Hz,'b')
colormap('jet')
% Delta(t) Run 1
nexttile([1 2])
plot(Temps,Delta,'blue')
xlim ([0 M_temps])
ylabel('Delta [rpm]')
grid on
% Accelero(t) Run 1
nexttile([1 2])
plot(Temps,Accelero,'magenta')
xlabel('temps [s]')
xlim ([0 M_temps])
ylabel('acceleromètre [V]')
grid on
19 Comments
Mathieu NOE
on 9 May 2023
Hello
I played a bit with your code and data and tried to make a better spectrogram rendering
I added some high pass filetring on your accel data, otherwise the spectrogram has a very strong output at DC 'not very interesting !)
now you can play with NFFT and OVERLAP if you want to have different frequency and time resolution
also I noticed that inside your code you hardcoded Fs = 5000, whereas from the dat file it's 50 Hz, what is the correct value ?
I could not find any visual correlation between your accel spectrogram and the other signals (converted to Hz) but I wonder if it's correct to compare this way as it assumes that you are only looking at the first engine order (when you convert from RPM to Hz).
For an IC engine, the vibration are higher at multiple engine orders (2nd and above typically)
code a bit tweaked :
clear all;
clc ;
close all;
% Valeurs = xlsread('C:\Users\tlam\Desktop\Run 80.xlsx',"Data1");
Valeurs = xlsread('Run 1.xlsx',"Data1");
Temps = Valeurs(:,1);
dt = mean(diff(Temps));
Fs = 1/dt;
Accelero = Valeurs(:,2);
Moteur_rpm = Valeurs(:,3);
BV_rpm = Valeurs(:,4);
Delta = Moteur_rpm - BV_rpm;
M_temps = max(Temps);
Moteur_Hz = Moteur_rpm /60;
BV_Hz = BV_rpm /60 ;
Delta_Hz = Delta /60;
% high pass filtering of accel signal
f_high = 1; % Hz
order = 3
[b,a] = butter(3,2/Fs*f_high,'high');
Accelero_filtered = filtfilt(b,a,Accelero);
%%
% F = figure( 'Position', [50 50 1600 900]) ; % Visible 'off' pour ne pas afficher
F = figure() ; % Visible 'off' pour ne pas afficher
t = tiledlayout(3,3,'TileSpacing','Compact','Padding','Compact');
title(t,"Spectro FULL ")
% Mot_rpm(t) et BV_rpm(t) Run 1
nexttile([1 2])
plot(Temps,Moteur_rpm,'red',Temps,BV_rpm,'green')
xlim ([0 M_temps])
ylabel('Moteur & BV [rpm]')
legend ('Moteur','Boite de Vitesse')
grid on
%Fs = 50000; %Hz ???
% Spectrogramme Run 1
nexttile([2 1])
NFFT =200;
NOVERLAP = round(0.95*NFFT);
spectrogram(Accelero_filtered,hanning(NFFT),NOVERLAP,NFFT,Fs,'yaxis','power','MinThreshold',-80)
hold on
plot(Temps,Moteur_Hz,'r',Temps,BV_Hz,'g',Temps,Delta_Hz,'b')
colormap('jet')
% Delta(t) Run 1
nexttile([1 2])
plot(Temps,Delta,'blue')
xlim ([0 M_temps])
ylabel('Delta [rpm]')
grid on
% Accelero(t) Run 1
nexttile([1 2])
plot(Temps,Accelero,'magenta')
xlabel('temps [s]')
xlim ([0 M_temps])
ylabel('acceleromètre [V]')
grid on
Mathieu NOE
on 9 May 2023
as you increase the sampling rate , the frequency resolution df will decrease if you keep NFFT the same (because df = Fs/NFFT)
so you may have to increase it significantly
now for overlap , i rarely go above 75% of overlap if the transients are not super brutal (and you need a refined time resolution) , so you shouldn't have to change it all the time
Mathieu NOE
on 9 May 2023
try also to shift the high pass filter cut off frequency to a higher value (~ 1 kHz ?)
Mathieu NOE
on 9 May 2023
I'm curious
why did you apply a factor 1000 on the green and red frequency curves (that is like you are plotinng engine 1000th order frequency ?)
Mathieu NOE
on 9 May 2023
which brings me another remark from your spectrogram , as the energy sems to be vanishing above 1 kHz
what frequency range are you indeed looking for ?
Théo
on 9 May 2023
In the Excel file, which I am working currently the acquisition frequency is 50k Hz. I ploted the raw data and divided by 60 to get the value in Hz. I didn't apply a factor 1000, did I ?
I am not looking for a specific frequency, In fact I want to observe the judder on the spectrogram to determine if it's a geometrical judder :
or a "normal" one
Sorry I don't explain very well, I see what I am supposed to get but I can't explain it very well.
I tried to shift the high pass filter cut off frequency to a higher value :
% high pass filtering of accel signal
f_high = 1000; % Hz
Mathieu NOE
on 9 May 2023
I am a bit lost why you have either 50 Hz or 50 kHz sampling rate
are you measuring two distinct engines (one slow and one fast) ? or are you measuring the same engine with two different recording parameters ?
I just don't get why in the second case (Fs = 50 kHz), on the spectrogram , the green, blue and red curves shows values in the kHz range , whereas the RPM plots (data) seems the same as in the first case (Fs = 50 Hz)
I'm talking about these data :
Moteur_Hz = Moteur_rpm /60;
BV_Hz = BV_rpm /60 ;
Delta_Hz = Delta /60;
anyway it seems the spectrogram shows no energy that is close more or less to the 3 curves (only background noise)
Théo
on 10 May 2023
That's the same engine with two recording parameters..
The values are supposed to be the same. And you are right I dont really understand, suddenly I have got values in the kHz range on the spectrogram..
That's my code :
Valeurs = xlsread('C:\Users\tlam\Desktop\Run 80.xlsx',"Data1");
Temps = Valeurs(:,1);
Accelero = Valeurs(:,2);
Moteur_rpm = Valeurs(:,3);
BV_rpm = Valeurs(:,4);
Delta = Moteur_rpm - BV_rpm;
M_temps = max(Temps);
Moteur_Hz = Moteur_rpm /60;
BV_Hz = BV_rpm /60 ;
Delta_Hz = Delta /60;
I didn't apply any factor and the values between 50Hz and 50kHz are the same...
It works well for 50 Hz but not for 50 kHz strange ..
Mathieu NOE
on 10 May 2023
Edited: Mathieu NOE
on 10 May 2023
can you share your last version of the code with an extract of the 50 kHz data ? I don't need the full length of data
Théo
on 10 May 2023
Yea sure :
clear all;
clc ;
close all;
Valeurs = xlsread('C:\Users\tlam\Desktop\Run 80bis.xlsx',"Data1");
Temps = Valeurs(:,1);
Accelero = Valeurs(:,2);
Moteur_rpm = Valeurs(:,3);
BV_rpm = Valeurs(:,4);
Delta = Moteur_rpm - BV_rpm;
M_temps = max(Temps);
Moteur_Hz = Moteur_rpm /60;
BV_Hz = BV_rpm /60 ;
Delta_Hz = Delta /60;
dt = mean(diff(Temps));
Fs = 1/dt; %Hz
% high pass filtering of accel signal
f_high = 1000; % Hz
order = 3;
[b,a] = butter(3,2/Fs*f_high,'high');
Accelero_filtered = filtfilt(b,a,Accelero);
%%
F = figure( 'Position', [50 50 1600 900]) ; % Visible 'off' pour ne pas afficher
t = tiledlayout(3,3,'TileSpacing','Compact','Padding','Compact');
title(t,"Spectro FULL ")
% Mot_rpm(t) et BV_rpm(t) Run 1
nexttile([1 2])
plot(Temps,Moteur_rpm,'red',Temps,BV_rpm,'green')
xlim ([0 M_temps])
ylabel('Moteur & BV [rpm]')
legend ('Moteur','Boite de Vitesse')
grid on
nexttile([2 1])
NFFT =200;
NOVERLAP = round(0.95*NFFT);
spectrogram(Accelero_filtered,hanning(NFFT),NOVERLAP,NFFT,Fs,'yaxis','power','MinThreshold',-80)
hold on
plot(Temps,Moteur_Hz,'r',Temps,BV_Hz,'g',Temps,Delta_Hz,'b')
colormap('jet')
% Delta(t) Run 1
nexttile([1 2])
plot(Temps,Delta,'blue')
xlim ([0 M_temps])
ylabel('Delta [rpm]')
grid on
% Accelero(t) Run 1
nexttile([1 2])
plot(Temps,Accelero,'magenta')
xlabel('temps [s]')
xlim ([0 M_temps])
ylabel('acceleromètre [V]')
grid on
I linked a file with the first values at 50 kHz rate.
Thanks again Mathieu, taking time to help me that's realy kind :)
Mathieu NOE
on 10 May 2023
hello
try this new code
I let you discover the new spectrogram (with frequency axis in log scale to have good rendering both low and high freqs)
now the engine frequency trace is to be seen
by again I see no special event that coincide with the engine freq line
clear all;
clc ;
close all;
Valeurs = xlsread('Run 80bis.xlsx',"Data1");
Temps = Valeurs(:,1);
Accelero = Valeurs(:,2);
Moteur_rpm = Valeurs(:,3);
BV_rpm = Valeurs(:,4);
Delta = Moteur_rpm - BV_rpm;
M_temps = max(Temps);
Moteur_Hz = Moteur_rpm /60;
BV_Hz = BV_rpm /60 ;
Delta_Hz = Delta /60;
dt = mean(diff(Temps));
Fs = 1/dt; %Hz
% high pass filtering of accel signal
f_high = 5; % Hz
order = 2;
[b,a] = butter(3,2/Fs*f_high,'high');
Accelero_filtered = filtfilt(b,a,Accelero);
%%
F = figure( 'Position', [50 50 1600 900]) ; % Visible 'off' pour ne pas afficher
t = tiledlayout(3,3,'TileSpacing','Compact','Padding','Compact');
title(t,"Spectro FULL ")
% Mot_rpm(t) et BV_rpm(t) Run 1
nexttile([1 2])
plot(Temps,Moteur_rpm,'red',Temps,BV_rpm,'green')
xlim ([0 M_temps])
ylabel('Moteur & BV [rpm]')
legend ('Moteur','Boite de Vitesse')
grid on
nexttile([2 1])
NFFT =20000;
NOVERLAP = round(0.95*NFFT);
[S,F,T] = spectrogram(Accelero_filtered,hanning(NFFT),NOVERLAP,NFFT,Fs,'yaxis','power');
idf = (F>0); % remove zero (dc) freq
F = F(idf);
S = S(idf,:);
SdB = 20*log10(abs(S)); % convert S in dB
SdBmax = max(SdB,[],'all');
SdB_range = 80; % range in dB to display
SdB(SdB<SdBmax-SdB_range) = NaN; % remove data below threshold (corresponding to max value - SdB_range)
imagesc(T,F,SdB);
set(gca, 'YDir','normal'); % to have the y = 0 at the bottom and not at the top
set(gca, 'YScale','log'); %
colorbar('vert');
xlabel('temps [s]')
xlim ([0 M_temps])
ylabel('Freq (Hz)')
hold on
plot(Temps,Moteur_Hz,'r',Temps,BV_Hz,'g',Temps,Delta_Hz,'b')
colormap('jet')
% Delta(t) Run 1
nexttile([1 2])
plot(Temps,Delta,'blue')
xlim ([0 M_temps])
ylabel('Delta [rpm]')
grid on
% Accelero(t) Run 1
nexttile([1 2])
plot(Temps,Accelero,'magenta')
xlabel('temps [s]')
xlim ([0 M_temps])
ylabel('acceleromètre [V]')
grid on
Théo
on 11 May 2023
Edited: Théo
on 11 May 2023
Hello Mathieu,
I tried and get this :
I was supposed to work on the files with the acquisition freq of 50 kHz but the process time and the results are not convincing.
May I keep on the 50Hz even if the resolution of the spectrogram isn't the best.
Anyway thanks Mathieu for taking all this time for me, I have one last question to which I think I have the answer :
I tried the first code you did on an other vehicle test (freq : 50 Hz) And the curves on the spectrogram dont fit in the window.
I guess it's normal because the input signal in the spectrogram function is the accelerometer, then the engine curve wont necessarily fit. The windowing depends exclusively of the accelerometer.
Can you confirm I am right? Thanks
Mathieu NOE
on 11 May 2023
to you last question :
yes, the spectrogram max displayed frequency is Fs / 2 = 25 Hz
so any RPM trace above 25*60 = 1500 RPM will not appear
all the best for the future !
Mathieu NOE
on 11 May 2023
My pleasure !
Answers (0)
See Also
Categories
Find more on Multirate Signal Processing in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom(English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)