Integral average of y values with respect to x values

조회 수: 7 (최근 30일)
matnewbie
matnewbie 2015년 10월 4일
댓글: Star Strider 2015년 10월 7일
I have two vectors (x and y) of the same size with x having more than one value in the y direction. For example, x is [0.02 0.02 0.03 0.03 0.03 0.04 0.04 0.05 0.05 0.05] and y is [0.23 0.40 0.12 0.16 0.09 0.50 0.02 0.33 0.10 0.08]. I want to convert the attached scattered plot of the two vectors to a simple continuous curve, by evaluating the integral average of y values with respect to x values. Take into account that the x vector is already ordered from small to large values. How could I do this? Thank you in advance.
  댓글 수: 1
Image Analyst
Image Analyst 2015년 10월 5일
"I have two vectors (x and y) of the same size with x having more than one value in the y direction" Huh? Does it have one more, or are x and y the same length? Not sure what you're saying. Are you saying that for any given x value, it may be in the x array more than once, with each occurrence of that x value having a different y value stored at corresponding indexes in the y array?

댓글을 달려면 로그인하십시오.

채택된 답변

Star Strider
Star Strider 2015년 10월 4일
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
  댓글 수: 7
matnewbie
matnewbie 2015년 10월 5일
편집: matnewbie 2015년 10월 5일
You can find the "input.csv" file with the data attached, and here is the code I use to load and sort the data.
clc,close all,clear all
fileID=fopen('input.csv');
C=textscan(fileID,'%f,%f,%f,%f,%f,%f,%f,%f');
A=cell2mat(C); % loads the values in A
fclose(fileID);
N=size(A,1);
grid_size=3e-4; %size of grid element
Rcol=0.0105; %radius
X=A(:,7)-Rcol-grid_size; %translation of x coordinate
Y=A(:,8)-Rcol-grid_size; %translation of y coordinate
r=sqrt(X.^2+Y.^2); %radial position
axial_vel=A(:,2); %velocity component in the axial direction
% Cutting off values outside the region of interest
tol=1e-8; %tolerance for velocity cutoff
j=1;
for i=1:N
if axial_vel(i)>tol
R(j)=r(i);
V(j)=axial_vel(i);
j=j+1;
end
end
figure
plot(R,V,'r.')
[Rsorted, SortIndex] = sort(R);
Vsorted = V(SortIndex);
figure
plot(Rsorted,Vsorted,'b.')
The idea is to obtain a smoothed continuous curve for the integral average... I tried the filter function, but I can't get what I need.
Star Strider
Star Strider 2015년 10월 5일
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.

댓글을 달려면 로그인하십시오.

추가 답변 (1개)

arich82
arich82 2015년 10월 5일
편집: arich82 2015년 10월 5일
[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)
  댓글 수: 4
matnewbie
matnewbie 2015년 10월 7일
Yes, you're right: with filtfilt the results are more similar to those obtained by Star Strider (magenta line in the plot).
Star Strider
Star Strider 2015년 10월 7일
In my latest Comment (5 Oct 2015 at 21:09), I used filtfilt.

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 Digital Filtering에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by