How would you draw this in MATLAB?

Hi,

I have an equation for the absorption coefficient in underwater acoustic communication channel (UWA) as:

alpha=0.11.*(f.^2./(1+f.^2))+44.*(f.^2./(4100+f))+2.75*10^-4.*f.^2+0.003;

where f is the frequency in kHz, and the absorption coefficient is given by dB/Km. I have a paper that plots this relation versus the frequency. It is straightforward as it appears at first, but when I tried to do it myself I have not gotten the same curve, even though the vertical axis in the paper is in dB/Km, which means it is just a direct substitution in the above equation. The question is why?

Thanks in advance

[G A]

The link to the paper or the figure will help us to answer your question.

[Rick Rosson]

Please post your MATLAB code.

[Jan]

All we know is the formula and that you do not get the "same curve". It might be a typo in the formula, another scaling of log-scaling, different colors or line styles, a post-script problem in the paper, etc. I do not see a chance to guess the source of the differences yet. Please post the necessary details.

Answer 1

[S. David]

This is my code:

clear all;
clc;

f=10^3.*(0:1000);

alpha=(0.11.*((f.^2)./(1+f.^2))+44.*(f.^2./(4100+f))+(2.75.*10^(-4)).*(f.^2)+0.003);

plot(f,alpha)

The figure is Fig 1 at http://www.mit.edu/~millitsa/resources/pdfs/bwdx.pdf

Thanks

Answer 2

[Jan]

Unfortunately you still did not describe, which difference you mean. I guess it is the small wobble near to 60Hz.

When I read the explanation about the Figure 1 in this paper, I find at first the above formula and then:

 This formula is generally valid for frequencies above a
few hundred Hz. For lower frequencies, the following
formula may be used:

 10 log a(f) = 0.002+ 0.11*(f.^2 ./ (1 + f.^2)) + 0.011 * f.^2

I guess, that the author used this formula for the low frequency part.

[EDITED] Check your formula again. It is "4100+f.^2" instead of "4100+f".

[S. David]

First of all, the scaling in the vertical axis is not the same as the one in the paper. If I divide the result by 10^6 the scaling will be fine, but I do not have a justification of doing so.

Despite that fact, even at high frequencies the curves are not the same. For example, in the 400 kHz the absorption coefficient in the paper is around 90 dB/Km, while it is about 60 dB/Km (after dividing by 10^6) using the above code!!

Thanks

[Jan]

See [EDITED]. Btw, 2.75*10^-4 consumes much more time than 2.75e-4.

[S. David]

Not much difference when I wrote f.^2 instead of f. Now, I need the program to work. Later I can consider more efficient one.

[Jan]

Ok, not much, but at least enough?

I suggest the standard notation of numbers mainly, because this is faster to write and to read. This is more important than the runtime.

Answer 3

[G A]

This gives the same picture as in the paper:

clear;
f=0:1000;%kHz
alpha=0.11.*f.^2./(1+f.^2)+44.*f.^2./(4100+f.^2)+2.75e-4.*f.^2+0.003;%dB/km
figure(1)
clf
plot(f,alpha,'-b')

[Jan]

Is there another difference except for the 4100+f.^2, as I've said already?

[G A]

multiplying f-array by 1e3 is not necessary

[Honglei Chen]

In your original question, it says f is in kHz so that's why you don't need to multiply 1e3.

Answer 4

[Saed]

That is great. What did you do? I mean I tried all of these stuff. Anyway, it is fine now. Thanks a lot.

[G A]

1)I have used equation given in the paper. Jan has spotted your mistake. 2)I have removed frequency multiplication by 1e3 you have done.

[S. David]

Thanks