samer@4: function y=frint(r,x,dc)
samer@4: % frint - Fractional integration
samer@4: %
samer@4: % y=frint(r,x[,dc])
samer@4: % Returns the rth order integral of x, where r can be fractional
samer@4: % if dc=0 (default is 1) dc term is not added back in
samer@4: 
samer@4: x=x(:);
samer@4: n=length(x);
samer@4: z=fft(x); z0=z(1); z(1)=0;
samer@4: K=(0:n-1)';
samer@4: D=exp(2*pi*i*K/n)-1;
samer@4: D(1)=1; % to avoid divide by zero
samer@4: z2=z./D.^r;
samer@4: y=real(ifft(z2)); % integral minus DC term
samer@4: 
samer@4: % put DC term  back as a ramp
samer@4: if nargin<3 || dc==1,
samer@4: 	y=y+sign(z0)*(K*abs(z0)/n).^r;
samer@4: end
samer@4: 
samer@4: