Solving Sea Surface Height with FFT in MATLAB

[nickthequick]

Hi,

I have the following problem that I would like to solve by using the 'fft' function in matlab.

Some background on the problem: The evolution of sea surface height is given, to first order, by

$$\eta_o(x,t)=A\cos\left(kx-\omega t\right)$$

From this, we see that there are three parameters, A, k, \omega, that dictate \eta(x,t). According to linear theory, in general \eta will be a super position of these waves and will be given by

$$\eta(x,t) = \sum_n A_n \cos\left(k_nx-\omega_nt \right)$$
Each of these waves obey the deep water dispersion relationship

$$\omega_n^2=gk_n$$

Finally, to find the A_n, we note that we know \eta(0,t). The A_n are constants, therefore we can find them by inverting the relation

$$\eta(0,t)= \sum_n A_n \cos\left(\omega_nt \right)$$

which is done in MATLAB by taking a fft of \eta(0,t). The last pieces of information I need are the \omega_n s. I'm not sure what sets the frequencies of the system and this is what is giving me trouble. Of course there might be a significantly easier way to do all of this.

Any help would be appreciated,Nick

 

[Pythagorean]

not quite sure where you're stuck.

from the top, you'd define the function numerically by first defining the range. I.e., you have to define all the variables and constants in your equation. A variable is defined such that

t = 0:.1:10

would make a vector called "t" that has values form 0 to 10 in steps of .1, the units will be seconds if you choose the units of [tex]\omega[\tex] to be radians per second.

now you can define w, x, k as constants. If any of those are variables as well, this becomes a higher dimension problem.

now you can define eta:

eta = sin(k*x-w*t)

Sometimes, to get the proper vector operation on an expression, you have to put a period (.) in front of the operator, like "w.*t" instead of "w*t". So keep that in mind if you get errors.

Finally, the fft function is called:

fft_eta = fft(eta)

This is the vector of amplitudes for your fft. To construct the frequency vector, "f", you have to base it on the length of elements in eta and your sampling frequency (1/dt)

f = 0:fs/length(eta):fs/2 - fs/length(eta)

fs/2 is referred to as the nyquist frequency.

 

[nickthequick]

Pythagorean,

Thanks for the help! I guess my main issue was the seemingly opaque way that MATLAB defined the independent variable frequency for the Fourier transform. I found a nice tutorial http://blinkdagger.com/matlab/matlab-fourier-transform-with-freqz/".My code is, however, very inefficient. Does anyone see a better way to do this?

g=9.81; %gravity
Fs=50; %sampling frequency

t=(0:length(WG_all_resync.WG_mean(1,-1)/50; %time length all data was sampled
y1=packet49; %this is the output of my wave maker and is the signal I'm processing
y=[y1;zeros(7731,1)]; %pads the output to be the (temporal) length that the data is observed
[Y,f]=positiveFFT(y,Fs); %this function calculates the positive part of the fft as well as the associate frequencies
k=(2*pi*f).^2/g; %defining wavenumbers based on deep water disp relation
x=linspace(0,169,length(k)); %is this the proper way to define the x-ordinate?

eta=zeros(length(f),length(k));
for i =1:length(f)
for j=1:length(k)
eta(i,j)=real(sum(Y'.*cos(k*x(j)-2*pi*f*t(i))));
end
end

Also, I think a lot of the subtleties associated with defining the independent variables are lost on me. For instance, I'm not sure if I'm defining x and t properly. Are they dependent on the way I define f and k?

Any help would be appreciated. Thanks!Nick

 

[nickthequick]

I have a further question regarding these results.

My results are shown below.

There are two things that are bothering me. For one, MATLAB has indexed my time axis backwards, so in order to plot my time values appropriately i have to make them negative. Is there an easy way to fix this?

The second thing is that because this solution is periodic, I'm seeing artificial "noise" from the other solutions. Is there an easy way to remove these unwanted solutions so that I just display, in effect, the primary branch?

Thanks,

Nick

 

[Pythagorean]

g=9.81; %gravity
Fs=50; %sampling frequency

you can do this on one line:

g=9.81;Fs=50; %gravity and sampling rate

t=(0:length(WG_all_resync.WG_mean(1,-1)/50; %time length all data was sampled
y1=packet49; %this is the output of my wave maker and is the signal I'm processing
y=[y1;zeros(7731,1)]; %pads the output to be the (temporal) length that the data is observed
[Y,f]=positiveFFT(y,Fs); %this function calculates the positive part of the fft as well as the associate frequencies...

...

x=linspace(0,169,length(k)); %is this the proper way to define the x-ordinate?

I wouldn't hard code the 50, 7731, 169 or the 1 that designates which row of WG_mean to pick from; I'd call theme by variables and define the variables at the beginning of the code. you might have already defined the 50 as Fs (sampling frequency).

for x though, I would use a mathematical definition as I did for the frequency vector instead of relying on matlab's linspace function.

eta=zeros(length(f),length(k));
for i =1:length(f)
for j=1:length(k)
eta(i,j)=real(sum(Y'.*cos(k*x(j)-2*pi*f*t(i))));
end
end
Nick

I feel like you could do all this one line, but I have no clue how.

For instance, I'm not sure if I'm defining x and t properly. Are they dependent on the way I define f and k?

x is dependent on the length of k which is dependent on the length of f, so it's really comes down to length of y as the input to the fft function. You can pad y with zeros, for instance, to make it longer.

 

[Pythagorean]

There are two things that are bothering me. For one, MATLAB has indexed my time axis backwards, so in order to plot my time values appropriately i have to make them negative. Is there an easy way to fix this?

flip the matrix:
http://www.mathworks.com/help/techdoc/ref/flipud.html

The second thing is that because this solution is periodic, I'm seeing artificial "noise" from the other solutions. Is there an easy way to remove these unwanted solutions so that I just display, in effect, the primary branch?

not sure what you're talking about exactly, but you can reduce noise by either smoothing (matlab's smooth function is god-awefully slow though, I'd suggest the fastsmooth on the open file exchange) or by just low-pass filtering (assuming the noise is a reasonably higher frequency than your data-of-interest).

 

[nickthequick]

Thanks for the responses, Pythagorean.

My first problem was solved using set(gca,'YDir','normal');

As for the second problem, I'll try to quantify what I mean a bit more. I believe what we are observing in the previously appended plot are three solutions. Call them f(x,t), f(x,t-tau) and f(x,t+tau) where f(x,t) is the solution that I'm interested in and is the one that originates around t = 40-60 seconds. I believe that these other solutions, f(x,t-tau) and f(x,t+tau), are simply an artifact of the periodicity of solutions when using the Fourier transform. To me, these solutions represent noise and I would like to remove them so that all I have left is f(x,t).

I have been experimenting with creating filters in Fourier space to try and remove this unwanted information but so far have been unsuccessful.

As usual any help would be greatly appreciated!

Nick

 

[Pythagorean]

Well, filters will take out frequency components; but the way you've stated it, you're trying to remove temporal components? Is this correct?

you can move between temporal shifts and frequency via the "shift theorem"

 

[nickthequick]

what i ended up doing was padding my initial signal with a sufficient number of zeros such that the period becomes extremely large - this way the unwanted periodic noise is no longer an issue in the domain that i care about.

thanks for all of your help,

nick

 

[AlephZero]

The second thing is that because this solution is periodic, I'm seeing artificial "noise" from the other solutions. Is there an easy way to remove these unwanted solutions so that I just display, in effect, the primary branch?

Find out about using a "window function" on your raw data before you do the FFT. Two standard "general purpose" windows are called the Hamming and Hanning windows (invented by two different people with confusingly similar names!). There are others (Kaiser-Bessel, exponential, flat-top, etc, etc) but unless you want to get deep into this, either of the Hamming or Hanning should work pretty well and you probably won't see much difference between using one or the other.

It would be amazing if Matlab doesn't have these already implemented (but I don't use it much) - check the documentation.

 

[Pythagorean]

yeah MATLAB has a whole digital signal processing toolbox with it that includes windowing funcitons and such.