octave-dev

[OctDev] pwelch bugs

[Dmitri A. S. <das...@gm...>]

The new version of pwelch (by  Peter V. Lanspeary) seems to have problem
parsing argument list. E.g. "demo pwelch" would fail:

pwelch example 1: failed
error: pwelch: arg 3 (overlap) must be real from 0 to 0.950000

apparently on the line pwelch(x,[],[],[],[],0.9);

It also gets easily confused when pwelch called with two channels, e.g.:

pwelch(x, y, [], [], [], 1000, [], [], "mean", [], "trans")
error: pwelch: can give confidence interval for x power spectrum only
error: evaluating if command near line 507, column 3
error: evaluating if command near line 279, column 1
error: called from `pwelch' in file
`/usr/share/octave/site/m/octave-forge/signal/pwelch.m'

In fact it appears that i cannot specify any parameter between
"sampling frequency" and "result" parameters...

Sorry, I do not have a patch, the source of pwelch.m is too
unwieldy for me.

Sincerely,

Dmitri.
--

Re: [OctDev] pwelch bugs

[dbateman <dba...@fr...>]

Dmitri A. Sergatskov wrote:
> 
> The new version of pwelch (by  Peter V. Lanspeary) seems to have problem
> parsing argument list. E.g. "demo pwelch" would fail:
> 
> pwelch example 1: failed
> error: pwelch: arg 3 (overlap) must be real from 0 to 0.950000
> 
> apparently on the line pwelch(x,[],[],[],[],0.9);
> 
> It also gets easily confused when pwelch called with two channels, e.g.:
> 
> pwelch(x, y, [], [], [], 1000, [], [], "mean", [], "trans")
> error: pwelch: can give confidence interval for x power spectrum only
> error: evaluating if command near line 507, column 3
> error: evaluating if command near line 279, column 1
> error: called from `pwelch' in file
> `/usr/share/octave/site/m/octave-forge/signal/pwelch.m'
> 
> In fact it appears that i cannot specify any parameter between
> "sampling frequency" and "result" parameters...
> 
> Sorry, I do not have a patch, the source of pwelch.m is too
> unwieldy for me.
> 
> Sincerely,
> 
> Dmitri.
> 

The demo ran through fine for me with the 20070528 release of Octave forge
and 2.9.12. Could you give another simpler example that fails?

D.

-- 
View this message in context: http://www.nabble.com/pwelch-bugs-tf3854366.html#a10920639
Sent from the octave-dev mailing list archive at Nabble.com.

Re: [OctDev] pwelch bugs

[Peter V. L. <pv...@me...>]

Hi Dmitri,

On Fri, Jun 01, 2007 at 03:38:06PM -0500, Dmitri A. Sergatskov wrote:
> The new version of pwelch (by  Peter V. Lanspeary) seems to have problem
> parsing argument list. E.g. "demo pwelch" would fail:
> 
> pwelch example 1: failed
> error: pwelch: arg 3 (overlap) must be real from 0 to 0.950000
> 
> apparently on the line pwelch(x,[],[],[],[],0.9);

"demo pwelch" works for me.
"pwelch(x,[],[],[],[],0.9)" gives 90% confidence limits, as expected.
To obtain "error: pwelch: arg 3 (overlap) must be real from 0 to 0.950000"
I need to do something like "pwelch(x,[],0.99)".

> It also gets easily confused when pwelch called with two channels, e.g.:
> 
> pwelch(x, y, [], [], [], 1000, [], [], "mean", [], "trans")
> error: pwelch: can give confidence interval for x power spectrum only
> error: evaluating if command near line 507, column 3
> error: evaluating if command near line 279, column 1
> error: called from `pwelch' in file
> `/usr/share/octave/site/m/octave-forge/signal/pwelch.m'
> 
> In fact it appears that i cannot specify any parameter between
> "sampling frequency" and "result" parameters...

I have added an explanation of your problem in the documentation string:

