Hi - thanks Ludwig. I don't think that a major reworking of the logic of the function is needed. Simply replacing the line you mentioned with:
Pxy *= 1 / (np.abs(windowVals)**2).sum()
Pxy[1:-1] *= scaling_factor
if scale_by_freq:
Pxy[[0,-1]] /= Fs
seems to solve the 0 and NFFT/2 problem for the 'onesided' case (that is,
the result is very similar to the Matlab result, patch attached) and still
seems to work with twosided (because then the scaling factor is equal to 1
and it doesn't really do anything, except scale by frequency, if that is
required). In the 'twosided' case, the output in matlab is fftshift-ed
relative to the mlab output, but that's not crucial.
What does become more apparent when I do that is that in frequency bands in
which the power is rather small, the ratio discrepancies between the mlab
result and the matlab result can be rather large, on the order of a factor
of 2-2.5, even when the differences are tiny. Similarly, when the power is
rather large, there can be non-negligible differences between the two
results. Is there anything to do about that?
Cheers,
Ariel
On Fri, Feb 5, 2010 at 11:53 PM, Ludwig Schwardt
<ludwig.schwa...@gmail.com>wrote:
> Hi,
>
> > From: Ariel Rokem <aro...@berkeley.edu>
> > However - two elements are off by a factor of approximately 2 - the
> > very first element and the very last. ... Does anyone have any idea
> > why this would be the case?
>
> >From a quick look at the mlab code, it looks like a bug in
> mlab._spectral_helper.
>
> The default spectrum is 'onesided' (same as for Matlab's cpsd). A
> single-sided spectrum of a real signal has double the magnitude of a
> double-sided spectrum, *except* at the origin (frequency index n = 0)
> and Nyquist frequency (n = NFFT / 2), where it is the *same* as the
> double-sided one [1]_.
>
> In the mlab code, all the spectral values are simply scaled by a
> factor of 2 (among other factors) in this line:
>
> # Scale the spectrum by the norm of the window to compensate for
> # windowing loss; see Bendat & Piersol Sec 11.5.2. Also include
> # scaling factors for one-sided densities and dividing by the sampling
> # frequency, if desired.
> Pxy *= scaling_factor / (np.abs(windowVals)**2).sum()
>
> This should be easy to fix (although the function probably needs a
> little rework).
>
> Regards,
> Ludwig
>
> Quick reference from my bookshelf:
> ---------------------------------------------------
> .. [1] W. L. Briggs, V. E. Henson, "The DFT: An Owner's Manual for the
> Discrete Fourier Transform," Section 1.3, Problem 6 (a), p. 13.
>
>
> ------------------------------------------------------------------------------
> The Planet: dedicated and managed hosting, cloud storage, colocation
> Stay online with enterprise data centers and the best network in the
> business
> Choose flexible plans and management services without long-term contracts
> Personal 24x7 support from experience hosting pros just a phone call away.
> http://p.sf.net/sfu/theplanet-com
> _______________________________________________
> Matplotlib-devel mailing list
> Matplotlib-devel@lists.sourceforge.net
> https://lists.sourceforge.net/lists/listinfo/matplotlib-devel
>
--
Ariel Rokem
Helen Wills Neuroscience Institute
University of California, Berkeley
http://argentum.ucbso.berkeley.edu/ariel
one_sided_spectrum_normalization.diff
Description: Binary data
------------------------------------------------------------------------------ The Planet: dedicated and managed hosting, cloud storage, colocation Stay online with enterprise data centers and the best network in the business Choose flexible plans and management services without long-term contracts Personal 24x7 support from experience hosting pros just a phone call away. http://p.sf.net/sfu/theplanet-com
_______________________________________________ Matplotlib-devel mailing list Matplotlib-devel@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/matplotlib-devel
