7
$\begingroup$

This is a sample data

data = {{-6., 0.}, {-1.85, 0.}, {-1.84195, 13.04}, {-1.80547, 
212.84}, {-1.73252, 46.4}, {-1.69605, 16.025}, {-1.65957, 
35.36}, {-1.6231, 45.265}, {-1.58663, 7.58}, {-1.55015, 
45.96}, {-1.44073, 69.41}, {-1.40426, 8.2}, {-1.07599, 
8.84}, {-1.00304, 8.97}, {-0.930091, 9.11}, {-0.893617, 
15.05}, {-0.857143, 17.38}, {-0.820669, 9.745}, {-0.784195, 
2.045}, {-0.74772, 22.4267}, {-0.711246, 2.44}, {-0.674772, 
6.075}, {-0.601824, 6.24}, {-0.56535, 19.67}, {-0.528875, 
10.01}, {-0.492401, 29.565}, {-0.455927, 10.24}, {-0.419453, 
10.37}, {-0.382979, 31.514}, {-0.346505, 15.682}, {-0.31003, 
6.995}, {-0.273556, 7.105}, {-0.237082, 14.192}, {-0.200608, 
24.2917}, {-0.164134, 30.242}, {-0.12766, 12.268}, {-0.0911854, 
10.232}, {-0.0547113, 6.2}, {-0.0182371, 25.63}, {0.0182371, 
29.188}, {0.0547113, 8.1675}, {0.0911854, 7.06667}, {0.12766, 
7.16}, {0.164134, 28.7383}, {0.200608, 18.3625}, {0.237082, 
8.585}, {0.273556, 8.415}, {0.31003, 15.9475}, {0.346505, 
23.682}, {0.382979, 12.0267}, {0.419453, 29.295}, {0.455927, 
16.36}, {0.492401, 9.13}, {0.528875, 22.9867}, {0.56535, 
26.43}, {0.601824, 51.0913}, {0.638298, 27.8514}, {0.674772, 
29.108}, {0.711246, 16.8533}, {0.74772, 15.112}, {0.784195, 
31.7615}, {0.820669, 23.5644}, {0.857143, 39.652}, {0.893617, 
39.1586}, {0.930091, 34.349}, {0.966565, 25.388}, {1.00304, 
28.89}, {1.03951, 29.1918}, {1.07599, 42.6718}, {1.11246, 
28.6156}, {1.14894, 60.143}, {1.18541, 58.38}, {1.22188, 
53.6367}, {1.25836, 21.2467}, {1.29483, 23.8542}, {1.33131, 
19.1442}, {1.36778, 34.9167}, {1.40426, 21.45}, {1.44073, 
27.0567}, {1.4772, 41.0813}, {1.51368, 41.2724}, {1.55015, 
29.635}, {1.58663, 23.2117}, {1.6231, 50.5792}, {1.65957, 
47.3894}, {1.69605, 33.3813}, {1.73252, 39.548}, {1.769, 
37.9856}, {1.80547, 34.9631}, {1.84195, 40.6659}, {1.87842, 
45.345}, {1.91489, 16.9433}, {1.95137, 20.1307}, {1.98784, 
37.7388}, {2.02432, 37.7213}, {2.06079, 22.66}, {2.09726, 
69.326}, {2.13374, 53.78}, {2.17021, 35.7221}, {2.20669, 
25.4244}, {2.24316, 20.7723}, {2.27964, 38.6847}, {2.31611, 
9.05125}, {2.35258, 37.4248}, {2.38906, 42.3067}, {2.42553, 
56.1165}, {2.46201, 32.8183}, {2.49848, 37.9471}, {2.53495, 
33.1789}, {2.57143, 34.9625}, {2.6079, 39.9675}, {2.64438, 
42.0679}, {2.68085, 40.0265}, {2.71733, 37.4567}, {2.7538, 
46.3014}, {2.79027, 29.1707}, {2.82675, 47.1314}, {2.86322, 
49.8656}, {2.8997, 24.9844}, {2.93617, 44.8233}, {2.97264, 
18.9875}, {3.00912, 31.8763}, {3.04559, 63.9644}, {3.08207, 
29.475}, {3.11854, 44.6591}, {3.15502, 32.72}, {3.19149, 
26.3475}, {3.22796, 25.315}, {3.26444, 2.38}, {3.30091, 
5.1}, {3.33739, 7.425}, {3.37386, 29.54}, {3.41033, 
70.9567}, {3.44681, 0.45}, {3.48328, 13.355}, {3.51976, 
2.215}, {3.55623, 8.135}, {3.59271, 19.035}, {3.62918, 
9.895}, {3.66565, 1.17}, {3.70213, 0.39}, {3.7386, 0.39}, {3.77508, 
22.3}, {3.81155, 1.79}, {3.84802, 1.08}, {3.8845, 0.36}, {3.92097, 
3.25}, {3.95745, 22.04}, {4.0304, 0.99}, {4.06687, 0.33}, {4.10334, 
0.33}, {4.24924, 0.91}, {4.28571, 0.3}, {4.32219, 2.72}, {4.46809, 
1.93}, {4.50456, 0.83}, {4.54103, 0.27}, {4.57751, 3.04}, {4.75988, 
1.25}, {4.79635, 0.25}, {4.83283, 0.25}, {5.08815, 1.12}, {5.12462, 
0.22}, {5.16109, 21.83}, {5.45289, 1.}, {5.48936, 0.2}, {5.52584, 
3.}, {5.89058, 0.88}, {5.92705, 0.17}}

