[SOLVED] AMATH482 HW1-Tracking a Submarine

Abstract

In this paper, Fast Fourier Transform(FFT) is used to find the center frequency and location of a moving submarine by analyzing noisy acoustic data. Each recorded time’s frequency spectrum is averaged to reduce the white noise and obtain the center frequency by applying FFT. I obtained the path of the submarine in time dimension after denoise the data with a Gaussian filter.

1 Introduction

Acoustic data detection has been widely used in security and monitoring. A new submarine technology also utilizes acoustic data to emit an unknown acoustic frequency. By applying FFT, the acoustic data is analyzed to determine a submarine’s center frequency through averaging the spectrum. Besides, the acoustic data is denoised by a Gaussian filter to locate its path. Once figuring out the trajectory of the submarine by receiving the acoustic data, P-8 Poseidon subtracking aircrafts would be sent to track submarines.

2 Theoretical Background

Fourier Series:

Any function in the interval 𝑥 ∈ [−𝐿, 𝐿] can be expressed in terms of linear combinationofharmonicallysinusoidalfunctions sin(2𝜋𝑓𝑥+𝜙),𝑐𝑜𝑠(2𝜋𝑓𝑥+𝜙) of different amplitudes(A), frequencies(f) and phases(𝜙), known as Fourier series. We write the periodic function as the infinite sum of these sines and cosines:

𝑓(𝑥)=!! +∑# (𝑎 cos(𝑘𝑥)+𝑏 sin(𝑘𝑥)),𝑥∈[−𝐿,𝐿] (1) “$%&# #

