1
$\begingroup$

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");

enter image description here

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

$\endgroup$
  • $\begingroup$ For any two points, you can define a line. Then all you have to do is solve for when that line hits $\mathrm{Amplitude} = 4$. $\endgroup$ – AnonSubmitter85 Sep 22 '14 at 23:16
  • $\begingroup$ just do (abs(conv(((amp_inv-amp_orig)>0)-0.5,[1,-1],'valid'))>0.999) the positions with 1 means "crossing in between these frequencies". If you need higher resolution than "yes somewhere in between these frequencies" you can just do an interpolation before you do that line. $\endgroup$ – mathreadler May 23 '17 at 9:21
1
$\begingroup$

(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.

sample output from script

$\endgroup$
  • $\begingroup$ No my data will not always oscillate about amplitude=4? $\endgroup$ – Rick T Aug 5 '15 at 4:15
  • $\begingroup$ Ok, your amp_orig and amp_inv just happen to invert about the line amp=4 then? will they always be mirrors of each other like in the example you show, but possibly about a number different than 4? or could they be very distinct curves you're trying to find crossings for? the approach is fairly simple if they are always mirrored about the constant value you're looking for. $\endgroup$ – Nick J Aug 5 '15 at 18:38
  • 1
    $\begingroup$ how you get the y-value depends on a couple things. 1st: your initial problem said "find where they cross 4". do you expect that they'll always be crossing at the same value whenever they cross, such that we can always look for it to cross that value? or do you expect the value might be different for each crossing, and you need to account for that? $\endgroup$ – Nick J Aug 7 '15 at 15:53
  • 1
    $\begingroup$ actually, looking back on the code I think my suggestion after the code, to change the line setting amp_shift, will solve for them crossing anywhere, not just at a fixed value '4'. the subtraction will basically make a new line segments that = zero where the two lines have the same value. so, the rest of the code (solving for the zero crossing) should still work. $\endgroup$ – Nick J Aug 7 '15 at 15:58
  • 1
    $\begingroup$ ok, last comment... I updated the code to put the subtraction in, and then to also crab the yvals at the intersection. it will now work for any overlap. note that when the curves overlap on a segment (example data between first two frequency values) it is undefined because there are an infinite number of over lap points. the next segment will put a point where they separate (like at freq=30). two commented out lines generate random data to test the code, use them instead of the following lines setting amp_orig1 and 2, and you should be able to run the code a few times and see the output. $\endgroup$ – Nick J Aug 7 '15 at 16:50
1
$\begingroup$

If your troble is only about extracting intersecting point and not on mathematical formulation, how about try "Curve intersection" function in Matlab File Exchange? I've been using it for postprocessing for years, and it has been worked fine (for thousands of data)

try visit http://www.mathworks.com/matlabcentral/fileexchange/22441-curve-intersections/content/InterX.m

$\endgroup$
0
$\begingroup$

The code provided worked will with Octave. I tried to figure out how to display the coordinates of the intersections for a few hours and didn't find anything. Any suggestions.

clear all;
xaxis=[-10:10];

#points_1 = floor(10*rand(1,7))
#points_2 = floor(10*rand(1,7))

points_1 = xaxis.^2;
points_2 = 4*xaxis -4;



numpoints = numel(xaxis);
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 = points_1-points_2; #solve for zero crossing in rest of code

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

    xs = xaxis(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 = points_1(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(xaxis,points_1,xaxis,points_2,xaxis,amp_shift,crossings_xvals,crossings_yvals,'ok',[xaxis(1),xaxis(end)],[0,0],'--');
xlabel('xaxis');ylabel('yaxis');

legend('orig1','orig2','orig_{diff}','intersect');
$\endgroup$
  • 2
    $\begingroup$ You might want to format your post a little. $\endgroup$ – StubbornAtom May 23 '17 at 10:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.