Where is my mistake in this MATLAB code?

[EngWiPy]

Hi,

I am trying to plot the AN product sec II.C in the attached paper as:

clear all
clc

s=0.5;
w=0;
d=5; 

f=0:20;

NdB1=(17-30.*log10(f));
N1=10.^(NdB1./10);
NdB2=40+20*(s-0.5)+26.*log10(f)-60.*log10(f+0.03);
N2=10.^(NdB2./10);
NdB3=50+7.5*w^0.5+20.*log10(f)-40.*log10(f+0.4);
N3=10.^(NdB3./10);
NdB4=-15+20.*log10(f);
N4=10.^(NdB4./10);

N=N1+N2+N3+N4;

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

A=(d^(1.5)).*a.^d;
A_N=1./(A.*N);
 
A_NdB=10*log10(A_N);
 
plot(f,A_NdB)

The details of these equations are highlighted in the attached file.

I appreciate if someone points to me my mistake, and the figure should be like Fig 3

Thanks

 

[NemoReally]

What problem are you seeing? I'm not that familiar with Matlab, but one potential problem is trying to take log(0) ... how does Matlab handle this (eg, error, NaN)?

I've implemented your equations in Mathcad and the shape of the resultant plot seems to agree with the document's plot, seemingly differing by an offset and/or scaling factor - the peaks seem to agree and the overall size looks OK, so I suspect an offset of about 45 dB.

 

[EngWiPy]

What is the f range in your plot? Can you post your code please?

 

[NemoReally]

What is the f range in your plot? Can you post your code please?

I'm (hopefully) plotting the same range as you, on the assumption that your code was specialized for kHz. Here's an image of the (attached) Mathcad code.

 

[EngWiPy]

Are you using log10 or natural logarithm in your code?

It looks like very much the same as my code!

 

[NemoReally]

Yes, it does look pretty much like your code, so that's why I was wondering what problem you're seeing.

Mathcad's "log" is log10 by default; Mathcad uses "ln" for loge. Whereas, it seems, Matlab's "log" is loge. Using loge results in plots that look markedly different from the article's plots - far more peaky.

 

[uart]

Are you using log10 or natural logarithm in your code?

It looks like very much the same as my code!

There's no real problem with your code, it's just the normalizing constant $A_0$ causing the dB offset. That and the fact that you are plotting at too few values of "f", especially at low frequencies, to get a very accurate picture of what's going on.

Just replace the line

"f = 0:20"

with

"f=[0.1 : 0.1 : 20]" (or even better with f = logspace(-1,log10(20),100)

and you'll get exactly the same graph as figure 3 in the paper (k=1.5, w=0, s = 0.5 and d=5km), except for about a 45 dB offset, presumably due to A0.

BTW. You could have made it a lot easier for people to see what your problem was if you had of actually described the errors or anomalies you were concerned about (or even better posted an image of the plot you were obtaining).

 

[EngWiPy]

There's no real problem with your code, it's just the normalizing constant $A_0$ causing the dB offset. That and the fact that you are plotting at too few values of "f", especially at low frequencies, to get a very accurate picture of what's going on.

Just replace the line

"f = 0:20"

with

"f=[0.1 : 0.1 : 20]" (or even better with f = logspace(-1,log10(20),100)

and you'll get exactly the same graph as figure 3 in the paper (k=1.5, w=0, s = 0.5 and d=5km), except for about a 45 dB offset, presumably due to A0.

BTW. You could have made it a lot easier for people to see what your problem was if you had of actually described the errors or anomalies you were concerned about (or even better posted an image of the plot you were obtaining).

I am sorry if I were not clear, and thanks for your reply.

For d=5 Km the curves are very similar except for a certain scale, but for d=100 Km there is a huge difference! Is that because of A0, too?

Thanks

 

[uart]

I am sorry if I were not clear, and thanks for your reply.

For d=5 Km the curves are very similar except for a certain scale, but for d=100 Km there is a huge difference! Is that because of A0, too?

Thanks

Ok I don't know, I only checked the code you posted and it gave the correct results (except for the offset) for d=5. The paper refers to a normalizing constant A0, which is not included in your equations, so certainly that could account for the dB offset. The other problem was that the result is -inf dB at f=0, so that result wasn't being shown on the graph, and plotting at just the integer frequencies 1:20 was not enough to get a good picture of the graph.

 

[uart]

Ok I just checked the result at d=100 and it also agrees exactly with the graph (fig3) in the paper (except for the offset due to whatever normalizing constant they're using).

For d=100 I used the freq range: f = logspace(-2, log10(7), 100) to reproduce the figure in the paper (same shape exactly, but again with an offset of about 45 dB).

 

[EngWiPy]

Ok I just checked the result at d=100 and it also agrees exactly with the graph (fig3) in the paper (except for the offset due to whatever normalizing constant they're using).

For d=100 I used the freq range: f = logspace(-2, log10(7), 100) to reproduce the figure in the paper (same shape exactly, but again with an offset of about 45 dB).

Thank you so much, it is working now. Just the scaling factor is different which is of minor importance.