Fast Fourier Transform – Intelligent Additive Manufacturing

The Fourier transforms are used to analyze noisy acoustic data. It could describe any periodic function by sums of cosines and sines. It has a good application in data analysis and noise elimination. For example, we can use it to eliminate the noise of an image.

On this page, we will introduce its theoretical background and give an example of its application which is to detect a submarine from an area. For more details and questions, please feel free to contact us.

Theoretical Background

Fourier Transform

The Fourier Transform has equation:

(1) $\begin{equation*} f(x)=\frac{a_0}{2} + \sum_{i=1}^{\infty} \left( a_n \cos nx + b_n \sin nx \right) \qquad \left( -\pi, \pi \right) \end{equation*}$

In complex representation, we have:

(2) $\begin{equation*} f(x)= \sum_{-\infty}^{+\infty} C_n e^{ \frac{ik \pi x}{L}} \qquad x \in \left( -L, L \right)\end{equation*}$

(3) $\begin{equation*} C_n =\frac{1}{2L} \int_{-L}^{L} f(x) e^{ \frac{ik \pi x}{L}} dx\end{equation*}$

In this page, we need an expansion of equation (1) to make the interval from $x\in \left( -\pi ,\pi \right)$ to $x\in \left( -\infty ,\infty \right)$ . So, we use the following:

(4) $\begin{equation*} F \left( k\right) =\frac{1}{\sqrt{2\pi } } \int^{\infty }_{-\infty } {f(x)e^{-ikx}} dx\end{equation*}$

And

(5) $\begin{equation*} f\left( x\right) =\frac{1}{\sqrt{2\pi } } \int^{\infty }_{-\infty } {F(x)e^{ikx}} dx\end{equation*}$

2. Fast Fourier Transform

Also, in the example of this page, Fast Fourier Transform is used to handle noisy data. It has many advantages. It not only has good accuracy but also has a low operation count. In FFT, the points are $2^n$ , so the operation only is $O(N \log N)$ .

3. Gaussian Filter

The Gaussian filter plays a key role in the work of the location of a submarine. Through averaging of the spectrum, we can determine the frequency signature. By filter the data around the center frequency, it will have a good denoise and shows the path of the submarine.

The function of the filter is:

(6) $\begin{equation*} G= e^{-nk^2}\end{equation*}$

Result

Figure 1 shows the Initial noisy data.

Figure 1: Initial noisy data

Figure 2 shows the average of Fourier transformed measurements. The large green mass is the maximum value.

Figure 2: averaging of the spectrum

Figure 3: Submarine Path

Figure 3 shows the path of detected submarine.

Summary and Conclusion

5 Summary and conclusion

By using Fast Fourier Transform and Gaussian Filter, we successfully locate the frequency signature. After the filter of inverse transform of the data, the path of the submarine is found. Besides, we could use this skill to analyze SEM images and identify the growth of defects.

Code & Raw Data

% Clean workspace 
clear all; close all; clc 
load subdata.mat % Imports the data as the 262144x49 (space by time) matrix called subdata

L = 10; % spatial domain 
n = 64; % Fourier modes
x2 = linspace(-L,L,n+1); x = x2(1:n); y =x; z = x;
k = (2*pi/(2*L))*[0:(n/2 - 1) -n/2:-1]; ks = fftshift(k);
[X,Y,Z]=meshgrid(x,y,z); 
[Kx,Ky,Kz]=meshgrid(ks,ks,ks); 

%% Reshape data by resolution
Un(:,:,:)=reshape(subdata(:,1),n,n,n); 
M = max(abs(Un),[],'all');
close all, isosurface(X,Y,Z,abs(Un)/M,0.5) 
axis([-15 15 -15 15 -15 15]), grid on, drawnow 
xlabel('X'), ylabel('Y'), zlabel('Z')

%% Averaging of the spectrum
Uave=zeros(n,n,n);

for jj=1:49
    U=reshape(subdata(:,jj),n,n,n);
    Ut=fftn(U);
    Uave=Uave+Ut;
end
Uave=abs(fftshift(Uave))./49;

figure(2)
isosurface(X,Y,Z,abs(Uave)./max(max(max(abs(Uave)))), 0.3)
axis([-15 15 -15 15 -15 15]), grid on, drawnow 
title('Wavenumber')
xlabel('ut'),ylabel('ut'),zlabel('ut') 
%% Guassian filter
[M,I]=max(Uave(:));
[R,C,V]=ind2sub(size(Uave),I);%find center frequency
kx=Kx(R,C,V); ky=Ky(R,C,V); kz=Kz(R,C,V);
filter=exp(-0.2*((Kx-kx).^2 + (Ky-ky).^2 + (Kz-kz).^2)); 
% Apply filter
xp=zeros(1,jj); yp=zeros(1,jj); zp=zeros(1,jj);
for jj=1:49
     Un2(:,:,:)=fftn(reshape(subdata(:,jj),n,n,n));
     Unft=filter.*fftshift(Un2);
     Unus=ifftshift(Unft);
     Unf=ifftn(Unus);
     [M,I]=max(abs(Unf(:)));
     [R2,C2,V2]=ind2sub(size(Unf),I);
    xp(jj)=X(R2,C2,V2);yp(jj)=Y(R2,C2,V2);zp(jj)=Z(R2,C2,V2);
end
figure(3)
plot3(xp,yp,zp,'-o','MarkerFaceColor','r','Linewidth',[1]), grid on;
title('Trajectory of Submarine')
xlabel('X')
ylabel('Y')
zlabel('Z')

data download

Welcome to Intelligent Additive Manufacturing Lab

Theoretical Background

Result

Summary and Conclusion

Code & Raw Data