view general/numerical/frint.m @ 16:db7f4afd27c5

Rearranging numerical toolbox.

function y=frint(r,x,dc)
% frint - Fractional integration
%
% y=frint(r,x[,dc])
% Returns the rth order integral of x, where r can be fractional
% if dc=0 (default is 1) dc term is not added back in

x=x(:);
n=length(x);
z=fft(x); z0=z(1); z(1)=0;
K=(0:n-1)';
D=exp(2*pi*i*K/n)-1;
D(1)=1; % to avoid divide by zero
z2=z./D.^r;
y=real(ifft(z2)); % integral minus DC term

% put DC term  back as a ramp
if nargin<3 || dc==1,
	y=y+sign(z0)*(K*abs(z0)/n).^r;
end