How to measure the peaks in my plot?

5 次查看(过去 30 天)
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?
  2 个评论
Rik
Rik 2017-2-21
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
Denny Muttathil 2017-2-21
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?

请先登录,再进行评论。

采纳的回答

Rik
Rik 2017-2-21
[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)));

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Spectral Measurements 的更多信息

标签

尚未输入任何标签。

Community Treasure Hunt

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

Start Hunting!

Translated by