I am trying to write a function to calculate frequency of a signal by comparing values to a threshold.

12 views (last 30 days)
I am trying to create a function which calculates the frequency of a given signal by comparing it to a threshold. In the function, x is the threshold, y is the hysteresis window length, and z is the number of cycles to run the function to find an average frequency. There is an array (arrayWithLatestValues) of size 500 (arraySize) which contains the signal values and reading backward through the array, if the current value is less than the threshold and the next value is greater than the threshold, this is considered a "falling edge". I am trying to calculate the index of this falling edge and check of the hysteresis condition holds true, i.e., that the subsequent values after this falling edge stay above the threshold. I check this condition by checking the next y values after the initial falling edge and making sure they stay above this threshold. I am trying to repeat this to calculate another falling edge, and the difference between these indices can then allow me to calculate the signal period and frequency.
Is this a good way to go about achieving my outlined functionality? What changes can I make to the code logic to make it work more effectively?
I am currently getting this error message when running this function in its current state:
"Not enough input arguments.
Error in fn_detectFrequency (line 12)
for c = 0:z
"
function frequency = fn_detectFrequency(x,y,z)
global arraySize;
global arrayWithLatestValues;
global samplingRate;
v = 0;
w = 0;
hysteresisWorked1 = 0;
hysteresisWorked2 = 0;
indexDifference = 0;
period = 0;
for c = 0:z
for i = arraySize:-1:1
currentValue = arrayWithLatestValues(1,i);
nextValue = arrayWithLatestValues(1,i-1);
if (currentValue < x) && (nextValue > x)
v = i;
for d = v:-1:v-y
nextValue1 = arrayWithLatestValues(1,v-1);
if (nextValue1>x)
hysteresisWorked1 = 1;
else
hysteresisWorked1 = 0;
end
end
end
currentValue1 = arrayWithLatestValues(1,v-y);
nextValue2 = arrayWithLatestValues(1,v-y-1);
if (currentValue1 < x) && (nextValue2 > x)
w = i;
for e = w:-1:w-y
nextValue3 = arrayWithLatestValues(1,w-1);
if(nextValue3>x)
hysteresisWorked2 = 1;
else
hysteresisWorked2 = 0;
end
end
end
if (hysteresisWorked1) && (hysteresisWorked2)
indexDifference = v - w;
period = indexDifference/samplingRate;
frequency = 1/period;
end
end
end
end
  3 Comments
Richard Liu
Richard Liu on 26 Jun 2022
Edited: Richard Liu on 26 Jun 2022
The sign function looks promising, but is there a way to do the comparison around a value other than 0?
We cannot use FFT or PSD for this problem because the goal is to crunch out the frequency by finding the time difference between these falling edges, essentially trying to find the "period" of the function since we assume that the time between falling edges should be the same for a periodic signal.
dpb
dpb on 26 Jun 2022
Edited: dpb on 26 Jun 2022
" but is there a way to do the comparison around a value other than 0?"
As written it's about the threshold value. sign is applied to the difference, not the raw signal.
Also note that the hysterisis thingie can be implemented with diff, too, again there's no loop needed to find where the crossings are a minimum distance apart.

Sign in to comment.

Answers (1)

Mathieu NOE
Mathieu NOE on 27 Jun 2022
hello
FYI
demo with constant and non constant signals
period are computed by getting the threshold crossing time data
hope it helps
clc
clearvars
% dummy data
n=1000;
x= 10*(0:n-1)/n;
y1 = sin(6*x -0.5);
y2 = sin(x.^2 -0.5);
threshold = max(y1)/2; % half amplitude
[t0_pos1,s0_pos1,t0_neg1,s0_neg1]= crossing_V7(y1,x,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"
[t0_pos2,s0_pos2,t0_neg2,s0_neg2]= crossing_V7(y2,x,threshold,'linear'); % positive (pos) and negative (neg) slope crossing points
% periods
d1 = diff(t0_pos1)
period1 = [d1(1) d1];
d2 = diff(t0_pos2)
period2 = [d2(1) d2];
figure(1)
subplot(3,1,1),plot(x,y1,'b',t0_pos1,s0_pos1,'dr',t0_neg1,s0_neg1,'dg','linewidth',2,'markersize',12);grid on
legend('signal','signal positive slope crossing points','signal negative slope crossing points');
subplot(3,1,2),plot(x,y2,'b',t0_pos2,s0_pos2,'dr',t0_neg2,s0_neg2,'dg','linewidth',2,'markersize',12);grid on
legend('signal','signal positive slope crossing points','signal negative slope crossing points');
subplot(3,1,3),plot(t0_pos1,period1,t0_pos2,period2,'linewidth',2,'markersize',12);grid on
legend('signal 1 period','signal 2 period');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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

Products


Release

R2022a

Community Treasure Hunt

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

Start Hunting!