So, I created an example using the original un-modified code for fcoeffs_fft (located here: http://pastebin.com/001id7SB) and was able to get it to work. It looks like you have to define your cylinder bore on the interval from -pi to +pi (if you define it from 0 to 2*pi you get incorrect results). I only made a slight modification to the code (I removed some negative signs that effectively cancel from the bk terms):
function [a0 ak bk]=fcoeffs_fft(profile,angles,Ncoeffs)
%FCOEFFS_FFT Fourier coefficients (trigonometric) via Fast Fourier Transform (FFT)
% [a0,ak,bk]=fcoeffs_fft(profile,angles,Ncoeffs) determines the Fourier coefficients of
% irregulary spaced angles and profile.
% Ncoeffs is the number of Fourier coefficients to return
% Oversamples input via a (periodic) cubic spline interpolation
% Returns DC component a0, cosine terms ak and sine terms bk
%Oversampling the data points so the total number of points is Ntheta
Ntheta = 1024;
theta = linspace(-pi,pi,Ntheta);
delta_theta = 2*pi/Ntheta;
%This defines a temporary vector that "wraps around" pi to -pi
theta_over = -pi:delta_theta:(pi + 5*delta_theta);
%Creating the function to which the fft can be applied
pp = spline(angles,profile);
zeta_temp= ppval(pp,theta_over);
%Overwrite here
zeta = zeros(size(theta));
zeta(1:Ntheta) = zeta_temp(1:Ntheta);
zeta(1:6)=zeta_temp(Ntheta:Ntheta +5);
zetahat=fft(zeta);
a0 = zetahat(1)/Ntheta; %DC component
ak = zeros(1,Ncoeffs);
bk = ak;
for ik = 1:Ncoeffs
ij = ik +1; %index for fft
ak(ik)=2*real(zetahat(ij))/(Ntheta); %cosine terms
bk(ik)=2*imag(zetahat(ij))/(Ntheta); %sine terms
end
Once you save the above in a *.m file. You can run the following example from the command line:
% Define Sampling information
Theta_Inc = 2*pi/360;
N_Theta = 360;
Theta_Range = linspace(-180,180,N_Theta+1)*Theta_Inc;
Theta_Range = Theta_Range(1:N_Theta);
% Define Radius information for borehole shape and plot
rA = 4; rB = 2;
Radius = (rA*rB)./sqrt((rB*cos(Theta_Range)).^2+(rA*sin(Theta_Range)).^2);
Radius = max(Radius,...
(rA*rB)./sqrt((rA*cos(Theta_Range)).^2+(rB*sin(Theta_Range)).^2));
figure;plot(Radius.*cos(Theta_Range),Radius.*sin(Theta_Range)); axis equal;
% Determine Fourier coefficients
Ncoeffs = 4;
profile = Radius;
angles = Theta_Range;
[a0 ak bk]=fcoeffs_fft(profile,angles,Ncoeffs);
% Determine new radius information from Fourier coefficients
% of order Ncoeffs
Radius_new = ones(size(Radius))*a0;
for ii = 1 : Ncoeffs
Radius_new = Radius_new + ak(ii)*cos(angles*ii) + bk(ii)*sin(angles*ii);
end
% Plot comparison
figure;plot(Radius.*cos(Theta_Range),Radius.*sin(Theta_Range)); hold on;
plot(Radius_new.*cos(Theta_Range),Radius_new.*sin(Theta_Range),'r--');
hold off; axis equal;
