calculating the x1 and x2 of a peak automatically for every data ?

Question

VG 2022-3-30

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1684289-calculating-the-x1-and-x2-of-a-peak-automatically-for-every-data

编辑： VG 2022-4-4

Hello Everyone,

I would like to know how to calculate the x1 and x2 (width) of a maximum peak in frequency domain? I have tried findpeaks and intep but it did not work. please someone help me to find out the x1 and x2 (marked in red in the below image)of the maximum peak automatically for every data (raw data when I upload as csv.file) . I tried many methods and examples but its not working for my data. I could sucessufully calculate and findout the maximum peak for every data but not the x1 and x2.

Below I am attaching the picture and my code for you reference. Hoping for a postive reply.

clc ;

clearvars;

close all;

format compact ;

fontSize = 24;

%load in the data

data = readtable('data.csv');

%convert x and sensor from table to ARRAY

time(:,1) =table2array(data(:,1)) ;

sensor(:,2) =table2array(data(:,2));

% converting the double array to single array

sensor2 = sensor(:,2);

%time domain

plot(time,sensor2);

t =time(29:600); % trimming the vector

s= sensor2(29:600);

plot(t,s);

% frequency domain

df = 1/(t(end)-t(1)) ; % damping factor

N = numel(t); % signal length (number of sample )

f = 0:df:(N-1)*df ; %frequency vector

S = fft(s) ; %fourier transform

% plot(f,abs(fft(s))); % frequency domain

plot(f,abs(S)); % plotting in frequency domain

S(2:4) =0;

plot(f,abs(S))

S(1:6) =0 ;

plot(f,abs(S));

% Q = 0.01173 /(0.01243-0.01103)

[x ,y] = max(abs(S));

% x =

% 0.0783

% y =

% 571

fc = S(y); % center frequency

% fc =

% -0.0643 - 0.0446i

Sabs = abs(S);

fc = Sabs(y);

% fc =

% 0.0783

% f1 = f(1:N/2)

f1= f(1:N/2) ;

S1 = Sabs(1:N/2); % changing to the total length

B = smoothdata(S1);

plot(f1,S1);

plot(f1,B);

xlabel("FrequencY(Hz) ") ;

ylabel("Amplitude")

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Mathieu NOE 2022-3-30

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1684289-calculating-the-x1-and-x2-of-a-peak-automatically-for-every-data#answer_930659

在 MATLAB Online 中打开

hello

tweaked and completed your code so that the Q factor is now automatically computed

x1 is flow in my code, x2 is fhigh , and Q = fc/(fhigh - flow)

Q = 21.1344

the center peak is taken on the raw fft while it's better to pick x1 and x2 on the smoothed curve (to get unique values)

hope it helps

clc ;
clearvars;
close all;
format compact ;
fontSize = 24;
%load in the data
data = readtable('data.csv');
%convert x and sensor from table to ARRAY 
time(:,1) =table2array(data(:,1)) ;
sensor(:,2) =table2array(data(:,2));
% converting the double array to single array
sensor2 = sensor(:,2);
%time domain 
plot(time,sensor2);
t =time(29:600); % trimming the vector 
s= sensor2(29:600);
plot(t,s);
%  frequency domain
df = 1/(t(end)-t(1)) ; % damping factor
N = numel(t);  % signal length (number of sample )
f = 0:df:(N-1)*df ;   %frequency vector
Sabs = abs(fft(s)) ;          %fourier transform
% upper half frequency range
f= f(1:N/2) ;
Sabs = Sabs(1:N/2); % changing to the total length 
% plot(f,abs(fft(s))); % frequency domain 
plot(f,Sabs); % plotting in frequency domain 
Sabs(1:20) =0;
[Smax ,ind_max] = max(Sabs);
fc = f(ind_max); % peak frequency 
B = smoothdata(Sabs,'gaussian',10);
plot(f,Sabs,f,B); 
xlabel("FrequencY(Hz) ") ;
ylabel("Amplitude")
% find Q from the -3 dB rule
Smax_3dB = 0.707*Smax;
threshold = Smax_3dB; % your value here
[flow,s0_pos1,fhigh,s0_neg1]= crossing_V7(B,f,threshold,'linear'); % positive (pos) and negative (neg) slope crossing points 
% ind => time index (samples)
% t0 => corresponding time (x) values 
% s0 => corresponding function (y) values , obviously they must be equal to "threshold"
plot(f,Sabs,f,B,fc,Smax,'dk',flow,s0_pos1,'dr',fhigh,s0_neg1,'dr'); 
xlabel("FrequencY(Hz) ") ;
ylabel("Amplitude")
Q = fc/(fhigh - flow)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [t0_pos,s0_pos,t0_neg,s0_neg] = crossing_V7(S,t,level,imeth)
% [ind,t0,s0,t0close,s0close] = crossing_V6(S,t,level,imeth,slope_sign) % older format
% CROSSING find the crossings of a given level of a signal
%   ind = CROSSING(S) returns an index vector ind, the signal
%   S crosses zero at ind or at between ind and ind+1
%   [ind,t0] = CROSSING(S,t) additionally returns a time
%   vector t0 of the zero crossings of the signal S. The crossing
%   times are linearly interpolated between the given times t
%   [ind,t0] = CROSSING(S,t,level) returns the crossings of the
%   given level instead of the zero crossings
%   ind = CROSSING(S,[],level) as above but without time interpolation
%   [ind,t0] = CROSSING(S,t,level,par) allows additional parameters
%   par = {'none'|'linear'}.
%	With interpolation turned off (par = 'none') this function always
%	returns the value left of the zero (the data point thats nearest
%   to the zero AND smaller than the zero crossing).
%
% check the number of input arguments
error(nargchk(1,4,nargin));
% check the time vector input for consistency
if nargin < 2 | isempty(t)
	% if no time vector is given, use the index vector as time
    t = 1:length(S);