%%  [spectra,Pxx_ci,freq] = pwelch(x,window,overlap,Nfft,Fs,conf,
%%                                 range,plot_type,detrend,sloppy)
%%  [spectra,Pxx_ci,freq] = pwelch(x,y,window,overlap,Nfft,Fs,conf,
%%                                 range,plot_type,detrend,sloppy,results)
%%     Estimates confidence intervals for the spectral density.
%%     See Hint (7) below for compatibility options.  Confidence level "conf"
%%     is the 6th or 7th numeric argument.  If "results" control-string 
%%     arguments are used, one of them must be "power" when the "conf"
%%     argument is present; pwelch can estimate confidence intervals only for
%%     the power spectrum of the "x" data.  It does not know how to estimate
%%     confidence intervals of the cross-power spectrum, transfer function or
%%     coherence; if you can suggest a good method, please send a bug report.

and I have committed the changes.

pwelch(x, y, [], [], [], 1000, [], [], "mean", [], "trans", "power" )
returns with 
 a) the 95% confidence interval for the power spectral density of "x", and 
 b) the transfer function without a confidence interval.

> Sorry, I do not have a patch, the source of pwelch.m is too
> unwieldy for me.

That's fine.  I think that what you apparently found was normal
(expected) behaviour, not a bug.

Thanks,
Peter

Re: [OctDev] pwelch bugs

[Dmitri A. S. <das...@gm...>]

On 6/4/07, Peter Vernon Lanspeary <pv...@me...> wrote:
> Hi Dmitri,
>
> On Fri, Jun 01, 2007 at 03:38:06PM -0500, Dmitri A. Sergatskov wrote:
> > The new version of pwelch (by  Peter V. Lanspeary) seems to have problem
> > parsing argument list. E.g. "demo pwelch" would fail:
> >
> > pwelch example 1: failed
> > error: pwelch: arg 3 (overlap) must be real from 0 to 0.950000
> >
> > apparently on the line pwelch(x,[],[],[],[],0.9);
>
> "demo pwelch" works for me.
> "pwelch(x,[],[],[],[],0.9)" gives 90% confidence limits, as expected.
> To obtain "error: pwelch: arg 3 (overlap) must be real from 0 to 0.950000"
> I need to do something like "pwelch(x,[],0.99)".
>

I guess I forgot to CC my reply to David to the list. I was having the problem
because, apparently,  I had two "pwelch.m" in my path. For some reason
octave picked up demo code from the wrong file (while "which pwelch",
"help pwelch" , "type pwelch", etc... all looked at the correct one).

> > It also gets easily confused when pwelch called with two channels, e.g.:
> >
> > pwelch(x, y, [], [], [], 1000, [], [], "mean", [], "trans")
> > error: pwelch: can give confidence interval for x power spectrum only
> > error: evaluating if command near line 507, column 3
> > error: evaluating if command near line 279, column 1
> > error: called from `pwelch' in file
> > `/usr/share/octave/site/m/octave-forge/signal/pwelch.m'
> >
> > In fact it appears that i cannot specify any parameter between
> > "sampling frequency" and "result" parameters...
>
> I have added an explanation of your problem in the documentation string:
>
> %%  [spectra,Pxx_ci,freq] = pwelch(x,window,overlap,Nfft,Fs,conf,
> %%                                 range,plot_type,detrend,sloppy)
> %%  [spectra,Pxx_ci,freq] = pwelch(x,y,window,overlap,Nfft,Fs,conf,
> %%                                 range,plot_type,detrend,sloppy,results)
> %%     Estimates confidence intervals for the spectral density.
> %%     See Hint (7) below for compatibility options.  Confidence level "conf"
> %%     is the 6th or 7th numeric argument.  If "results" control-string
> %%     arguments are used, one of them must be "power" when the "conf"
> %%     argument is present; pwelch can estimate confidence intervals only for
> %%     the power spectrum of the "x" data.  It does not know how to estimate
> %%     confidence intervals of the cross-power spectrum, transfer function or
> %%     coherence; if you can suggest a good method, please send a bug report.
>
> and I have committed the changes.
>
> pwelch(x, y, [], [], [], 1000, [], [], "mean", [], "trans", "power" )
> returns with
>  a) the 95% confidence interval for the power spectral density of "x", and
>  b) the transfer function without a confidence interval.
>
> > Sorry, I do not have a patch, the source of pwelch.m is too
> > unwieldy for me.
>
> That's fine.  I think that what you apparently found was normal
> (expected) behaviour, not a bug.
>

