43 views (last 30 days)

Show older comments

Image Analyst
on 5 Oct 2015

Star Strider
on 4 Oct 2015

I’m not exactly certain what you want, but one approach would be to use the cumtrapz function:

x = [0.02 0.02 0.03 0.03 0.03 0.04 0.04 0.05 0.05 0.05];

y = [0.23 0.40 0.12 0.16 0.09 0.50 0.02 0.33 0.10 0.08];

I = cumtrapz(x,y); % Integrated Vector Of ‘y’ w.r.t. ‘x’

figure(1)

plot(x, I)

grid

Star Strider
on 5 Oct 2015

I would experiment with the filter order and use filtfilt if you have the Signal Processing Toolbox. This gave what appeared to me to be close to what you want:

N = 50; % Length Of Filter

b = ones(1,N); % Filter Coefficients

a = N; % Filter Coefficients

sf = filtfilt(b, a, Vsorted); % Filter Data

I experimented with filter lengths ranging from 50 to 200 and it seemed to give the result you illustrated. It is not possible to design a typical low-pass filter because your data are not regularly sampled. It would be necessary to interpolate your data (using interp1 if your data meet its requriements) to a regularly-sampled time base (here Rsorted) to design a filter for it. That could make the curve smoother, but I can’t promise that it would be.

The first step after interpolation would be to do a fft in order to determine what frequencies you want to exclude as noise. When you decide on a cutoff frequency for your low-pass filter, I would then use the filter design procedure I outlined here.

arich82
on 5 Oct 2015

Edited: arich82
on 5 Oct 2015

[edited to fix typo in comments, and clean-up useless clutter]

Perhaps try combining the filter solution with the accumarray (or filtfilt, if you have the appropriate toolbox, which I do not):

x = Rsorted(:).';

y = Vsorted(:).';

% compute average y for each unique x

[c, ic, ia] = unique(x);

sums = accumarray(ia, y);

wgts = accumarray(ia, ones(size(y))); % weights, i.e. number of reps

avgs = sums./wgts;

% interpolate data on regular grid

n = numel(c); % or some other number

xq = linspace(c(1), c(end), n);

yq = interp1(c, avgs, xq);

% filter gridded data

windowSize = 50;

a = windowSize;

b = ones(1, windowSize);

yf = filter(b, a, yq);

% plot results

hold all;

plot(xq, yf, 'r-', 'LineWidth', 3)

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!