find multiples of real number in sampled array

2 次查看(过去 30 天)
Hi,
I have an array t with a serie of values with equal spacing dt. I would like to find the values of t which are multiples of a given real number b, for different values of a parameter a.
Since i am dealing with real numbers, I don't expect to find the exact multiples. So I came up with the following solution:
dt=1e-4;
t=0:dt:1;
b=1;
a=24;
multiples=t(mod(a*t,b)<dt)
I find that this solution works well for some combination of numbers but for others, such as those in the example, I have some rounding errors and i miss some points. Am I doing something wrong?
The ultimate goal would be to have a function to find the nodes and antinodes of a sine with an arbitrary frequency and phase. this method works well as long as the number of periods in t is small (2,3) but fails and is inconsistend as the frequency grows.
Thank you so much for your help!

回答(1 个)

Mathieu NOE
Mathieu NOE 2022-4-8
hello
I had some code laying around and found it could be used for your purpose (the sine wave problem , and btw i could not really make the connection with the posted code).
so demo with an exponential decaying sine wave - you can try with your own signal
  • nodes = points of zero crossing (here shown both positive slope and negative slope zero crossing points)
  • antinodes = local maxima (here I just picked the positive ones)
clc
clearvars
n=1000;
x=linspace(0,2*pi,n);
y = sin(3*x-1).*exp(-x/10);
% antinodes = find positive local maxima
[tf, P] = islocalmax(y,'MinProminence',max(y)/3);
x_antinodes = x(tf);
y_antinodes = y(tf);
% nodes = find points crossing the zero y value
threshold = 0; % your value here
[t0_pos,s0_pos,t0_neg,s0_neg]= crossing_V7(y,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"
figure(1)
plot(x,y,x_antinodes,y_antinodes,'dk',x,threshold*ones(size(x)),'k--',t0_pos,s0_pos,'dr',t0_neg,s0_neg,'dg','linewidth',2,'markersize',12);grid on
legend('signal','antinodes','nodes level threshold','positive slope nodes','negative slope nodes');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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
  4 个评论
Tommaso Seresini
Tommaso Seresini 2022-4-12
thank you for the suggestion. I don't want to filter out the background because it is not noise and it is relevant for the analysis I need to make.
One possible solution could be: isolate the sine of interest with bandpass -> find antinodes with findpeaks -> use the indexes in the original signal.
So far, I am using another workaround which uses only the fact that I know phase and frequency of the desired signal.
Instead of looking for when the remainder "mod(omega*t+phase,pi/2)" reaches a treshold it looks when its derivative changes sign. To distinguish antinodes (odd multiples of pi/2) and nodes (even multiples of pi/2), take the set difference between the set of multiples of pi and of pi/2. it goes as follows:
function [nodes_idx] = FindAntinodes(t,omega,phase)
%FindAntinodes find the nodes of a sine with angular frequency omega and a given phase
rest_90=mod(omega*t+phase,pi/2);
rest_180=mod(omega*t+phase,pi);
sgn_90=sign(diff(rest_90));
sgn_180=sign(diff(rest_180));
nodes_idx_90=find(sgn_90==-1);
nodes_idx_180=find(sgn_180==-1);
nodes_idx=setdiff(nodes_idx_90,nodes_idx_180)+1;
% figure(1)
% hold on
% plot(t,sin(omega*t+phase))
% plot(t,rest_90/pi*2);
% plot(t(2:end),sgn_90);
% plot(t(2:end),sgn_180);
% plot(t,rest_180/pi);
%
% plot(t(nodes_idx),sin(omega*t(nodes_idx)+phase),'o')
end
Mathieu NOE
Mathieu NOE 2022-4-12
hello
glad to see you have a solution that works !
have a good day
Mathieu

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Multirate Signal Processing 的更多信息

产品


版本

R2019a

Community Treasure Hunt

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

Start Hunting!

Translated by