and this is the corresponding list plot

L0 = ListPlot[data, Frame -> True, Axes -> False, Joined -> True, 
PlotStyle -> {Black}, AspectRatio -> 1, PlotRange -> {0, 100}]

enter image description here

My question is how can I reduce the "noise" of the data so as to obtain something similar with red line which has been drawn by hand.

Improve this question
$\endgroup$
2
  • $\begingroup$ It really depends on what the data represent and what your requirements are for the result. The Savitzky-Golay filter is a generalized moving average and can be reasonably good for this sort of thing, under the right circumstances. $\endgroup$
    – Oleksandr R.
    Jul 14, 2015 at 1:16
  • $\begingroup$ Have a look at the Gaussian Processes below! $\endgroup$
    – Dr Xorile
    Jul 14, 2015 at 14:32

6 Answers 6

Reset to default
3
$\begingroup$

So it's hard to say from your data what is signal and what is noise. But one way to get what you are looking at is to plot a moving average.

Here you plot the average of three points,

data3pointaverage = 
    Table[{data[[n, 1]], Mean[data[[n - 1 ;; n + 1, 2]]]}, {n, 2,Length@data - 1}];
ListLinePlot[{data, data3pointaverage}, PlotStyle -> {Automatic, Red}]

enter image description here

I think doing a 5-point average gets more to what you want

data5pointaverage=
  Table[{data[[n,1]],Mean[data[[n-2;;n+2,2]]]},{n,3,Length@data-2}];
ListLinePlot[{data,data5pointaverage},PlotStyle->{Automatic,Red}]

enter image description here

Improve this answer
$\endgroup$
12
$\begingroup$

Mathematica has a lot of utilities for smoothing.

If your data is evenly sampled, a simple MovingAverage filter may suit your needs, but your abscissa values jump around a bit:

Part[Differences@data, All, 1] // ListPlot

enter image description here

In this case, you can get a more accurate smoothed curve with MovingMap, which can deal with irregularly spaced data:

GraphicsGrid[{{
ListPlot[{data, MovingAverage[data, 7]},
    Joined -> True, PlotLabel -> "Pointwise Moving Average"],
ListPlot[{data, MovingMap[Mean, data, {0.3, Center}, "Reflected"]},
    Joined -> True, PlotLabel -> "Mapped Moving Average"]}}, 
 ImageSize -> 650]

enter image description here

For another useful approach, you may want to take a look at the answers to this question. Here's @Szabolcs answer applied to your data:

ListPlot[
  Mean[Last /@ #] & /@ 
    GatherBy[data, Ceiling[First[#], 0.25] &],
      Joined -> True]

enter image description here

Just remember the abscissa in that result is in bin units.

Improve this answer
$\endgroup$
8
$\begingroup$

You may want to explore the use of LowpassFilter for smoothing.

fd = Transpose@{data[[All, 1]], LowpassFilter[data[[All, 2]], 0.5]};
Improve this answer
$\endgroup$
4
$\begingroup$

It's a little blocky in the "lump" (-1 to +4 on the x-axis), but you might also like

Show[
    ListPlot[data, 
        Joined -> True, PlotRange -> All, PlotStyle -> Black],
    ListLinePlot[Transpose[{
        data[[All, 1]],
        InverseWaveletTransform[
          WaveletThreshold[
            DiscreteWaveletTransform[data[[All, 2]]]]]
        }], 
        Joined -> True, PlotRange -> All, PlotStyle -> Red]
]

If I had any idea how to post images, I would. The "benefit" over some of the other solutions is that the spike at ~-2 is much better preserved.

The actual computational part is 

        InverseWaveletTransform[
          WaveletThreshold[
            DiscreteWaveletTransform[data[[All, 2]]]]]

the rest just splits out the data into parts and glues them back together for plotting. There is much "automagic" in the version as written, however, there are many opportunities for specializing the thresholding method in WaveletThreshold.

I was led to this method because the hand drawn solution wants to retain features whose widths differ by more than an order of magnitude. Most (simple) kernel based methods have a hard time doing this.

Since I have something to add to the discussion (the above), I now have some complaints about the data...

  • The x-axis values are far from uniformly spaced
  • There appears to be missing data
  • There is no segment of data that can be used to establish a noise model. The segment from -6 to -2 should have contained "noise" so that one could estimate the nature of the noise. Ideally, the baseline behaviour in that region would have varied enough to discriminate between additive and multiplicative noise. In the absence of a noise model, it is impossible to determine whether the spike at ~-3.4 is noise or is signal. Your diagram rejects it, but it's the third most energetic feature in the trace, so it's not so clear that the diagram is correct.
  • The y values have been rectified. This makes most eigensystem fitting/feature methods behave poorly. If at all possible, don't throw half the information in your signal (the signs) away.
Improve this answer
$\endgroup$
1
$\begingroup$

Gaussian Processes are a very useful approach for such things.

I wrote about them here (where it wasn't what the OP wanted), but it's basically a text book example for your situation:

Curve Fitting to Represent Any Data

It's computationally more expensive than a rolling average, but mathematically far superior. Basically a function-free fit. You can also set some meta parameters which essentially set the length scale to distinguish between noise and data.

Here's the approach using python code for this example for different length scales:

enter image description here

enter image description here

enter image description here

What's interesting, is that it does good interpolation. For any value of $x$, you can calculate a corresponding $f^*(x)$ value, based on the matrix (harder to do with rolling average for example). You can also change how closely it sticks to the signal depending on the noise (sigmaf, sigman and l in the code).

This is the (python) code:

from scipy import *
from bisect import *
from pylab import *
import scipy.linalg as la
from scipy import interpolate


data = [[-6., 0.], [-1.85, 0.], [-1.84195, 13.04], [-1.80547, 
212.84], [-1.73252, 46.4], [-1.69605, 16.025], [-1.65957, 
35.36], [-1.6231, 45.265], [-1.58663, 7.58], [-1.55015, 
45.96], [-1.44073, 69.41], [-1.40426, 8.2], [-1.07599, 
8.84], [-1.00304, 8.97], [-0.930091, 9.11], [-0.893617, 
15.05], [-0.857143, 17.38], [-0.820669, 9.745], [-0.784195, 
2.045], [-0.74772, 22.4267], [-0.711246, 2.44], [-0.674772, 
6.075], [-0.601824, 6.24], [-0.56535, 19.67], [-0.528875, 
10.01], [-0.492401, 29.565], [-0.455927, 10.24], [-0.419453, 
10.37], [-0.382979, 31.514], [-0.346505, 15.682], [-0.31003, 
6.995], [-0.273556, 7.105], [-0.237082, 14.192], [-0.200608, 
24.2917], [-0.164134, 30.242], [-0.12766, 12.268], [-0.0911854, 
10.232], [-0.0547113, 6.2], [-0.0182371, 25.63], [0.0182371, 
29.188], [0.0547113, 8.1675], [0.0911854, 7.06667], [0.12766, 
7.16], [0.164134, 28.7383], [0.200608, 18.3625], [0.237082, 
8.585], [0.273556, 8.415], [0.31003, 15.9475], [0.346505, 
23.682], [0.382979, 12.0267], [0.419453, 29.295], [0.455927, 
16.36], [0.492401, 9.13], [0.528875, 22.9867], [0.56535, 
26.43], [0.601824, 51.0913], [0.638298, 27.8514], [0.674772, 
29.108], [0.711246, 16.8533], [0.74772, 15.112], [0.784195, 
31.7615], [0.820669, 23.5644], [0.857143, 39.652], [0.893617, 
39.1586], [0.930091, 34.349], [0.966565, 25.388], [1.00304, 
28.89], [1.03951, 29.1918], [1.07599, 42.6718], [1.11246, 
28.6156], [1.14894, 60.143], [1.18541, 58.38], [1.22188, 
53.6367], [1.25836, 21.2467], [1.29483, 23.8542], [1.33131, 
19.1442], [1.36778, 34.9167], [1.40426, 21.45], [1.44073, 
27.0567], [1.4772, 41.0813], [1.51368, 41.2724], [1.55015, 
29.635], [1.58663, 23.2117], [1.6231, 50.5792], [1.65957, 
47.3894], [1.69605, 33.3813], [1.73252, 39.548], [1.769, 
37.9856], [1.80547, 34.9631], [1.84195, 40.6659], [1.87842, 
45.345], [1.91489, 16.9433], [1.95137, 20.1307], [1.98784, 
37.7388], [2.02432, 37.7213], [2.06079, 22.66], [2.09726, 
69.326], [2.13374, 53.78], [2.17021, 35.7221], [2.20669, 
25.4244], [2.24316, 20.7723], [2.27964, 38.6847], [2.31611, 
9.05125], [2.35258, 37.4248], [2.38906, 42.3067], [2.42553, 
56.1165], [2.46201, 32.8183], [2.49848, 37.9471], [2.53495, 
33.1789], [2.57143, 34.9625], [2.6079, 39.9675], [2.64438, 
42.0679], [2.68085, 40.0265], [2.71733, 37.4567], [2.7538, 
46.3014], [2.79027, 29.1707], [2.82675, 47.1314], [2.86322, 
49.8656], [2.8997, 24.9844], [2.93617, 44.8233], [2.97264, 
18.9875], [3.00912, 31.8763], [3.04559, 63.9644], [3.08207, 
29.475], [3.11854, 44.6591], [3.15502, 32.72], [3.19149, 
26.3475], [3.22796, 25.315], [3.26444, 2.38], [3.30091, 
5.1], [3.33739, 7.425], [3.37386, 29.54], [3.41033, 
70.9567], [3.44681, 0.45], [3.48328, 13.355], [3.51976, 
2.215], [3.55623, 8.135], [3.59271, 19.035], [3.62918, 
9.895], [3.66565, 1.17], [3.70213, 0.39], [3.7386, 0.39], [3.77508, 
22.3], [3.81155, 1.79], [3.84802, 1.08], [3.8845, 0.36], [3.92097, 
3.25], [3.95745, 22.04], [4.0304, 0.99], [4.06687, 0.33], [4.10334, 
0.33], [4.24924, 0.91], [4.28571, 0.3], [4.32219, 2.72], [4.46809, 
1.93], [4.50456, 0.83], [4.54103, 0.27], [4.57751, 3.04], [4.75988, 
1.25], [4.79635, 0.25], [4.83283, 0.25], [5.08815, 1.12], [5.12462, 
0.22], [5.16109, 21.83], [5.45289, 1.], [5.48936, 0.2], [5.52584, 
3.], [5.89058, 0.88], [5.92705, 0.17]]

xdata=[]
ydata=[]
for d in data:
    xdata.append(d[0])
    ydata.append(d[1])
x = linspace(-6,6,400)
xdata = array(xdata)
ydata = array(ydata)


#covariance function
sigmaf=1   #sigma_f  (noise in f)
sigman=1      #sigma_n   (noise in x)
l=0.1           #length scale
def k(x1,x2,sigmaf,l):
    "Covariance function"
    return (sigmaf**2)*exp(-((x1-x2)**2)/(2*(l**2)))

K=zeros((len(xdata),len(xdata)))
for i,x1 in enumerate(xdata):
    for j,x2 in enumerate(xdata):
        K[i][j]=k(x1,x2,sigmaf,l)
K=matrix(K)

K2=K+(sigman**2)*eye(len(xdata))
L=la.cholesky(K2,lower=True)
y2=la.solve(L,ydata)
alpha=la.solve(L.T,y2)

def fstar(xstar,xdata,alpha,sigmaf,l):
    """calculate the expected value of f for a given x"""
    kstar=zeros(len(xdata))
    for i,x in enumerate(xdata):
        kstar[i]=k(xstar,x,sigmaf,l)
    return dot(kstar.T,alpha)


ycalc=[]
for xstar in x:
    ycalc.append(fstar(xstar,xdata,alpha,sigmaf,l))


#plot(x,ycalc-y)
plot(xdata,ydata,"k-")
plot(x,ycalc,"r-")
title("sigmaf=%d, sigman=%d"%(sigmaf,sigman))
savefig("graph%d-%d.png"%(sigmaf,sigman))
Improve this answer
$\endgroup$
5
  • $\begingroup$ Please demonstrate your suggested approach on the actual data, if possible. $\endgroup$
    – bbgodfrey
    Jul 13, 2015 at 23:06
  • $\begingroup$ I could do so in python (I think), but don't have mathematica currently. $\endgroup$
    – Dr Xorile
    Jul 13, 2015 at 23:14
  • $\begingroup$ Done! One other thing to note, is that there can be "ringing" at places where there's a sharp change. In such instances, you're probably better of defining different fits before and after the inflection. $\endgroup$
    – Dr Xorile
    Jul 14, 2015 at 0:06
  • 2
    $\begingroup$ Wonders if I've ended up on the signal processing stack exchange by mistake ? $\endgroup$
    – image_doctor
    Jul 14, 2015 at 15:42
  • 1
    $\begingroup$ Anyone up for translating this into Mathematica? Everything he used in Python should also be usable here apparently. $\endgroup$
    – J. M.'s eventual burnout ♦
    Jul 15, 2015 at 16:49
1
$\begingroup$

One way is to smear the samples, that is, convolve with something fairly well localized. I show this using a Gaussian. The code is nothing special, and with some work one might be able to use Fourier transform methods in a similar way, with perhaps better results. Anyway, here goes. We start with data as given in the original post.

{otimes, ovals} = Transpose[data];
times = otimes - Min[otimes];
vals = ovals - Mean[ovals];

Here is the convolution. At the expense of computation time, one might use different Gaussians in different regions of the samples (e.g. based on time gaps between neighbors). As is, this is coded to allow a specification for the sigma argument, and also a weighting to say how large in the middle the Gaussian should be. Possibly I should just force this to be unity, to get the overall weight correct.

smoothed[t_?NumberQ, sig_, wgt_] := 
 Total[Table[
    vals[[j]]*PDF[NormalDistribution[0, sig], t - times[[j]]], {j, 
     Length[times]}]]/(wgt*PDF[NormalDistribution[0, sig], 0])

The time gaps are mostly around .035. Setting sigma to .072, roughly two such gaps, seems to give a reasonable smoothing. Possibly with more distortion than is desired, I'm not sure.

Show[Plot[smoothed[t, .072, 1], {t, 0, 12}, PlotRange -> All], 
 ListPlot[Transpose[{times, vals}], PlotRange -> All, 
  PlotStyle -> Red]]

enter image description here

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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