elseif length(t) ~= length(S)
	% if S and t are not of the same length, throw an error
    error('t and S must be of identical length!');    
end
% check the level input
if nargin < 3
	% set standard value 0, if level is not given
    level = 0;
end
% check interpolation method input
if nargin < 4
    imeth = 'linear';
end
% make row vectors
t = t(:)';
S = S(:)';
% always search for zeros. So if we want the crossing of 
% any other threshold value "level", we subtract it from
% the values and search for zeros.
S   = S - level;
% first look for exact zeros
ind0 = find( S == 0 ); 
% then look for zero crossings between data points
S1 = S(1:end-1) .* S(2:end);
ind1 = find( S1 < 0 );
% bring exact zeros and "in-between" zeros together 
ind = sort([ind0 ind1]);
% and pick the associated time values
t0 = t(ind); 
s0 = S(ind);
if ~isempty(ind)
    if strcmp(imeth,'linear')
        % linear interpolation of crossing
        for ii=1:length(t0)
            %if abs(S(ind(ii))) >= eps(S(ind(ii)))    % MATLAB V7 et +
            if abs(S(ind(ii))) >= eps*abs(S(ind(ii)))    % MATLAB V6 et -    EPS * ABS(X)
                % interpolate only when data point is not already zero
                NUM = (t(ind(ii)+1) - t(ind(ii)));
                DEN = (S(ind(ii)+1) - S(ind(ii)));
                slope =  NUM / DEN;
                slope_sign(ii) = sign(slope);
                t0(ii) = t0(ii) - S(ind(ii)) * slope;
                s0(ii) = level;
            end
        end
    end
    % extract the positive slope crossing points 
    ind_pos = find(sign(slope_sign)>0);
    t0_pos = t0(ind_pos);
    s0_pos = s0(ind_pos);
    % extract the negative slope crossing points 
    ind_neg = find(sign(slope_sign)<0);
    t0_neg = t0(ind_neg);
    s0_neg = s0(ind_neg);
else
    % empty output
    ind_pos = [];
    t0_pos = [];
    s0_pos = [];
    % extract the negative slope crossing points 
    ind_neg = [];
    t0_neg = [];
    s0_neg = [];
    
end
end

7 个评论
显示 5更早的评论隐藏 5更早的评论

VG 2022-3-31

data_4.1.csv

Hello Mathieu NOE

Thank you very much for your help.I have two final question regarding the Q value.

Question:1 When I am uploading different data set , somehow the Q value is not calculated automatically and its picking several x values from the next maxima. But it worked perfectly for the above mentioned data set. I can understand that the threshold was set to 0.707, if I am not wrong .Now I would like to calculate the Q value for any data set below the threshold , will that be possible?

Note :From the above data set the time was in micro seconds and I have changed it into seconds.

