How can I find a threshold of amplitude and (in the same time) of frequency inside a spectrogram?

I am analyzing 2 vectors of acceleration, ax, az in the frequency domain. They are vectors with sample frequency of 5 ms and gone through a pass band filter of 40 Hz.

My aim is to detect the spectrogram and identify the index of vector ax and/or az that has a certain threshold in amplitude and frequency both in the final spectrogram. I would like to put an alarm (identified with an asterisk) when I have a threshold of the power in a certain range of frequency.

Acc=[ax,az];
Fs=1/0.005;
N=length(az);
%p = nextpow2(N)
figure
nfft=80;
win=hamming(nfft);
nOvl=nfft*0.95;
[s,f,t,pxx] = spectrogram(az,win,nOvl,nfft,Fs,'psd');
waterfall(f,t,pxx')
xlabel('frequency(Hz)')
ylabel('time(sec)')
zlabel('PSD')

[Mathieu NOE]

hello

I modified a bit your code so it matches my visual preferences . Now looking at the spectrogram , what kind of thing do you want to implement when you say my aim is to detect the spectrogram and identify the index of vector ax and/or az that has a certain threshold in amplitude and frequency both.

data = readmatrix('ax_az.xlsx');
ax = data(:,1);
az = data(:,2);

signal = az; % select ax or az

Fs=1/0.005;
N=length(signal);
%p = nextpow2(N)

figure(1)
plot(signal);
xlabel('samples')
ylabel('amplitude')

nfft=256;
win=hanning(nfft);
nOvl=fix(nfft*0.99);
[s,f,t,pxx] = spectrogram(detrend(signal),win,nOvl,nfft,Fs,'psd'); % NB  usage of detrend to avoid spectrogram plot being dominated by DC component of signal

% % FFT normalisation and conversion amplitude from linear to dB (peak)
s_dBpeak = 20*log10(abs(s))+20*log10(2/length(f));     % NB : X=fft(x.*hanning(N))*4/N; % hanning only
% 
% saturation of the dB range : 
spectrogram_dB_scale = 80;  % dB range scale (means , the lowest displayed level is XX dB below the max level)
min_disp_dB = round(max(max(s_dBpeak))) - spectrogram_dB_scale;
s_dBpeak(s_dBpeak<min_disp_dB) = min_disp_dB;

% % plots spectrogram
figure(2)

imagesc(t,f,s_dBpeak);
set(gca,'YDir','normal')
xlabel('time(sec)')
ylabel('frequency(Hz)')

zlabel('dB peak')
colormap('jet');
colorbar('vert');

[matteo bottoni]

Thank you for the interest. Firstly i ask u if in your opinion is good to take an hanning window. With which kind of criteria I can choose the input parameter for spectrogram?

By the way, if I run your code, that returns only a different shape of my same spectrogram, and i plot in a figure(2) the vector ax for example, what i want is to highlight the various coordinates of ax in figure(2) corresponding at least at 4 as power value in a range of 10-30Hz (they are simple values only to explain my aim)

[Mathieu NOE]

hello Matteo

hanning windows are a good start , but the choice of the window depends on the nature of the signal.

Windows are used to make the signal more "periodic" at each buffer , as this is a key specification for make proper fft analysis. there is a lot of information on which window is more appropriate for which type of signal - like the doc in the attachement.

After that , the buffer size (= the size of the fft) will give the frequency resolution (df = Fs / nfft), but the time resolution of the spectrogram is inversely proportionnal to nfft. So increasing nfft improves the freq resolution, but at the price of a reduced time resolution. there is no free lunch here.

sometimes , if your signal contains sharp peaks, and you are looking at precise time events, it's better to work on the time data rather than the spectrogram itself.

for further processing, this is a code example to extend what I already sent you :

% selection 
% max  value of spectrogram dB amplitude in freq range 10 to 30 Hz
ind = find(f> 10 & f<30);
s_dBpeak_10_30Hz_max = max(s_dBpeak(ind,:),[],1);

threshold = 0; % in dB
ind = find(s_dBpeak_10_30Hz_max>=threshold);
time_selected = t(ind);
s_dBpeak_10_30Hz_max_selected = s_dBpeak_10_30Hz_max(ind);


figure(3);plot(t,s_dBpeak_10_30Hz_max,'b',time_selected,s_dBpeak_10_30Hz_max_selected,'.r');
xlabel('time(sec)')
ylabel('max dB amplitude in 10/30 Hz')

[Mathieu NOE]

hello again

I wonder why we do it so complicated...

why not simply apply the band pass filter (10 - 30 hz) to your time data and look at the positive and negative spikes ?

this would be certainly sipmler and more accurate in time than going through the spectrogram stuff

[matteo bottoni]

In the end I figured out with this script. I preferred to continue with the script that gives me a surface like output.

figure
Fs=1/0.005;
nfft=80;
win=hamming(nfft);
nOvl=nfft*0.95;
signal=ax;
[s,f,t,pxx] = spectrogram(signal,win,nOvl,nfft,Fs,'psd');
waterfall(f,t,pxx')
xlabel('frequency(Hz)')
ylabel('time(sec)')
zlabel('PSD')
xlim([0,50])
index=find(f>=10 & f<=40);      %btw 10-40Hz
pxx40Hz=pxx(index(1):index(end),:);
t_pxx40Hz=[t;pxx40Hz];
Times_stumbles= time_of_stumbles (t_pxx40Hz)

Using this function

function time_stumble=time_of_stumbles(Mat)
j=1;
time_stumble=[];
for i=1:size(Mat,2)
    thr= find(Mat(2:end,i)>=1.5);
    if ~isempty (thr)
        time_stumble(j)=Mat(1,i);
        j=j+1;
    end
end
end
        

And then:

time_vector=0:0.005:size(signal,1);
time_vector=time_vector(1:size(signal,1));
plot(time_vector,KinMatrix(:,7))
hold on
for j=1:size(Times_stumbles,2)
    xline(Times_stumbles(j),'-.r','LineWidth',1)
end