Using interpolation to find the zeros of a data set that forms a sinusoidal profile

Hello Matlab friends,

I ran into a roadblock when trying to find the zero values of a pair of sine waves that are each defined by a discrete set of points.

In the image I uploaded, I highlighted the nature of the problem. My goal is to isolate the two points that are closest to the x-axis (above and below) for each instance where the line crosses the x-axis. Then, I would like to interpolate between each of these two-point sets to find the x-axis values for which y=0. (These x-axis, or zero values should be in an array).

I was able to isolate several index values for points that are above and below a certain y-axis value. Here is the code for how I did this:

"f" defines the set of data for the displacement waveform, and "d" defines the data for the force waveform. "lvdt_nearaxis" defines the index values of points in "f" that are between -2000 and 2000, or close to the x-axis, and so on. "bone" and "fish" define the consecutive points in "lvdt_nearaxis" and "force_nearaxis". The two-point-sets between which I would like to interpolate will be hidden somewhere in "bone" and "fish".

lvdt_nearaxis = find( f >= -2000 & f <= 2000 );
force_nearaxis = find( d >= -2000 & d<= 2000 ); 

dog = find(diff( lvdt_nearaxis ) ==1 );
bark = [ dog; dog + 1 ];
bone = lvdt_nearaxis( bark  );

cat = find(diff( force_nearaxis ) ==1 );
meow = [ cat; cat + 1 ];
fish = force_nearaxis( meow );

I am not quite sure how to proceed from here. I would appreciate any help! :)

To further elucidate what I am trying to do, I gathered both data sets from a mechanical spectrometer I built. One sine wave is generated by a displacement sensor that measures the motions of a plate that cyclically compresses a sample, and the other waveform is from a force sensor that measures the transmitted force.

Answer 1

[Star Strider]

Try this code with your data:

x = linspace(0,10*pi,1E+5);                                             % Create Data
y = sin(1./(x+0.01));                                                   % Create Data

zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0);                    % Returns Approximate Zero-Crossing Indices Of Argument Vector

dy = zci(y);                                                            % Indices of Approximate Zero-Crossings
for k1 = 1:size(dy,1)-1
    b = [[1;1] [x(dy(k1)); x(dy(k1)+1)]]\[y(dy(k1)); y(dy(k1)+1)];      % Linear Fit Near Zero-Crossings
    x0(k1) = -b(1)/b(2);                                                % Interpolate ‘Exact’ Zero Crossing
    mb(:,k1) = b;                                                       % Store Parameter Estimates (Optional)
end
freq = 1./(2*diff(x0));                                                 % Calculate Frequencies By Cycle

figure(1)
plot(x, y)
hold on
plot(x0, zeros(size(x0)), '+r')
hold off
grid
axis([0  1    ylim])

You may have to tweak it slightly to use your variables, or you can create a function from it.

To use it as a function, first save this function in a separate file on your MATLAB user path as ‘data_zeros.m’:

      function [x0,freq] = data_zeros(x,y)
      zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0);                    % Returns Approximate Zero-Crossing Indices Of Argument Vector
      dy = zci(y);                                                            % Indices of Approximate Zero-Crossings
      for k1 = 1:size(dy,1)-1
          b = [[1;1] [x(dy(k1)); x(dy(k1)+1)]]\[y(dy(k1)); y(dy(k1)+1)];      % Linear Fit Near Zero-Crossings
          x0(k1) = -b(1)/b(2);                                                % Interpolate ‘Exact’ Zero Crossing
          mb(:,k1) = b;                                                       % Store Parameter Estimates (Optional)
      end
      freq = 1./(2*diff(x0));                                                 % Calculate Frequencies By Cycle
      end

---------------------------------------------------------------------

The function call code would be:

x = linspace(0,10*pi,1E+5);                                             % Create Data
y = sin(1./(x+0.01));                                                   % Create Data

[xz,frq] = data_zeros(x,y);

figure(1)
plot(x, y)
hold on
plot(xz, zeros(size(xz)), '+r')
hold off
grid
axis([0  1    ylim])

The ‘freq’ variable gives the half-cycle frequencies, because I needed that when I wrote the code. If you don’t need it, don’t return it in the function call.