Question 2 : I also want to check if this raw data is fitting with a gaussian curve and compare the center frequency with the raw data. Gaussian function formaula that I want to use A = a exp( - b (f-fc)^2) and later find out the center frequency of the gaussian curve.

Below I am uploading the figures for different data sets ands also the raw data for your reference. Hoping for a postive reply.

Thank you

code for your reference ;

clc ;

clearvars;

close all;

format compact ;

fontSize = 24;

%load in the data

data = readtable('data4.1.csv');

%convert x and sensor from table to ARRAY

time(:,1) =table2array(data(:,1)) ;

sensor(:,2) =table2array(data(:,2));

% converting the double array to single array

sensor2 = sensor(:,2);

%time domain

plot(time,sensor2);

t = time/10000 ;

T = t(29:600); % trimming the vector

s = sensor2(29:600);

% frequency domain

df = 1/(T(end)-T(1)) ; % damping factor

N = numel(T); % signal length (number of sample )

f = 0:df:(N-1)*df ; %frequency vector

Sabs = abs(fft(s)) ; %fourier transform

% upper half frequency range

f= f(1:N/2) ;

Sabs = Sabs(1:N/2); % changing to the total length

% plot(f,abs(fft(s))); % frequency domain

plot(f,Sabs); % plotting in frequency domain

Sabs(1:20)= 0;

[Smax ,ind_max] = max(Sabs);

fc = f(ind_max); % peak frequency

B = smoothdata(Sabs,'gaussian',10);

plot(f,Sabs,f,B);

xlabel("FrequencY(Hz) ") ;

ylabel("Amplitude");

% find Q from the -3 dB rule

Smax_3dB = 0.707*Smax;

threshold = Smax_3dB; % your value here

[flow,s0_pos1,fhigh,s0_neg1]= crossing_V7(B,f,threshold,'linear'); % positive (pos) and negative (neg) slope crossing points

% ind => time index (samples)

% t0 => corresponding time (x) values

% s0 => corresponding function (y) values , obviously they must be equal to "threshold"

figure(1);

subplot(2,1,1);

plot(T,s);

title({'Outer ring : Air';'Time domain'});

xlabel('Time (sec)');

ylabel('Amplitude');

subplot(2,1,2);

plot(f,Sabs,f,B,fc,Smax,'dk',flow,s0_pos1,'dr',fhigh,s0_neg1,'dr'); hold on

xlabel("FrequencY(Hz) ");

ylabel("Amplitude");

%Calculating the quality factor

Q = fc/(fhigh - flow)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [t0_pos,s0_pos,t0_neg,s0_neg] = crossing_V7(S,t,level,imeth)

% [ind,t0,s0,t0close,s0close] = crossing_V6(S,t,level,imeth,slope_sign) % older format

% CROSSING find the crossings of a given level of a signal

% ind = CROSSING(S) returns an index vector ind, the signal

% S crosses zero at ind or at between ind and ind+1

% [ind,t0] = CROSSING(S,t) additionally returns a time

% vector t0 of the zero crossings of the signal S. The crossing

% times are linearly interpolated between the given times t

% [ind,t0] = CROSSING(S,t,level) returns the crossings of the

% given level instead of the zero crossings

% ind = CROSSING(S,[],level) as above but without time interpolation

% [ind,t0] = CROSSING(S,t,level,par) allows additional parameters

% par = {'none'|'linear'}.

% With interpolation turned off (par = 'none') this function always

% returns the value left of the zero (the data point thats nearest

% to the zero AND smaller than the zero crossing).

%

% check the number of input arguments

error(nargchk(1,4,nargin));

% check the time vector input for consistency

if nargin < 2 | isempty(t)

% if no time vector is given, use the index vector as time

t = 1:length(S);

elseif length(t) ~= length(S)

% if S and t are not of the same length, throw an error

error('t and S must be of identical length!');

end

% check the level input

if nargin < 3

% set standard value 0, if level is not given

level = 0;

end

% check interpolation method input

if nargin < 4

imeth = 'linear';

end

% make row vectors

t = t(:)';

S = S(:)';

% always search for zeros. So if we want the crossing of

% any other threshold value "level", we subtract it from

% the values and search for zeros.

S = S - level;

% first look for exact zeros

ind0 = find( S == 0 );

% then look for zero crossings between data points

S1 = S(1:end-1) .* S(2:end);

ind1 = find( S1 < 0 );

% bring exact zeros and "in-between" zeros together

ind = sort([ind0 ind1]);

