Can I smoothen the cdfplot() result?
8 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
okoth ochola
el 11 de Mzo. de 2024
Comentada: Alexander
el 26 de Mzo. de 2024
Hi, suppose I have a sample data x, and i want to plot the the cdf function on the data and plot it, ofcoz I will use cdfplot(x) and I will a graph similar to the one attached below. My question is, can i smoothen it, is there a way I can make the plot smoother and not look like a a step funtion plot?
3 comentarios
Mathieu NOE
el 11 de Mzo. de 2024
even with smoothing , I doubt you can make those large flats disappear
wuld probably rather try to fit either a polynomial or exponential function
Respuesta aceptada
Alexander
el 25 de Mzo. de 2024
@Mathieu NOE, thank you for remaindig me. This question kept nagging me, but I just forgotten it. Here are some rudimentary aproches. As the data were artificially generated, with real measurement the solution are only a rough direction what can be done and needed a lot of refinement.
commandwindow;
y = [];
for ii = 1:20
y = [y ones(1,100*ii)*ii];
end
n = length(y);
x = linspace(1,25,n);
plot(x,y);grid minor;
% 1. Polynomial Fit, 2. Order
% As suggested by @Mathieu NOE
P = polyfit(x,y,2);
yNew = polyval(P,x);
plot(x,y,x,yNew);grid minor; title('Polyfit');pause;
% 2. Logarithmic Fit
% Could be done much easier with lsqcurvefit, if you have optimization
% toobox. I don't have it, hence, it's handcrafted.
yLogSum = sum(y.*log(x)); ySum = sum(y); sumLnx = sum(log(x));
Zaeler = yLogSum - (ySum*sumLnx)/n;
Nenner = sum(log(x).^2) - sumLnx^2/n;
b = Zaeler/Nenner;
a = ySum/n - b*sumLnx/n;
yLogRes = a + b*log(x);
plot(x,y,x,yLogRes);grid minor; title('Logarithmic Fit');pause;
% 3. Filtering
% I just used a butterworth filter. As more longer the steps are, as higher
% is the weaviness of the output. The filter parameters has to be adjusted
% to fit the input curve. Another filter type might fit better.
[B,A] = butter(2,0.0005);
yFF = filtfilt(B,A,y);
plot(x,y,x,yFF);grid minor; title('Butterworth filtered');pause;
% 4. Step/Edge Evaluation
% The input was generated. With measured data the following lines must be
% adjusted.
% Also the data has to be interpolated between the x-values. yDiff has now
% as many data values as there are steps in the measurement. Maybe with
% interp().
yDiff = diff(y);
yFind = find(yDiff > 0.5);
Edge = 1; %Upper Edge = 1; lower edge = 0;
plot(x,y,x(yFind+Edge),y(yFind+Edge));grid minor; title('Edge Evaluation');
legend('Original',['Edge ' num2str(Edge)], 'location','northwest')
4 comentarios
Más respuestas (1)
Mathieu NOE
el 26 de Mzo. de 2024
some more code to look at - I kept only the best solutions , even though there are probably lot of other solutions
also I introduce some randomness in the data to see how robust is the code
enjoy !
y = [];
for ii = 1:20
y = [y ones(1,round(10*ii+25*rand))*ii];
end
n = length(y);
x = linspace(1,25,n);
% add some noise on y
y = y + 0.05*randn(size(y));
% 1. Polynomial Fit,
% As suggested by @Mathieu NOE
P = polyfit(x,y,6);
yNew = polyval(P,x);
figure,plot(x,y,x,yNew);grid minor; title('Polyfit');
% 2. Step/Edge Evaluation
% The input was generated. With measured data the following lines must be
% adjusted.
% Also the data has to be interpolated between the x-values. yDiff has now
% as many data values as there are steps in the measurement. Maybe with
% interp().
yDiff = diff(y);
yFind = find(yDiff > max(yDiff)/2);
Edge = 0; %Upper Edge = 1; lower edge = 0;
xl = x(yFind+Edge);
yl = y(yFind+Edge);
Edge = 1; %Upper Edge = 1; lower edge = 0;
xu = x(yFind+Edge);
yu = y(yFind+Edge);
% Median curve
xm = (xu+xl)/2;
ym = (yu+yl)/2;
figure,plot(x,y,xl,yl,xu,yu,xm,ym);grid minor; title('Edge Evaluation');
legend('Original','Lower Edge','Upper Edge','Median curve' , 'location','northwest')
% 3. Envelop (secant method)
N = 200;
[yu] = env_secant(x, y, N, 'top'); % (see function attached)
[yl] = env_secant(x, y, N, 'bottom'); % (see function attached)
% Median curve
ym = (yu+yl)/2;
figure,plot(x,y,x,yl,x,yu,x,ym);grid minor; title('Envelop (secant method)');
legend('Original','Lower Envelop','Upper Envelop','Median curve' , 'location','northwest')
%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [env] = env_secant(x_data, y_data, view, side)
% Function call: env_secant(x_data, y_data, view, side)
% Calculates the top envelope of data <y_data> over <x_data>.
% Method used: 'secant-method'
% env_secant() observates the max. slope of about <view> points,
% and joints them to the resulting envelope.
% An interpolation over original x-values is done finally.
% <side> ('top' or 'bottom') defines which side to evolve.
% Author: Andreas Martin, Volkswagen AG, Germany
if nargin == 0
test( false );
test( true );
return
end
if nargin < 4
error( '%s needs at least 4 input arguments!', mfilename );
end
assert( isnumeric(view) && isscalar(view) && view > 1, ...
'Parameter <view> must be a value greater than 1!' );
assert( isvector(x_data) && isnumeric(x_data) && all( isfinite( x_data ) ), ...
'Parameter <x_data> has to be of vector type, holding finite numeric values!' );
assert( isvector(y_data) && isnumeric(y_data) && all( isfinite( y_data ) ), ...
'Parameter <y_data> has to be of vector type, holding finite numeric values!' );
assert( isequal( size(x_data), size(y_data) ), ...
'Parameters <x_data> and <y_data> must have same size and dimension!' );
assert( ischar(side), ...
'Parameter <side> must be ''top'' or ''bottom''!' );
switch lower( side )
case 'top'
side = 1;
case 'bottom'
side = -1;
otherwise
error( 'Parameter <side> must be ''top'' or ''bottom''!' );
end
sz = size( x_data );
x_data = x_data(:);
x_diff = diff( x_data );
x_diff = [min(x_diff), max(x_diff)];
assert( x_diff(1) > 0, '<x_data> must be monotonic increasing!' );
y_data = y_data(:);
data_len = length( y_data );
x_new = [];
y_new = [];
if diff( x_diff ) < eps( max(x_data) ) + eps( min(x_data) )
% x_data is equidistant
search_fcn = @( y_data, ii, i ) ...
max( ( y_data(ii) - y_data(i) ) ./ (ii-i)' * side );
else
% x_data is not equidistant
search_fcn = @( y_data, ii, i ) ...
max( ( y_data(ii) - y_data(i) ) ./ ( x_data(ii) - x_data(i) ) * side );
end
i = 1;
while i < data_len;
ii = i+1:min( i + view, data_len );
[ m, idx ] = search_fcn( y_data, ii, i );
% New max. slope: store new "observation point"
i = i + idx;
x_new(end+1) = x_data(i);
y_new(end+1) = y_data(i);
end;
env = interp1( x_new, y_new, x_data, 'linear', 'extrap' );
env = reshape( env, sz );
function test( flagMonotonic )
npts = 100000;
y_data = cumsum( randn( npts, 1 ) ) .* cos( (1:npts)/50 )' + 100 * cos( (1:npts)/6000 )';
if flagMonotonic
x_data = (1:npts)' + 10;
else
x_diff = rand( size( y_data ) );
x_data = cumsum( x_diff );
end
view = ceil( npts * 0.01 ); % 1 Percent of total length
env_up = env_secant( x_data, y_data, view, 'top' );
env_lo = env_secant( x_data, y_data, view, 'bottom' );
figure
plot( x_data, y_data, '-', 'Color', 0.8 * ones(3,1) );
hold all
h(1) = plot( x_data, env_up, 'b-', 'DisplayName', 'top' );
h(2) = plot( x_data, env_lo, 'g-', 'DisplayName', 'bottom' );
grid
legend( 'show', h )
end
end
0 comentarios
Ver también
Categorías
Más información sobre Get Started with Curve Fitting Toolbox en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!