The thing that confused me here is that if i used [] in place for
string arguments,
I get an error. E.g:

octave:1> x=rand(1,1024); y=x;
octave:2> pwelch(x,y,[],[],[],100,[],[],[],"mean", [], "trans")
error: pwelch: can give confidence interval for x power spectrum only
error: evaluating if command near line 507, column 3
error: evaluating if command near line 279, column 1
error: called from `pwelch' in file
`/usr/share/octave/site/m/octave-forge/signal/pwelch.m'
octave:2>

If I call it instead like that:

pwelch(x,y,[],[],[],100,"mean", "trans")

it would work fine. I guess the issue is that pwelch for (say) "trans"
would fail even with the defaul value for "range". In my opinion
it would be better to issue a warning:

"Do not know how to estimate the confifdence intervals of the cross-power
spectrum. Parameter ignored"

and, if possible, have default value for "conf" is [] (?!)

> Thanks,
> Peter
>

Thank you for your help.

Sincerely,

Dmitri.
--

[OctDev] pwelch: confidence interval for power spectrum

[Peter V. L. <pv...@me...>]

> >%%  [spectra,Pxx_ci,freq] = pwelch(x,window,overlap,Nfft,Fs,conf,
> >%%                                 range,plot_type,detrend,sloppy)
> >%%  [spectra,Pxx_ci,freq] = pwelch(x,y,window,overlap,Nfft,Fs,conf,
> >%%                                 range,plot_type,detrend,sloppy,results)
> >%%     Estimates confidence intervals for the spectral density.
> >%%     See Hint (7) below for compatibility options. Confidence level "conf"
> >%%     is the 6th or 7th numeric argument.  If "results" control-string
> >%%     arguments are used, one of them must be "power" when the "conf"
> >%%     argument is present; pwelch can estimate confidence intervals only 
> >%%     for the power spectrum of the "x" data.  It does not know how to 
> >%%     estimate confidence intervals of the cross-power spectrum, transfer
> >%%     function or coherence; if you can suggest a good method, please send
> >%%     a bug report.
> >
> If I call it instead like that:
> pwelch(x,y,[],[],[],100,"mean", "trans")
> it would work fine. I guess the issue is that pwelch for (say) "trans"
> would fail even with the defaul value for "range". In my opinion
> it would be better to issue a warning:
> "Do not know how to estimate the confifdence intervals of the cross-power
> spectrum. Parameter ignored"
> and, if possible, have default value for "conf" is [] (?!)

I think the safest policy is to stop with an error if pwelch cannot
estimate the confidence intervals.  Then, at least the user/programmer
is forced to fix the call to pwelch.  If pwelch just continues with a
warning, it provides "freq" in place of the expected "Pxx_ci" return
value, which leads to yet more problems further downstream.
However, I'm happy to listen other opinions.

If a "conf" argument is not provided in a call to pwelch, pwelch makes
no attempt to estimate the confidence intervals. This maintains
compatibility with Matlab.  If an empty "conf" arguement is given, the
default confidence level is 95%.

Also, it might be worth pointing out that the method used for estimating
confidence interval in pwelch is not entirely rigorous.  The method
incorrectly assumes that the periodogram has a Gaussian probability
distribution.  The distribution is exponential, so
it is wrong to use erfinv() in the calculation,
   half_width_of_confidence_interval =
      erfinv(confidence_level) * standard_deviation.
However, if the number of FFTs is large, pwelch may still give a useful
approximation to the confidence interval of the power spectrum.  I observe
that, in Matlab R11 pwelch, the feature which estimates confidence interval
was discontinued in the R12 version.

The reasons for not attempting to estimate confidence intervals for
cross-power spectrum, transfer function or coherence is that I haven't
seen the necessary mathematical theory.  Part of the problem is that the
cross-power spectrum is complex (rather than real).  I'll leave the
interested reader to think about the consequences ...
and arrive at a solution. :)

Thanks,
Peter