Borrar filtros
Borrar filtros

Can I smoothen the cdfplot() result?

8 visualizaciones (últimos 30 días)
okoth ochola
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
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
as mentionned by @Alexander, sharing the data would be a plus...
Mathieu NOE
Mathieu NOE el 25 de Mzo. de 2024
hello
so have you solved your problem ?

Iniciar sesión para comentar.

Respuesta aceptada

Alexander
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
okoth ochola
okoth ochola el 26 de Mzo. de 2024
Thank you all fro these wonderful answers. You dont how much have been helped
Alexander
Alexander el 26 de Mzo. de 2024
Cheers and upvote Mathieu also.

Iniciar sesión para comentar.

Más respuestas (1)

Mathieu NOE
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

Categorías

Más información sobre Get Started with Curve Fitting Toolbox en Help Center y File Exchange.

Productos


Versión

R2018a

Community Treasure Hunt

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

Start Hunting!

Translated by