for loop jump an element of an array

Hello there,

it's my first post here, I'm writing a code to perform linear interpolation of a sinusoid and then a resampling with an higeher sampling frequency but i'm facing some trouble since the foor loop seems to jump one value of the sampling time vector and that's a problem when i'm walking through the array to find the overlapping times of the two sampling time vectors.

f0 = 10; % Hz
T0 = 1/f0; % period of sin [s]
fc = 100; %Hz
T = 1/fc; % first sampling period [s]
N = floor(10*(fc/f0)); % number of samples in 10 periods

samples_time = 0:T:(N-1)*T;
samples_value = zeros(1,N);

% first sampling
for i=1:N
    samples_value(i) = sin(2*pi*f0*samples_time(i));
end

fc2 = 200; % Hz
T2 = 1/fc2; % second sampling period [s]
N2 = floor((N-1)*(T/T2)+1);    % number of samples
samples_time2 = 0:T2:(N2-1)*T2;

for j=1:N2
    if samples_time2(j) == 0.64
        disp("found");
    end
    %disp(samples_time2(j));
end

The foor loop seems not to find the value 0.64 in the array, it works with other values. I cannot find a solution, it's very strange. (I'm running R2023b)

Hope someone can help me, thanks

[Stephen23]

"I cannot find a solution, it's very strange."

It is not strange at all: you made the mistake of testing for exact equality of binary floating point numbers.

The solution is also very simple: always compare the absolute difference against a tolerance:

abs(A-B)<tol

More information on this topic:

https://www.mathworks.com/matlabcentral/answers/57444-faq-why-is-0-3-0-2-0-1-not-equal-to-zero

https://www.mathworks.com/matlabcentral/answers/316889-why-two-equal-numbers-are-not-equal

http://www.mathworks.com/matlabcentral/answers/102419-how-do-i-determine-if-the-error-in-my-answer-is-the-result-of-round-off-error-or-a-bug

http://www.mathworks.com/matlabcentral/answers/140656-about-numerical-difficulty-of-equality-comparison-operation-between-two-double-real-numbers-in-matlab

This is worth reading as well:

https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html

Answer 1

[Voss]

From the documentation for the colon operator:

"x = j:i:k creates a regularly-spaced vector x using i as the increment between elements. The vector elements are roughly equal to [j,j+i,j+2*i,...,j+m*i] where m = fix((k-j)/i). However, if i is not an integer, then floating point arithmetic plays a role in determining whether colon includes the endpoint k in the vector, since k might not be exactly equal to j+m*i."

I'd expand that to add, if i is not an integer, then floating point arithmetic plays a role in determining whether colon includes any intended value beyond the first.

To avoid this problem, construct your time vectors using a spacing of 1 in colon, i.e.:

samples_time = (0:N-1)*T;

samples_time2 = (0:N2-1)*T2

instead of:

samples_time = 0:T:(N-1)*T;

samples_time2 = 0:T2:(N2-1)*T2;

(Yes, the same problem of missing values happens in samples_time as well.)

Some illustrations:

f0 = 10; % Hz
T0 = 1/f0; % period of sin [s]
fc = 100; %Hz
T = 1/fc; % first sampling period [s]
N = floor(10*(fc/f0)); % number of samples in 10 periods

Constructing samples_time both ways:

samples_time = 0:T:(N-1)*T

samples_time = 1×100

0 0.0100 0.0200 0.0300 0.0400 0.0500 0.0600 0.0700 0.0800 0.0900 0.1000 0.1100 0.1200 0.1300 0.1400 0.1500 0.1600 0.1700 0.1800 0.1900 0.2000 0.2100 0.2200 0.2300 0.2400 0.2500 0.2600 0.2700 0.2800 0.2900

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

samples_time_test = (0:N-1)*T

samples_time_test = 1×100

0 0.0100 0.0200 0.0300 0.0400 0.0500 0.0600 0.0700 0.0800 0.0900 0.1000 0.1100 0.1200 0.1300 0.1400 0.1500 0.1600 0.1700 0.1800 0.1900 0.2000 0.2100 0.2200 0.2300 0.2400 0.2500 0.2600 0.2700 0.2800 0.2900

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

They look the same, but they are not:

isequal(samples_time,samples_time_test)

ans = logical
   0

In particular, 0.64 is in the one constructed using a colon spacing of 1 but not the one using T spacing:

find(samples_time == 0.64)

ans =

  1×0 empty double row vector

find(samples_time_test == 0.64)

ans = 65

There is an element in samples_time that is very close to 0.64 (about 1e-16 away), but it's not exactly 0.64:

[d,idx] = min(abs(samples_time-0.64))

d = 1.1102e-16

idx = 65

whereas in samples_time_test that element is exactly 0.64:

[d,idx] = min(abs(samples_time_test-0.64))

d = 0

idx = 65

The same is true for the second set of times.

fc2 = 200; % Hz
T2 = 1/fc2; % second sampling period [s]
N2 = floor((N-1)*(T/T2)+1);    % number of samples
samples_time2 = 0:T2:(N2-1)*T2

samples_time2 = 1×199

0 0.0050 0.0100 0.0150 0.0200 0.0250 0.0300 0.0350 0.0400 0.0450 0.0500 0.0550 0.0600 0.0650 0.0700 0.0750 0.0800 0.0850 0.0900 0.0950 0.1000 0.1050 0.1100 0.1150 0.1200 0.1250 0.1300 0.1350 0.1400 0.1450

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

samples_time2_test = (0:N2-1)*T2

samples_time2_test = 1×199

0 0.0050 0.0100 0.0150 0.0200 0.0250 0.0300 0.0350 0.0400 0.0450 0.0500 0.0550 0.0600 0.0650 0.0700 0.0750 0.0800 0.0850 0.0900 0.0950 0.1000 0.1050 0.1100 0.1150 0.1200 0.1250 0.1300 0.1350 0.1400 0.1450

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

isequal(samples_time2,samples_time2_test)

ans = logical
   0

find(samples_time2 == 0.64)

ans =

  1×0 empty double row vector

find(samples_time2_test == 0.64)

ans = 129

And 0.64 isn't the only missing value either, in either vector:

% elements in samples_time_test that aren't in samples_time:
samples_time_test(~ismember(samples_time_test,samples_time))

ans = 1×13

0.6400 0.6500 0.6600 0.6700 0.6800 0.6900 0.7000 0.8200 0.8300 0.9200 0.9300 0.9400 0.9500

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

% elements in samples_time2_test that aren't in samples_time2:
samples_time2_test(~ismember(samples_time2_test,samples_time2))

ans = 1×28

0.5150 0.5250 0.5350 0.5450 0.5550 0.5650 0.5750 0.6400 0.6500 0.6600 0.6700 0.6800 0.6900 0.7000 0.7850 0.7950 0.8050 0.8150 0.8200 0.8250 0.8300 0.9200 0.9300 0.9350 0.9400 0.9450 0.9500 0.9550

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Davide]

Thanks for the clearly reply, but, even if I declare the vector as you showed I can get some strange behaviour by comparing the time samples, am I right?

I mean, the code basically sample a sinusoid with T as sampling period, then perform a resampling with T2 as sampling period by applying a linear interpolator to calculate the new samples. Since a linear interpolator calculate the value in, let's say, T2 by evaluating in that time the value of the equation of the straight line passing through 0 and T, I need to use a for loop and check for each time sample samples_time2(j) if it's between let's say samples_time(start_index) and samples_time(stop_index), so even using >= and <= can I have the same problem?

% samples_time array
% 0                 T                2T
% |-----------------|-----------------|--
%
% samples_time2 array 
% 0     T2   2T2   3T2   4T2   5T2   6T2
% |-----|-----|-----|-----|-----|-----|--

[Voss]

You're welcome!

I'm not sure what you mean. There is no interpolation going on.

Can you provide some simple code that illustrates the strange behavior you observe?

[Stephen23]

"...can I have the same problem?"

Yes. Using two different methods can result in two different values. That is nature of binary floating point numbers.

[Davide]

Sure, of course. As soon as I get home I will provide you the complete code

Answer 2

[the cyclist]

You should not use strict equality to compare floating point numbers.

Instead, check equality with a small tolerance, e.g.

if abs(samples_time2(j) - 0.64) < 1.e-6