find multiples of real number in sampled array

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!

Answer 1

[Mathieu NOE]

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

[Tommaso Seresini]

Thank you so much for sharing the code!

The connection between the code snippet and the final goal is that i know the phase (phi) and frequency (f) of the sine of interest. It seems logical to find the odd-multples of pi/2 in the vector (2*pi*f*t +phi) as they are where the antinodes happen.

Findpeaks would struggle in my case, because of a large background that is present and I can't really filter out.

[Mathieu NOE]

hello again

my pleasure

why can't you filter out the background ? have you tried with smoothdata ?

[Tommaso Seresini]

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]

hello

glad to see you have a solution that works !

have a good day

Mathieu