How to measure the peaks in my plot?

I will keep the explanation, how my codes works, very short. I advise you to try this code yourself, so perhaps you understand it better that way. I have an audio-file and read it in my code. Now i switch from the time domain to the frequency domain by using the function FFT. But the only difference is, that i am performing an STFT on my audio signal. I do it every 30ms, until to the length of my signal. I am aware, that there are many different function in matlab, which also can perform this easily, but there are not giving me the results i need. Now, i am plotting many different frequency spectrums every 30ms. But i split up my signal in three frequency bands. They are called LOW, MEDIUM and HIGH. Basically, this means I have 3 different spectrums plotting every 30ms. The next step I do, is summing all the magnitudes from ONE frequency spectrum together, this means I have ONE VALUE per frequency spectrum, which are being squared.

Now, i have the power from every spectrum ! And all of these values are being plotted in my code. I am only plotting the power values, otherwise my code performance time would be extremely slow. Btw, the code looks long, but there are two for loop. In the first, i read the low spectrum and when it is finished, the second starts with the medium and high spectrum. Basically they are the same. I am aware, i can probably do that with findpeaks or something similar. But how can i do that?

clear;
clc;
%%MATLAB
%%read file
%_________________________________________
[y,fs]=audioread('Undertale - Megalovania.wav');
% audioread = read wav -file
% y = contains the audio signal
% fs = 44100
% 'UnchainMyHeart' = name of the wav-file
%_________________________________________
%%PARAMETER FOR STFT
%_________________________________________ 
t_seg=0.03; % length of segment in ms
fftlen = 4096; %FFT-Points

% Defining size of frequency bands
f_low= 1:200;    %lower frequencies
f_medium= 201:600;  %medium frequencies
f_high= 601:1000; %higher frequencies
%__________________________________________
%%CODE

segl =floor(t_seg*fs); 
windowshift=segl/2; 
% defining the size of the window shift
window=hann(segl); 
% apply hann function on segment length (30 ms)
window=window.'; 
% transpose vector
si=1; 
% defining start index
ei=segl; 
% defining end index

N=floor( length(y)/windowshift - 1);
% Calculates the number, how often the window has to shift
% until to length of the audio signal

f1=figure;
    % Generating new window

      f=0:1:fftlen-1;
      f=f/fftlen*fs;
      % defining frequency vector

      Ya=zeros(1,fftlen);

      ValuesOfYc = NaN(1,N);
      ValuesOfYd = NaN(1,N);
      ValuesOfYe = NaN(1,N);

      x =(1:N)*windowshift/fs;
      % defining x-axis

for m= 1:1:N 

      y_a = y(si:ei);
      % a segment is taken out from audio signal length(30ms) 
      y_a= y_a.*window;
      % multiplying segment with window (hanning) 
      Ya=fft(y_a, fftlen);
      % Applying fft on segment
      Yb=abs(Ya(1:end/2)).^2;  
      % Squaring the magnitudes from one-sided spectrum

          drawnow; % Updating the graphical values

      figure(f1);
      % Showing the power values

      %%frequency bands

      y_low = Yb(f_low);  % LOW frequency spectrum

      Yc=sum(y_low);
      % Summing all the power values from one frequency spectrum together
      % so you get one power value from one spectrum 
      ValuesOfYc(m) = Yc;
      %Output values are being saved here, which are generated from the for
      %loop
      % m = start variable from for loop

      subplot(2,1,1)
      p=plot(x,ValuesOfYc,'r-');%,x, ValuesOfYd,'g-', x, ValuesOfYe,'b-' );  
      p(1).LineWidth =0.5;
      xlabel('time (Audio length)')
      ylabel('Power')
      grid on

      si=si+windowshift; 
      % Updating start index 
      ei=ei+windowshift; 
      % Updating end index

end

for o= 1:1:N

      y_a = y(si:ei);
      % a segment is taken out from audio signal length(30ms) 
      y_a= y_a.*window;
      % multiplying segment with window (hanning) 
      Ya=fft(y_a, fftlen);
      % Applying fft on segment
      Yb=abs(Ya(1:end/2)).^2;  
      % Squaring the magnitudes from one-sided spectrum

          drawnow; % Updating the graphical values

      figure(f1);
      % Showing the power values

      %%frequency bands

      y_medium = Yb(f_medium); % MEDIUM frequency spectrum
      y_high = Yb(f_high); % HIGH frequency spectrum

      Yd=sum(y_medium);
      Ye=sum(y_high);
      % Summing all the power values from one frequency spectrum together
      % so you get one power value from one spectrum 
      ValuesOfYd(o) = Yd;
      ValuesOfYe(o) = Ye;
      %Output values are being saved here, which are generated from the for
      %loop
      % m = start variable from for loop

      subplot(2,1,2)
      p=plot(x, ValuesOfYd,'g-', x, ValuesOfYe,'b-' );
      p(1).LineWidth =0.5;
      xlabel('time (Audio length)')
      ylabel('Power')
      grid on

      si=si+windowshift; 
      % Updating start index 
      ei=ei+windowshift; 
      % Updating end index
  end

I hope you can read the PDF-file. In my plot(the red one), i want to measure the peaks (basically the coordinates) and the distance between them. How can I do that?

[Rik]

Have you tried findpeaks yet? I just scrolled through the doc and I noticed you can even use the 'MinPeakProminence' switch to only get the larger peaks, so you get the true peaks. I would suggest setting that value to a percentage of the max-min.

[Denny Muttathil]

Hello, thank you for trying to help me. I already used findpeaks. It worked, but not exactly the result, I wanted. Btw, what exactly is "true peaks" meant to be? Is it perhaps possible, to measure the peaks, which are above at the 5000 line in my red plot. If yes, any suggestion, how to go on there?

Answer 1

[Rik]

[continued from comments]

With true peaks I meant the larger peaks, not just every local maximum. If you use the prominence you can filter them out easily, but you can use the strategy you propose.

[~,locs] = findpeaks(data)
real_peak_locs=intersect(find(data>5000),locs);
pks=data(real_peak_locs);

The variable pks will contain the value for all points that are a local maximum whose value is greater than 5000.

Although I still think this syntax is more elegant:

 [pks,real_peak_locs] = findpeaks(data,'MinPeakProminence',0.75*(max(data)-min(data)));

[Denny Muttathil]

Ok, thank you very much, now i got it!