[Math] Finding where plots may cross with octave / matlab

Question

I have several data points that are plotted below and I would like to find the frequency value when the amplitude value crosses 4. I've included an example note: this is an example and the amplitude will vary along with the data points in the example below. I've circled the answer graphically but I'm not sure how to compute it mathematically and get all the values for the frequencies I desire. How can I do this with octave / matlab? Also is there a mathematical term for what I'm trying to do?

In this example I'm trying to get 5 frequencies (but this is just an example) I know two answers are 30 and 80 but not sure how to get the rest. The full list could be thousands. I'm using octave 3.8.1

clear all,clf, clc,tic
%graphics_toolkit gnuplot %use this for now it's older but allows zoom
freq=[20,30,40,50,60,70,80];
amp_orig=[2,4,3,7,1,8,4];
amp_inv=[6,4,5,1,7,0,4];


plot(freq,amp_orig,'-bo')
hold on
plot(freq,amp_inv,'-r*')
xlabel ("Frequency");
ylabel ("Amplitude");

PS: My data will not always oscillate when the amplitude is equal to 4 nor will the data always be a mirror this is a simple example to show what i'm trying to accomplish.
Thanks

Answer 1

(quick comment regarding @Hayoung.Chung 's answer since I can't comment on it yet. While that code might prove to be a valid solution for Matlab, but he mentioned Octave and Matlab. Use of Matlab File Exchange code in Octave appears to be against The Mathworks's restrictive Terms of Use for the File Exchange. Specifically Section 2.a.iii, even though it appears to conflict with Section 2.c and Section 5. Octave has been avoiding using File Exchange code internally for that reason. Obviously, you are in control of your own actions in that regard.)

what you're looking for is a "zero" crossing of the interpolated line between points after you've shifted it by -4. if your data had many more data points between zero crossings you could do a polyline fit and use the matlab/octave tools to find zeros of a polyline. Here, you may just want to make a loop that generates the equation of a line between successive points, does a (shifted) zero-find with that line, and then has some clean up functions to account for points that hit y=0 exactly (hence both segments will report the same zero location giving you duplicate points), or don't hit it at all. It would then be simple to loop this through all of your collected points.

As your comments indicated that 'always crossing 4' was arbitrary, and wouldn't always even cross at a fixed value, the script below will take care of the situation your provided in the your question, as well as the case where they randomly cross. it takes the two curves, subtracts them, and finds where that equals zero. it then returns the intersection frequencies and the amplitudes at those points.

note that the code below could be vectorized to run much faster without a for loop, but i think the for loop is more instructive at this point. what I'm doing is really repeated linear interpolation, and there are built in functions to do that for an array of points. I believe even the direct method here could likely be vectorized as well. But this appears sufficient for your purpose.

clear all;
freq=[20,30,40,50,60,70,80];

#amp_orig1 = floor(10*rand(1,7))
#amp_orig2 = floor(10*rand(1,7))

amp_orig1 = [8 8 9 0 5 1 2]
amp_orig2 = [8 8 4 2 4 5 6]


numpoints = numel(freq);
crossings_xvals = NaN(numpoints-1,2); #ones with no crossing will stay as NaN
crossings_yvals = NaN(numpoints-1,2);

#crossing_amp_value = 4; #selected since your inversion is about amp=4.

amp_shift = amp_orig1-amp_orig2; #solve for zero crossing in rest of code

#first find the zero crossing for each segment
for fval = 1 : numpoints-1

    xs = freq(fval:fval+1);
    ys = amp_shift(fval:fval+1);

    slope = (ys(2)-ys(1))/(xs(2)-xs(1));
    yintercept = ys(1) - slope*xs(1);

    xintercept = -yintercept/slope;

    #exclude cases where it doesn't intersect in that segment
    if (xintercept>=xs(1) && xintercept<=xs(2))
      #get actual y-value at crossing
      ys_orig = amp_orig1(fval:fval+1);
      yval_at_crossing= ((ys_orig(2)-ys_orig(1))/(xs(2)-xs(1)))*(xintercept-xs(1)) + ys_orig(1);

      crossings_xvals(fval) = xintercept;
      crossings_yvals(fval) = yval_at_crossing;
    end

end

#drop NaNs to account for it not always oscillating/crossing amp=4;
crossings_xvals = crossings_xvals(~isnan(crossings_xvals));
crossings_yvals = crossings_yvals(~isnan(crossings_yvals));

#remove duplicates. floating point arithmetic, so check if they're closer than 
#numerical error eps
dups = find(diff(crossings_xvals)<=eps);

#final xvalues where it crosses/equals '4'
crossings_xvals(dups)=[];
crossings_yvals(dups)=[];
crossings = [crossings_xvals,crossings_yvals]

plot(freq,amp_orig1,freq,amp_orig2,freq,amp_shift,crossings_xvals,crossings_yvals,'ok',[freq(1),freq(end)],[0,0],'--');
xlabel('frequency');ylabel('amplitude');
legend('orig1','orig2','orig_{diff}','intersect');

if they're sampled at different frequency values, you'll have to do a lot more interpolating freq values before calculating the rest.