Re: [matplotlib-devel] Comparison of mlab.csd and Matlab's cpsd

 Re: [matplotlib-devel] Comparison of mlab.csd and Matlab's cpsd From: Ludwig Schwardt - 2010-02-06 07:54:01 ```Hi, > From: Ariel Rokem > 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. ```

 Re: [matplotlib-devel] Comparison of mlab.csd and Matlab's cpsd From: Ludwig Schwardt - 2010-02-06 07:54:01 ```Hi, > From: Ariel Rokem > 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. ```
 Re: [matplotlib-devel] Comparison of mlab.csd and Matlab's cpsd From: Ariel Rokem - 2010-02-08 03:25:57 Attachments: Message as HTML ```Hi Ludwig (responding also to list), I don't think that the cause of the discrepancy is because of the hamming/hanning window difference. I do set the window in the matlab part to also be a hanning window of length nfft. Cheers, Ariel On Sun, Feb 7, 2010 at 7:45 AM, Ludwig Schwardt wrote: > Hi, > > On Sun, Feb 7, 2010 at 12:46 AM, Ariel Rokem wrote: > > 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 > > I agree. I was hoping the first two lines above would be sufficient. > Then I saw scaling_factor also included Fs if scaling by frequency, > and Pxy is of shape (numFreqs,n), i.e. not a straightforward 1-D > array, which caused me to reserve my judgment a little bit... :-) > > > 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? > > Could this be because Matlab uses a Hamming window by default, while > mlab uses a *Hanning* window by default? Very similar-sounding names, > and also very similar windows numerically (but not exactly the > same)... > > Ludwig > -- Ariel Rokem Helen Wills Neuroscience Institute University of California, Berkeley http://argentum.ucbso.berkeley.edu/ariel ```
 Re: [matplotlib-devel] Comparison of mlab.csd and Matlab's cpsd From: Ariel Rokem - 2010-02-08 17:32:30 Attachments: Message as HTML ```Hi - the maximal difference in this case is of about 6 (units?), which is an approximately 50% difference. This is in one point in the spectrum which has a relatively small value - the maximal peaks in the spectra are on the order of 1200, so in the grand scheme of things, not that horrible. Other large differences are on the order of 4, which is approximately 5% or less of those points. Thanks - Ariel PS for some reason, matplotlib-devel will not get your email unless you hit "reply all", so they have been getting my emails, but not yours. I am not sure whether that is what you intended, so I thought I would mention it. On Mon, Feb 8, 2010 at 1:24 AM, Ludwig Schwardt wrote: > Hi, > > On Mon, Feb 8, 2010 at 5:25 AM, Ariel Rokem wrote: > > I don't think that the cause of the discrepancy is because of the > > hamming/hanning window difference. I do set the window in the matlab part > to > > also be a hanning window of length nfft. > > I suspected you gave the same window to both, but I was just > checking... :-) To find those smaller discrepancies might be a bit > harder then, requiring a careful comparison of the various steps > involved. Just for interest's sake, how big are the differences we are > talking about? > > Regards, > Ludwig > -- Ariel Rokem Helen Wills Neuroscience Institute University of California, Berkeley http://argentum.ucbso.berkeley.edu/ariel ```
 Re: [matplotlib-devel] Comparison of mlab.csd and Matlab's cpsd From: Ariel Rokem - 2010-02-10 00:59:15 Attachments: Message as HTML ```Hi everyone, Mystery solved (I think)! It turns out that matlab has two functions, one called "hann" and the other called "hanning". They both purport to produce the hanning window, but they produce two (slightly) different things. Only one of these (the one given by "hann") actually produces the hanning window as it is defined in the standard reference on this (Harris 1978. And also Wikipedia ;D). When using the "hann" function in Matlab in this little experiment, the difference between the python and the matlab results, is on the order of %10e-6. I think I am willing to live with that. The only mystery that remains (I think) is - what is Matlab's "hanning" supposed to be doing? Are there two different windows: hann and hanning (not "hamming" - I know that's a different one)? Cheers, Ariel On Mon, Feb 8, 2010 at 9:32 AM, Ariel Rokem wrote: > Hi - > > the maximal difference in this case is of about 6 (units?), which is an > approximately 50% difference. This is in one point in the spectrum which has > a relatively small value - the maximal peaks in the spectra are on the order > of 1200, so in the grand scheme of things, not that horrible. Other large > differences are on the order of 4, which is approximately 5% or less of > those points. > > Thanks - Ariel > > PS for some reason, matplotlib-devel will not get your email unless you hit > "reply all", so they have been getting my emails, but not yours. I am not > sure whether that is what you intended, so I thought I would mention it. > > > On Mon, Feb 8, 2010 at 1:24 AM, Ludwig Schwardt > wrote: > >> Hi, >> >> On Mon, Feb 8, 2010 at 5:25 AM, Ariel Rokem wrote: >> > I don't think that the cause of the discrepancy is because of the >> > hamming/hanning window difference. I do set the window in the matlab >> part to >> > also be a hanning window of length nfft. >> >> I suspected you gave the same window to both, but I was just >> checking... :-) To find those smaller discrepancies might be a bit >> harder then, requiring a careful comparison of the various steps >> involved. Just for interest's sake, how big are the differences we are >> talking about? >> >> Regards, >> Ludwig >> > > > > -- > Ariel Rokem > Helen Wills Neuroscience Institute > University of California, Berkeley > http://argentum.ucbso.berkeley.edu/ariel > -- Ariel Rokem Helen Wills Neuroscience Institute University of California, Berkeley http://argentum.ucbso.berkeley.edu/ariel ```
 Re: [matplotlib-devel] Comparison of mlab.csd and Matlab's cpsd From: Ariel Rokem - 2010-02-10 01:22:48 Attachments: Message as HTML ```Hi - mea culpa - as pointed out by a colleague (and documented in the help document for "hanning" in Matlab), this function produces the same window, but without the first and last zero-weighted samples. Cheers -- Ariel On Tue, Feb 9, 2010 at 4:58 PM, Ariel Rokem wrote: > Hi everyone, > > Mystery solved (I think)! It turns out that matlab has two functions, one > called "hann" and the other called "hanning". They both purport to produce > the hanning window, but they produce two (slightly) different things. Only > one of these (the one given by "hann") actually produces the hanning window > as it is defined in the standard reference on this (Harris 1978. And also > Wikipedia ;D). When using the "hann" function in Matlab in this little > experiment, the difference between the python and the matlab results, is on > the order of %10e-6. I think I am willing to live with that. > > The only mystery that remains (I think) is - what is Matlab's "hanning" > supposed to be doing? Are there two different windows: hann and hanning (not > "hamming" - I know that's a different one)? > > Cheers, > > Ariel > > > On Mon, Feb 8, 2010 at 9:32 AM, Ariel Rokem wrote: > >> Hi - >> >> the maximal difference in this case is of about 6 (units?), which is an >> approximately 50% difference. This is in one point in the spectrum which has >> a relatively small value - the maximal peaks in the spectra are on the order >> of 1200, so in the grand scheme of things, not that horrible. Other large >> differences are on the order of 4, which is approximately 5% or less of >> those points. >> >> Thanks - Ariel >> >> PS for some reason, matplotlib-devel will not get your email unless you >> hit "reply all", so they have been getting my emails, but not yours. I am >> not sure whether that is what you intended, so I thought I would mention it. >> >> >> >> On Mon, Feb 8, 2010 at 1:24 AM, Ludwig Schwardt < >> ludwig.schwardt@...> wrote: >> >>> Hi, >>> >>> On Mon, Feb 8, 2010 at 5:25 AM, Ariel Rokem wrote: >>> > I don't think that the cause of the discrepancy is because of the >>> > hamming/hanning window difference. I do set the window in the matlab >>> part to >>> > also be a hanning window of length nfft. >>> >>> I suspected you gave the same window to both, but I was just >>> checking... :-) To find those smaller discrepancies might be a bit >>> harder then, requiring a careful comparison of the various steps >>> involved. Just for interest's sake, how big are the differences we are >>> talking about? >>> >>> Regards, >>> Ludwig >>> >> >> >> >> -- >> Ariel Rokem >> Helen Wills Neuroscience Institute >> University of California, Berkeley >> http://argentum.ucbso.berkeley.edu/ariel >> > > > > -- > Ariel Rokem > Helen Wills Neuroscience Institute > University of California, Berkeley > http://argentum.ucbso.berkeley.edu/ariel > -- Ariel Rokem Helen Wills Neuroscience Institute University of California, Berkeley http://argentum.ucbso.berkeley.edu/ariel ```