Hi,
From: Ariel Rokem <arokem@...453...>
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.
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
one_sided_spectrum_normalization.diff (1.05 KB)
···
On Fri, Feb 5, 2010 at 11:53 PM, Ludwig Schwardt <ludwig.schwardt@…716…> wrote:
Hi,
From: Ariel Rokem <arokem@…453…>
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