FFT

Good Afternoon,

A question on the implementation of the FFT when applied to a distortion. I am using a set of polar coordinates for a cylinder that has been distorted. Background: the variation in geometry around the circumference of a distorted cylinder bore may be modeled by a general fourier series. R(theta)=sum(a_i cos(theta) +b_i sin(theta) from 0 to n where n is the order of distortion. This enables one to use the FFT algorithm. When implementing this idea: I first convert the nodal displacements into polar coordinates and then oversample said coordinates to use in an interpolation function to generate a function used for the fft. I’ll call this interpolation function p. So

%Calculate fft
Fourierfun=fft(p)
%Calculating DC component and Ntheta is number of points used
a0=2.*p(1)/Ntheta

%Calculating coefficients
for ik = 1:Ncoeffs
    ij = ik +1;       %index for fft

%Cosine Terms (ak)
ak(ik)=2.*(-1)^(ik).*real(p(ik))/Ntheta
%Sine Terms
bk(ik)=2.*(-1)^(j).*p(ij))/Ntheta

I think the order of distortion is sqrt(ak.^2+bk.^2) which will get me the order.

Does this seem like the correct methodology for amplitude of order of distortion?

Answer 1

[Dr. Seis]

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;

[Dr. Seis]

I used information from the following links:

http://en.wikipedia.org/wiki/Fourier_coefficient

http://en.wikipedia.org/wiki/Ellipse

However, there is a difference between how the Fourier coefficients are supposed to be used (see 3rd equation under definition) and how I implemented it to get it to work. I am unsure why I did not need to halve a0 when defining "Radius_new", but that is what worked.

[Dr. Seis]

Actually, it looks like using the FFT a0 is associated with the amplitude at 0 frequency, which is the mean of the signal [i.e., a0 should equal mean(Radius) or mean(profile)]. However, the way a0 is calculated using the equation using the wiki link, you need to do:

a0 = sum(Radius)*Theta_Inc/pi

or

a0 = sum(profile)*Theta_Inc/pi

If a0 is computed this way, then you do need to divide a0 by two in order to get the mean of the Radius.

[Melissa]

Oh wow so much information. Thank you so much for going out of your way to assist me. I am going to go through everything and try to comprehend it. Thank you!!!

[Melissa]

Once concern for trying to calculate the order, the Ncoeffs is the number of terms in the series so when I change the number of terms, it does not change the order. For example if I plot with Ncoeffs of 3 or 4 the picture is the same.

[Dr. Seis]

When I change Ncoeffs from 4 to 3, I get a much different image - it goes from a really close match to basically a circle. These are my coefficients (rounded to four decimal places):

a0 = 3.3238

ak = [0.000,0.000,0.000,0.6756];

bk = [0.000,0.000,0.000,0.0000];

So if you try to reconstruct the bore cylinder with anything less that Ncoeff = 4, you will basically just get a cylinder with radius approx. equal to a0.

Answer 2

[Dr. Seis]

Ok... let's start over. I found a link to the "fcoeffs_fft" code (located here: http://pastebin.com/001id7SB). Are you trying to modify this code to suite your needs? I noticed that the equations for "ak" are different between the website I listed and the one you listed before, which showed all but the last two lines of code (i.e., <http://www.mediafire.com/imageview.php?quickkey=439h3j30080hh2k&thumb=6>).

[Melissa]

Thats the one I'm trying to modify, the sampling points of 1-24 did not seem like enough to accurately generate the coeefficients but now I think that would be sufficient. I'm using a set polar coordinates as my profile and angles and then oversampling those points to use in the fft. My concerns are with the equations of calculating the fourier series coefficients, I don't understand how those were derived. I am attaching a link with the full code with what I think is the correct terms to use for the FFT. If you could tell me if I'm correct with this or tell me where I'm wrong so I can correctly learn the FFT it would be greatly appreciated.

http://www.mediafire.com/i/?9vv7djb2yur3m1x

[Dr. Seis]

I don't think you need the (-1)^(ij). I tried the original code and it seemed to work for me (it reconstructed an ellipse). I haven't tried anything more complicated through. I will write out an explanation of what is going on a little later tonight.