% and pick the associated time values

t0 = t(ind);

s0 = S(ind);

if ~isempty(ind)

if strcmp(imeth,'linear')

% linear interpolation of crossing

for ii=1:length(t0)

%if abs(S(ind(ii))) >= eps(S(ind(ii))) % MATLAB V7 et +

if abs(S(ind(ii))) >= eps*abs(S(ind(ii))) % MATLAB V6 et - EPS * ABS(X)

% interpolate only when data point is not already zero

NUM = (t(ind(ii)+1) - t(ind(ii)));

DEN = (S(ind(ii)+1) - S(ind(ii)));

slope = NUM / DEN;

slope_sign(ii) = sign(slope);

t0(ii) = t0(ii) - S(ind(ii)) * slope;

s0(ii) = level;

end

% extract the positive slope crossing points

ind_pos = find(sign(slope_sign)>0);

t0_pos = t0(ind_pos);

s0_pos = s0(ind_pos);

% extract the negative slope crossing points

ind_neg = find(sign(slope_sign)<0);

t0_neg = t0(ind_neg);

s0_neg = s0(ind_neg);

else

% empty output

ind_pos = [];

t0_pos = [];

s0_pos = [];

% extract the negative slope crossing points

ind_neg = [];

t0_neg = [];

s0_neg = [];

end

Mathieu NOE 2022-4-1

在 MATLAB Online 中打开

hello again

I think I have improved the code

try this one with your different data files and tell me how you feel about it

result for last data set :

code :

clc ;
clearvars;
close all;
format compact ;
fontSize = 24;
%load in the data
data = readmatrix('data_4.1.csv');
% data = readmatrix('data.csv');
%convert x and sensor from table to ARRAY 
t =data(:,1) ;
s =data(:,2);
samples = numel(t);
%time domain 
t = t(1:round(samples/3)); % trimming the vector 
s = s(1:round(samples/3));
[m,ind] = max(abs(s)); % find first peak (neg or pos)
t = t(ind-2:end); % trimming the vector 
s = s(ind-2:end);
% apply exponential decay window (will smooth the data) (optionnal)
win = exp(-0.01*(1:numel(t)));
s = s.*win';
% t = time/10000 ; % conversion to seconds
t = t/1e6 ; % conversion micro seconds to seconds
%  frequency domain
df = 1/(t(end)-t(1)) ; % damping factor
N = numel(t);  % signal length (number of sample )
f = 0:df:(N-1)*df ;   %frequency vector
Sabs = abs(fft(s)) ;   %fourier transform
% first half frequency range
f = f(1:N/2) ;
Sabs = Sabs(1:N/2); % changing to the total length 
Sabs(1:20)= Sabs(20);% remove very low frequency peaks
[Smax ,ind_max] = max(Sabs);
fc = f(ind_max); % peak frequency 
B = smoothdata(Sabs,'gaussian',10);
% find Q from the -3 dB rule
 Smax_3dB = 0.707*Smax;
threshold = Smax_3dB; % your value here
[flow,s0_pos1,fhigh,s0_neg1]= crossing_V7(B,f,threshold,'linear'); % positive (pos) and negative (neg) slope crossing points 
% ind => time index (samples)
% t0 => corresponding time (x) values 
% s0 => corresponding function (y) values , obviously they must be equal to "threshold"
% select "flow" value nearest to dominant peak (fc)
[~,ind] = min(abs(flow-fc));
flow = flow(ind);
% select "fhigh" value nearest to dominant peak (fc)
[~,ind] = min(abs(fhigh-fc));
fhigh = fhigh(ind);
figure(2);
subplot(2,1,1);
plot(t,s);
title({'Outer ring : Air';'Time domain'});
xlabel('Time (sec)'); 
ylabel('Amplitude');
subplot(2,1,2);
plot(f,Sabs,f,B,fc,Smax,'dk',flow,threshold,'dr',fhigh,threshold,'dr'); hold on
xlabel("FrequencY(Hz) ");
ylabel("Amplitude");
%Calculating the quality factor
Q = fc/(fhigh - flow)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [t0_pos,s0_pos,t0_neg,s0_neg] = crossing_V7(S,t,level,imeth)
% [ind,t0,s0,t0close,s0close] = crossing_V6(S,t,level,imeth,slope_sign) % older format
% CROSSING find the crossings of a given level of a signal
%   ind = CROSSING(S) returns an index vector ind, the signal
%   S crosses zero at ind or at between ind and ind+1
%   [ind,t0] = CROSSING(S,t) additionally returns a time
%   vector t0 of the zero crossings of the signal S. The crossing
%   times are linearly interpolated between the given times t
%   [ind,t0] = CROSSING(S,t,level) returns the crossings of the
%   given level instead of the zero crossings
%   ind = CROSSING(S,[],level) as above but without time interpolation
%   [ind,t0] = CROSSING(S,t,level,par) allows additional parameters
%   par = {'none'|'linear'}.
%	With interpolation turned off (par = 'none') this function always
%	returns the value left of the zero (the data point thats nearest
%   to the zero AND smaller than the zero crossing).
%
% check the number of input arguments
error(nargchk(1,4,nargin));
% check the time vector input for consistency
if nargin < 2 | isempty(t)
	% if no time vector is given, use the index vector as time
    t = 1:length(S);
elseif length(t) ~= length(S)
	% if S and t are not of the same length, throw an error
    error('t and S must be of identical length!');    
end
% check the level input
if nargin < 3
	% set standard value 0, if level is not given
    level = 0;
end
% check interpolation method input
if nargin < 4
    imeth = 'linear';
end
% make row vectors
t = t(:)';
S = S(:)';
% always search for zeros. So if we want the crossing of 
% any other threshold value "level", we subtract it from
% the values and search for zeros.
S   = S - level;
% first look for exact zeros
ind0 = find( S == 0 ); 
% then look for zero crossings between data points
S1 = S(1:end-1) .* S(2:end);
ind1 = find( S1 < 0 );
% bring exact zeros and "in-between" zeros together 
ind = sort([ind0 ind1]);
% and pick the associated time values
t0 = t(ind); 
s0 = S(ind);
if ~isempty(ind)
    if strcmp(imeth,'linear')
        % linear interpolation of crossing
        for ii=1:length(t0)
            %if abs(S(ind(ii))) >= eps(S(ind(ii)))    % MATLAB V7 et +
            if abs(S(ind(ii))) >= eps*abs(S(ind(ii)))    % MATLAB V6 et -    EPS * ABS(X)
                % interpolate only when data point is not already zero
                NUM = (t(ind(ii)+1) - t(ind(ii)));
                DEN = (S(ind(ii)+1) - S(ind(ii)));
                slope =  NUM / DEN;
                slope_sign(ii) = sign(slope);
                t0(ii) = t0(ii) - S(ind(ii)) * slope;
                s0(ii) = level;
            end
        end
    end
    % extract the positive slope crossing points 
    ind_pos = find(sign(slope_sign)>0);
    t0_pos = t0(ind_pos);
    s0_pos = s0(ind_pos);
    % extract the negative slope crossing points 
    ind_neg = find(sign(slope_sign)<0);
    t0_neg = t0(ind_neg);
    s0_neg = s0(ind_neg);
else
    % empty output
    ind_pos = [];
    t0_pos = [];
    s0_pos = [];
    % extract the negative slope crossing points 
    ind_neg = [];
    t0_neg = [];
    s0_neg = [];
    
end
end

Mathieu NOE 2022-4-1

在 MATLAB Online 中打开

Good news !

maybe I should have added a few comments on my last code modifications :

remove the leading zeros and cutting the signal length automatically
more important :added a smoothing exponential window to the data (has somehow the same effet as smoothdata) - this is because your signal ressembles a lot with a kind of impulse response and that noise and vibration engineers do apply a exp window to impulse response data to make sure they tends to zero as you fill the buffer. This adds a bit of damping to the peaks but avoids the ripples and oscillation. Of course a compromise must be found between too much damping and good looking data. This is driven by the coefficient (gain = 0.01) in the exp function

win = exp(-0.01*(1:numel(t)));

it is aslo assumed that we are focusing on the dominant peak only for the computation of Q

hope this helps understand what I did

VG 2022-4-4

编辑：VG 2022-4-4

Hello,

Thanks for your clear explanination now it make sense!

请先登录，再进行评论。

calculating the x1 and x2 of a peak automatically for every data ?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

7 个评论
显示 5更早的评论隐藏 5更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

calculating the x1 and x2 of a peak automatically for every data ?

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

7 个评论 显示 5更早的评论隐藏 5更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

WeChat

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

7 个评论
显示 5更早的评论隐藏 5更早的评论