spect = 2*abs(N)/length(N);
N is (likely) complex-valued. spect takes abs(N) so it contains only real non-negative values.
and the square of real non-negative values is real and non-negative.
PS_coeff is not used later in the code, but PS is used. If we assume that the actual code assigns to PS instead of to PS_coeff then we can continue analyzing.
Fs is read from the file. Presumably it is Sampling Frequency, and presumably it is non-negative. But could df be less than 1? Yes; for example sampling 2 seconds at 8000 Hz would be 16000 samples and 8000/16000 would be df = 1/2 . This will probably happen most of the time. Could df be more than 1? It can in theory: for example 1/2 second of 512k Hz would be 512000/256000 = 2. So we have to prepare for the possibility of df being less than 1 (often) or more than 1 (could happen in general, but you might have reason to know that it does not happen in any case you are interested in.) Worst case: length(n) could be 0 so df could be infinite.
Pnoise=sum(PS(1:round(0.67/df)))/round(0.67/df) + sum(PS(1+round(150/df):end))/(round(500/df)-(1+round(150/df))) + PS(round(50/df));
Each of the PS() is selecting a subset out of PS, each entry of which is non-negative. sum() of a set of non-negative values must be non-negative (well, except for the fact that it might be empty)
So, term by term:
sum(PS(1:round(0.67/df))) must be real and non-negative (but could be empty if df > 4/3 meaning that the signal was less than 3/4 of one period)
df is non-negative (because we assume that Fs is positive). round(0.67/df) will be 0 or positive -- in particular not negative. A positive value divided by a non-negative value is non-negative.
Have you considered, by the way, the possibility of using mean(PS(1:round(0.67/df))) ?
sum(PS(1+round(150/df):end))/(round(500/df)-(1+round(150/df)))
again sum() of non-negative values is non-negative. Let us suppose for a moment that the denominator is non-negative (we will return to that point later), then the result would be a positive scalar.
That would be non-negative again, and a scalar.
So you would be adding non-negative values, and getting either empty (depending on df) or non-negative.
Psignal=sum(PS(1+round(0.67/df):1+round(149/df)));
sum would be non-negative (or empty) again.
SNR=10*log10(Psignal/Pnoise)
so that would be log10() of two non-negative numbers. Without further analysis, we cannot know if the Psignal is greater than or less than Pnoise, so at the moment we cannot really guess whether the log10() is negative or positive.
But now let us go back to
sum(PS(1+round(150/df):end))/(round(500/df)-(1+round(150/df)))
I said earlier "Let us suppose for a moment that the denominator is non-negative" -- but could it be negative ?
Is it possible for round(500/df) < (1+round(150/df)) ? That would require round(500/df) - round(150/df) < 1 . With 500 being greater than 150, you would tend to think that 500/df must surely be greater than 150/df . But... remember the round() . Indeed, the expression drops to 0 between df = 200 (exclusive) and df = 300 (inclusive) then up to +1 again between 300 and 1000/3, then down to 0, and drops to -1 for df > 1000 . So you could be adding a negative term here.
Why isn't this more of a problem? Well, by then you would be dealing with empty arrays -- and the PS(round(50/df)) is going to give an indexing error when df exceeds 100.
It is not at all clear how you arrived at round(500/df)-(1+round(150/df)) but my suspicion is that, once more, you should be converting the sum and division into a mean() call.