Where the coefficient 𝑎’, 𝑎#, and 𝑏# are given by the formulae
𝑎’ = & ∫( 𝑓(𝑥)𝑑𝑥 𝑎# =&∫( 𝑓(𝑥)𝑐𝑜𝑠#*𝑥𝑑𝑥 𝑏# =&∫( 𝑓(𝑥)𝑠𝑖𝑛#*𝑥𝑑𝑥

“()( ()( ( ()( (

Using exponential form of the cos(𝑥) = +”#,+$”# and sin(𝑥) = +”#)+$”# in equation “”

(1) and coefficient formulae, we can express 𝑓(𝑥) in complex exponential form. 1 (.$*)0 1 ( )(.$*)0

𝑓(𝑥)=A𝑐$𝑒 ( ,𝑥∈[−𝐿,𝐿] 𝑤h𝑒𝑟𝑒𝑐$ = 2𝐿G 𝑓(𝑥)𝑒 ( 𝑑𝑥,𝑘∈Z )1 )(

Fourier Transform:

The Fourier Transform is an integral transform defined over the interval 𝑥 ∈ [−𝐿, 𝐿]. 1

Fourier transform and its inverse are defined as

1 1 ).$0 𝐹(𝑘) = G 𝑒

𝑓(𝑥)𝑑𝑥

𝑓(𝑘) =

1 1 .$0
G 𝑒 𝐹(𝑥)𝑑𝑥

√2𝜋 )1
From the definition, we set that the transform is over the area 𝑥 ∈ [−∞, ∞] and our

special domain is in the definite main 𝑥 ∈ [−𝐿, 𝐿].
Fast Fourier Transform:
The Fast Fourier Transform (FFT) is the algorithm that transform a function of time into a function of frequency. Comparing to Discrete Fourier Transform (DFT) which has a computing complexity 𝑂(𝑁”), FFT has its advantage that it has relatively low operation complexity of 𝑂(𝑁𝑙𝑜𝑔𝑁).
The FFT algorithm performs more efficiently with large data. Also, MATLAB has a built-in function to perform the FFT for data analysis.
Averaging:
The signal data are filled with large amount of white noise which prevent us from figuring out the signals. Since the white noise could be modeled in a normal distribution and added to each Fourier component of the spectrum. We could average over many realizations in frequency space. Therefore, the noise from each realization will cancel out and we could get the central frequency.
Filtering:
There are many ways to filter the signal data in order to denoise and better detect the signal. A spectral filter around the central frequency is used in this project to remove undesired frequencies and most of the white noise. The Gaussian filter is used here:

𝐹(𝑘) = exp (−𝜏(𝑘 − 𝑘’)”)
This is a Gaussian function with 𝜏 determining the width of the filter and the

constant 𝑘’ is the target signal.
3 Algorithm Implementation and Development

In order to find the center frequency generated by the submarine and determine the path of the submarine’s movement, four components have been conducted with algorithm implemented.

1. Data Acquisition and Setup

   After the acoustic data within 49 time points is loaded, the frequency components have been setup. Since the Fast Fourier Transform assumes 2𝜋 periodic signals

   and [−𝐿, 𝐿] is the domain of input data, the frequencies are rescaled by “*. The 3- “(

   D grid [X, Y, Z] in spatial-time discretized domain and the 3-D grid [Kx, Ky, Kz]

   in frequency domain have been setup.

2. Spectral averaging

   Averaging method with 49 realizations has been applied to average the acoustic data. I applied the 3-D Fourier transform (fftn function) to transform the data from spatial domain to frequency domain and rearrange the data by shifting the zero-frequency component to the center of the array(fftshift function). I got the sum of the data from 49 times of the loop and then average the sum to offset noise. I took the absolute

√2𝜋 )1

2

value of the averaged data, since the original results in the matrix contain both real

and imaginary parts.

3. Filtering and denoising

   First, the Gaussian filter is set with replacing the k0 with the center frequency values in the 3-D. A 3*49 matrix is created to record the 3-D locations in the x, y, and z- axis in each time’s measurement. For each time point, I applied the 3-D Fourier transform(fftn function) on the acoustic data and rearrange (fftshift function) it to the center at zero frequency. Then I applied the filter to the signal data in the frequency domain. Inverse 3-D Fourier transform (ifftn function) is used to acquire the signal data in the time domain. Then I found the indexes for peaks and obtained the 3-D locations. Lastly, I used the plot3 function to plot all of the coordinates recorded to figure out the submarine trajectory.

4. Obtaining the Final Locations
   By finding the last column of the locations that record all 49 positions of the submarine, I get the final location of the submarine. Since the subtracking aircraft cannot go underwater, I give the final x and y coordinates for the submarine.

4 Computational Results

The center frequency 𝐾2 = [5.3407, −6.9115, 2.1991]
The coordinates of the submarine at each instance of time were obtained and showed in Table 1. The figure 1 shows the trajectory of the submarine.

Figure 1: The trajectory of the submarine over a 24-hour period in half-hour increments.

3

Depending on different raw data we get, either averaging or various filters could be applied to denoise the noisy data. The denoising methods with the implication of Fourier transform in the data analysis could be used in many natural science research analyses.

5

Appendix A. MATLAB Functions

1. y = linespace(x1, x2, n) returns a row vector of n evenly spaced points between x1 and x2.
2. [X, Y, Z] = meshgrid(x, y, z) returns 3-D grid coordinates defined by the vectors x, y, and z.
3. B = reshape(A, sz) reshapes A using size vector sz into B
4. Y = fftn(X) returns the multidimensional Fourier transform of an N-D array using

   a fast Fourier transform algorithm

5. Y = ifftn(X) returns the multidimensional discrete inverse Fourier transform of an

   N-D array using a fast Fourier transform algorithm

6. Y = fftshift(X) rearranges a Fourier transform X by shifting the zero-frequency

   component to the center of the array

7. [row, col] = ind2sub(sz, ind) returns the array row and col containing the

   equivalent row and column subscripts corresponding to the linear indices ind for a matrix of size sz.

6

Appendix B. Matlab Code

%% 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); % grid in frequency domain

%% problem 1 - center frequency
Utave = zeros(64, 64, 64)

for j = 1:49
  Un(:,:,:)=reshape(subdata(:,j),n,n,n);
  Utn = fftshift(fftn(Un));
  Utave = Utave + Utn;

end
Utave = abs(Utave./49);
[maxVals, indices1] = max(Utave(:));
[kx0, ky0, kz0] = ind2sub([n,n,n], indices1);
Kx0 = Kx(kx0, ky0, kz0)
Ky0 = Ky(kx0, ky0, kz0)
Kz0 = Kz(kx0, ky0, kz0)
%% problem 2 - filter
locations = zeros(3, 49);
tau = 0.2;
filter = exp(-tau*((Kx-Kx0).^2+ (Ky-Ky0).^2+ (Kz-
Kz0).^2)); % Define the filter
for j = 1:49

   Un(:,:,:)=reshape(subdata(:,j),n,n,n);
   Utn = fftshift(fftn(Un));
   Unft = filter.*Utn; % Apply the filter to the signal

in frequency space
   Unf = ifftn(Unft); % the signal in time space

   M = max(abs(Unf),[],'all');
   indices2 = find(abs(Unf)==M);

7

   [a, b, c] = ind2sub([n, n, n], indices2);
   locations(1, j) = X(a,b,c);
   locations(2, j) = Y(a,b,c);
   locations(3, j) = Z(a,b,c);

end
figure(1)
plot3(locations(1, :), locations(2, :), locations(3,:),
'ko-', 'Linewidth', 1)
title('Submarine Movement trajectory', 'Fontsize', 15)
xlabel('x')
ylabel('y')
zlabel('z')
axis([-10 10 -10 10 -10 10]), grid on
x49 = locations(1, 49);
y49 = locations(2, 49);
z49 = locations(3, 49);
L = sprintf('The 49th location of the submarine is
at %s %d %f.', x49, y49, z49);
%% Problem 3 - get the table of the 49 2-D positions in
each time point
num = [1:49]
l = [locations(1,:);locations(2,:)]
XY_locations = table(num, l)
writetable(XY_locations, 'xy.csv')
type 'xy.csv'

8