3
$\begingroup$

Premise

I have a 3D dataset of a system responding to different frequencies https://pastebin.com/9MyQ6S74

(x: imposed frequency (f1) y: response frequency z: response amplitude)

I know that the system will respond to a frequency that always follows a standard intermodulation product response:

m f1 + n f2, where m,n integer (positive or negative) and f1 the stimulus and f2 a fixed number (here at 600 Hz).

data = Import["SE_example.txt", "TSV"];

ListDensityPlot[data, PlotRange -> Full, ImageSize -> 400, 
     FrameLabel -> {"Stimulus frequency [Hz]", "Response [Hz]"}]

enter image description here

Question

I am interested in finding a way to automatically detect these lines and be able to at least extract slope and offset.

enter image description here

What I have tried so far

The manual method

The most labour intensive method is to simply right click onto the plot and get the coordinate positions of two points along a line, and from there extract the coefficient terms. This is fair enough when the system is simple, but when you have 50+ 'lines' and multiple datasets, it can be a drag.

Using ImageLines

Another thing I tried was to raster the plot into an image and then use ImageLines to give me an idea as to where to look.

(* I unitise the data to make it a bit more clear for the image *)
data[[;; , 3]] = data[[;; , 3]] /. x_ /; x < 0.02 -> 0;
data[[;; , 3]] = data[[;; , 3]] /. x_ /; x > 0 -> 1;


(* then I raster the output of the ListDensityPlot without any frame to not confuse ImageLines *)
rasterImage = 
 Rasterize[
  ListDensityPlot[data, PlotRange -> Full, ImageSize -> 400, 
   Frame -> False, ImagePadding -> None, FrameTicks -> None, 
   PlotRangePadding -> None]]

Manipulate[
 lines = ImageLines[rasterImage, t, d, 
   Method -> {"Segmented" -> True}, MaxFeatures -> 100];
 HighlightImage[rasterImage, lines], {t, 0, 1}, {{d, 0.05}, 0, 1}, 
 Paneled -> False]

Using Manipulate, I played around with different thresholds to see what could work. The RANSAC and Hough methods are not particularly good as the data are not strictly speaking 'continuous' but there are some small gaps for a given line. As you can see below, it does identify the major points when the Method is set to Method->{"Segmented"->True} but because of the segmentation it is generating lines in all sorts of directions which is not helpful.

enter image description here

GradientFilter?

Having a look at this thread, I can appreciate that GradientFilter could be a useful approach, but I am not certain how I could use it in this 3D context.

A potential start could be to reduce the dimensionality based on a threshold:

ListPlot[rdata = Select[data, #[[3]] > 0.02 &][[;; , {1, 2}]], 
 PlotStyle -> AbsoluteThickness[2]]

enter image description here

Any suggestions would be very welcome.

Thank you in advance!

Improve this question
$\endgroup$

1 Answer 1

Reset to default
1
$\begingroup$

I have a slightly inelegant solution, but I will post it here hoping another person will yield a better one.

The answer is partially a brute-force method. It will pick a single datapoint based on their x-axis and then check all the points that form a straight line with a slope of "+1".

It will then repeat for all points and all slopes and at the end extract the unique slopes and offsets that come with such sets.

rdata = Select[data, #[[3]] > 0.03 &][[;; , {1, 2}]] (* threshold the data and reduce dimension *)

Clear[reference, meanSeparation]
reference[pointA_, pointB_, slope_] := (pointA[[2]] - pointB[[2]])/(
   pointA[[1]] - pointB[[1]]) == 
   slope (*&& Abs[pointA[[1]]-pointB[[1]]]<=100*); (* checks the \
slope between two points *)

meanSeparation[list_] := 
  Round[ Mean[
    Table[list[[i + 1, 1]] - list[[i, 1]], {i, 1, 
      Length@list - 
       1}]]]; (* checks the average separation between points*)

uniqueFreqs = 
  Union@rdata[[;; , 
     1]]; (* find the unique frequency points in the 2D-thresholded \
data *) 


dpSets = Table[
  Table[
   Select[
    DeleteCases[ (* delete empty sets*)
     Quiet@
      Table[
       Select[rdata, reference[#, point, slope] &], {point, 
        Select[rdata, #[[1]] == uFreq &]}], {}],(* 
    starting with points which contain one of the unique frequency \
points, select those who maintain a slope of 1 irrespective of offset \
*)
    Length[#] >= 5 && meanSeparation[#] <= 20 &], (* 
   only select sets which contain at least 5 elements and whose mean \
point separation does not exceed step size x2 *)
   {uFreq, uniqueFreqs}], (* repeat for all unique frequencies *)
  {slope, -3, 3, 1}]; (* repeat for all slopes between -3 and 3 *)
Show[ListPlot[rdata, PlotStyle -> AbsoluteThickness[2]], 
 ListPlot[Flatten[dpSets, 2], PlotStyle -> AbsolutePointSize[6]]]

enter image description here

It can now use a basic fitting function to find the slopes and offsets:

Union@(Round /@ ({a, b} /. 
     Chop[FindFit[#, a x + b, {a, b}, x] & /@ Flatten[dpSets, 2]]))
Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.