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

9 次查看(过去 30 天)
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')
  5 个评论
Mathieu NOE
Mathieu NOE 2021-1-20
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
matteo bottoni 2021-1-21
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

请先登录,再进行评论。

回答(0 个)

Community Treasure Hunt

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

Start Hunting!

Translated by