Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP

Loading…

ENH: signal: spectral: Add periodogram and welch functions #373

Merged
merged 17 commits into from

4 participants

Eric Moore Ralf Gommers Pauli Virtanen Thomas Hisch
Eric Moore
Collaborator

I think scipy.signal should offer functions to do spectral estimation. We have lombscargle, but no conventional DFT based spectral estimators. I've coded up a simple (modified) periodogram and Welch's method in this pull request. There are other methods that might be nice to include (multitaper, burg's method, etc) but I think that these make a good start.

The functions as given here don't allow cross spectra, but that could be easily remedied given what I have. (Although I think separate functions are warranted.)

There are no tests at present, but I'll add them soon. I suspect there might be discussion of the signature and exposed functionality.

I'll send a email to the discussion list, too.

ewmoore added some commits
Eric Moore ewmoore ENH: signal: spectral: Add periodogram and welch functions 0ae062a
Eric Moore ewmoore BUG: signal: spectral: fix various bugs in periodogram and welch.
1. various off by 1 issues in dealing with the output of fftpack.rfft.
2. the DC and nfft/2 points shouldn't be scaled by 2 for 'onesided'
3. error in computing frequency axis in 'onesided'

The output of welch exactly matches matplotlib.mlab.psd for both even
and odd length  one-sided spectra and even length two-sided spectra of
real data.  For odd length two-sided spectra of real data, welch matches
the spectrum but matplotlib.mlab.psd return an incorrect frequcny axis.
9bfb88b
Eric Moore ewmoore BUG: signal: spectral: Fix bug introduced in last commit
We want fftbins=True in signal.get_window.  matplotlib.mlab.psd
effectively uses fftbin=False since this option doesn't exist in numpy.
Swapped this in to do some testing.
0985cc0
Eric Moore ewmoore ENH: signal: spectral: Support detrending in welch. 1ffba7f
Eric Moore ewmoore MAINT: signal: update binary attrib on _spectral.c
This marking was lost when the file was renamed.
2bee282
Eric Moore ewmoore BUG: signal: welch: call detrend instead of doing it by hand.
To support detending via a polynomial, we could extend detrend.
44dadea
Eric Moore ewmoore TST: signal: add test for welch
Also, fix bug in welch.
47ebb63
Eric Moore ewmoore TST: signal: add tests for periodogram b171354
Eric Moore
Collaborator

I think that the welch and periodogram functions are complete now.

My welch supports everything that matplotlib's psd supports, with the exception of allowing padded FFTs. I think that matplotlib makes the wrong choice for scaling in that case. It preserves the absolute magnitude of the PSD instead of the area under the PSD. If someone speaks up about wanting this I'll add it.

scipy/signal/spectral.py
@@ -0,0 +1,384 @@
+"""Tools for spectral analysis.
+"""
+
+import numpy as np
+from scipy import fftpack
+import signaltools
+from windows import get_window
+from _spectral import *
Ralf Gommers Owner
rgommers added a note

from _spectral import lombscargle.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
scipy/signal/spectral.py
@@ -0,0 +1,384 @@
+"""Tools for spectral analysis.
+"""
+
+import numpy as np
+from scipy import fftpack
+import signaltools
+from windows import get_window
+from _spectral import *
+
+__all__ = ['periodogram', 'welch', 'lombscargle']
+
+
+def periodogram(x, fs=1.0, window=None, nfft=None, sides='default', scaling='density', axis=-1):
Ralf Gommers Owner
rgommers added a note

line break, 80 char lines

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Ralf Gommers
Owner

Some minor things about docstring formatting:

  • no blank lines in Parameters and Returns sections
  • type of string params should be str, not string
  • int or None should simply be int if None is already the default
Ralf Gommers
Owner

Doesn't work for py3k. In tools/py3tool.py, replace spectral by _spectral.

scipy/signal/spectral.py
@@ -0,0 +1,384 @@
+"""Tools for spectral analysis.
+"""
+
+import numpy as np
+from scipy import fftpack
+import signaltools
+from windows import get_window
+from _spectral import *
+
+__all__ = ['periodogram', 'welch', 'lombscargle']
+
+
+def periodogram(x, fs=1.0, window=None, nfft=None, sides='default', scaling='density', axis=-1):
+ """
+ Estimate power spectal density using a periodogram.
Ralf Gommers Owner
rgommers added a note

typo, "spectral"

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
scipy/signal/spectral.py
((12 lines not shown))
+
+def periodogram(x, fs=1.0, window=None, nfft=None, sides='default', scaling='density', axis=-1):
+ """
+ Estimate power spectal density using a periodogram.
+
+ Parameters
+ ----------
+ x : array_like
+ Time series of measurement values
+
+ fs : float, optional
+ Sampling frequency of the `x` time series in units of Hz. Defaults
+ to 1.0.
+
+ window : string or tuple of string and parameter values or array_like
+ Desired window to use. See scipy.signal.get_window for a list of
Ralf Gommers Owner
rgommers added a note

scipy.signal.get_window can be replaced by `get_window` with single backticks.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
scipy/signal/spectral.py
((14 lines not shown))
+ """
+ Estimate power spectal density using a periodogram.
+
+ Parameters
+ ----------
+ x : array_like
+ Time series of measurement values
+
+ fs : float, optional
+ Sampling frequency of the `x` time series in units of Hz. Defaults
+ to 1.0.
+
+ window : string or tuple of string and parameter values or array_like
+ Desired window to use. See scipy.signal.get_window for a list of
+ windows and required parameters. If `window` is an array it will
+ be used directly as the window. Defaults to 'None'; equivalent to
Ralf Gommers Owner
rgommers added a note

'None' --> None

Same in nfft description below.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Ralf Gommers rgommers commented on the diff
scipy/signal/spectral.py
((52 lines not shown))
+ Returns
+ -------
+ f : ndarray
+ Array of sample frequencies.
+
+ Pxx : ndarray
+ Power spectral density or power spectrum of `x`.
+
+ See Also
+ --------
+ welch: Estimate power spectral density using Welch's method
+ lombscargle: Lomb-Scargle periodogram for unevenly sampled data
+
+ Examples
+ --------
+ >>> from scipy import signal
Ralf Gommers Owner
rgommers added a note

Needs also import matplotlib.pyplot as plt.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Ralf Gommers rgommers commented on the diff
scipy/signal/spectral.py
((71 lines not shown))
+
+ >>> fs = 10e3
+ >>> N = 1e5
+ >>> amp = 2*np.sqrt(2)
+ >>> freq = 1234.0
+ >>> noise_power = 0.001 * fs / 2
+ >>> time = np.arange(N) / fs
+ >>> x = amp*np.sin(2*np.pi*freq*time)
+ >>> x += np.random.normal(scale=np.sqrt(noise_power), size=time.shape)
+
+ Compute and plot the power spectral density.
+
+ >>> f, Pxx_den = signal.periodogram(x, fs)
+ >>> plt.semilogy(f, Pxx_den)
+ >>> plt.xlabel('frequency [Hz]')
+ >>> plt.ylabel('PSD [V/sqrt(Hz)]')
Ralf Gommers Owner
rgommers added a note

Needs plt.show() here and in second example below.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
scipy/signal/spectral.py
((101 lines not shown))
+
+ The peak height in the power spectrum is an estimate of the RMS amplitude.
+
+ >>> np.sqrt(Pxx_spec.max())
+ 2.0077340678640727
+ """
+ x = np.asarray(x)
+
+ if axis != -1:
+ x = np.rollaxis(x, axis, len(x.shape))
+
+ if nfft is None:
+ nfft = x.shape[-1]
+
+ if window is not None:
+ if type(window) is str or type(window) is tuple:
Ralf Gommers Owner
rgommers added a note

Checking if something is a string is normally done with isinstance(window, basestring), this doesn't work for unicode strings.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
scipy/signal/spectral.py
((121 lines not shown))
+ raise ValueError('window must be 1-D')
+ elif win.shape[0] != x.shape[-1]:
+ raise ValueError('window length must match data length.')
+ x = win*x
+ s1 = win.sum()**2
+ s2 = (win*win).sum()
+ else:
+ s1 = nfft**2
+ s2 = nfft
+
+ if scaling == 'density':
+ scale = 1.0 / (fs * s2)
+ elif scaling == 'spectrum':
+ scale = 1.0 / s1
+ else:
+ raise ValueError('Unknown scaling "{0}".'.format(scaling))
Ralf Gommers Owner
rgommers added a note

str.format isn't valid Python 2.4 syntax

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
scipy/signal/spectral.py
((123 lines not shown))
+ raise ValueError('window length must match data length.')
+ x = win*x
+ s1 = win.sum()**2
+ s2 = (win*win).sum()
+ else:
+ s1 = nfft**2
+ s2 = nfft
+
+ if scaling == 'density':
+ scale = 1.0 / (fs * s2)
+ elif scaling == 'spectrum':
+ scale = 1.0 / s1
+ else:
+ raise ValueError('Unknown scaling "{0}".'.format(scaling))
+
+ if np.isrealobj(x) and (sides == 'default' or sides == 'onesided'):
Ralf Gommers Owner
rgommers added a note

The sides description in the docstring has a capital letter for Default. Put a line sides = sides.lower() somewhere at the top?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
scipy/signal/spectral.py
((130 lines not shown))
+
+ if scaling == 'density':
+ scale = 1.0 / (fs * s2)
+ elif scaling == 'spectrum':
+ scale = 1.0 / s1
+ else:
+ raise ValueError('Unknown scaling "{0}".'.format(scaling))
+
+ if np.isrealobj(x) and (sides == 'default' or sides == 'onesided'):
+ x = fftpack.rfft(x, nfft)
+ outshape = list(x.shape)
+ if nfft % 2 == 0: # even
+ outshape[-1] = nfft/2+1
+ Pxx = np.empty(outshape, x.dtype)
+ Pxx[..., (0,-1)] = x[..., (0,-1)]**2
+ Pxx[..., 1:-1] = x[..., 1:-1:2]**2 + x[..., 2::2]**2
Ralf Gommers Owner
rgommers added a note

This indexing notation is a little hard to decipher, could use a comment.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
scipy/signal/spectral.py
((147 lines not shown))
+ outshape[-1] = (nfft+1)/2
+ Pxx = np.empty(outshape, x.dtype)
+ Pxx[..., 0] = x[..., 0]**2
+ Pxx[..., 1:] = x[..., 1::2]**2 + x[..., 2::2]**2
+ Pxx[..., 1:-1] *= 2*scale
+ Pxx[..., (0,-1)] *= scale
+ f = np.arange(Pxx.shape[-1])*(fs/nfft)
+ elif (np.iscomplexobj(x) and sides == 'default') or sides == 'twosided':
+ x = fftpack.fft(x, nfft)
+ Pxx = (x * x.conj()).real
+ f = fftpack.fftfreq(nfft, 1.0/fs)
+ Pxx *= scale
+ elif np.iscomplexobj(x) and sides == 'onesided':
+ raise ValueError('Cannot compute onsided spectrum of complex data')
+ else:
+ raise ValueError('Unknown sides: "{0}".'.format(sides))
Ralf Gommers Owner
rgommers added a note

same as above, 2.4 syntax error

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
scipy/signal/spectral.py
((150 lines not shown))
+ Pxx[..., 1:] = x[..., 1::2]**2 + x[..., 2::2]**2
+ Pxx[..., 1:-1] *= 2*scale
+ Pxx[..., (0,-1)] *= scale
+ f = np.arange(Pxx.shape[-1])*(fs/nfft)
+ elif (np.iscomplexobj(x) and sides == 'default') or sides == 'twosided':
+ x = fftpack.fft(x, nfft)
+ Pxx = (x * x.conj()).real
+ f = fftpack.fftfreq(nfft, 1.0/fs)
+ Pxx *= scale
+ elif np.iscomplexobj(x) and sides == 'onesided':
+ raise ValueError('Cannot compute onsided spectrum of complex data')
+ else:
+ raise ValueError('Unknown sides: "{0}".'.format(sides))
+
+ if axis != -1:
+ Pxx = np.rollaxis(Pxx, len(Pxx.shape)-1, axis)
Ralf Gommers Owner
rgommers added a note

second argument can simply be -1.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
scipy/signal/spectral.py
((159 lines not shown))
+ elif np.iscomplexobj(x) and sides == 'onesided':
+ raise ValueError('Cannot compute onsided spectrum of complex data')
+ else:
+ raise ValueError('Unknown sides: "{0}".'.format(sides))
+
+ if axis != -1:
+ Pxx = np.rollaxis(Pxx, len(Pxx.shape)-1, axis)
+ return f, Pxx
+
+
+def welch(x, fs=1.0, window='hanning', nfft=256, noverlap=None, detrend='constant', sides='default', scaling='density', axis=-1):
+ """
+ Estimate power spectral density using Welch's method.
+
+ Welch's method [1]_ computes an estimate of the power spectral desnsity
+ by dividing the data into overlaping segments, computing a modified
Ralf Gommers Owner
rgommers added a note

typo "overlaping"

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
scipy/signal/spectral.py
((306 lines not shown))
+ raise ValueError('window must be 1-D')
+ nfft = win.shape[0]
+
+ if scaling == 'density':
+ scale = 1.0 / (fs * (win*win).sum())
+ elif scaling == 'spectrum':
+ scale = 1.0 / win.sum()**2
+ else:
+ raise ValueError('Unknown scaling "{0}".'.format(scaling))
+
+ if noverlap is None:
+ noverlap = nfft // 2
+
+ if not hasattr(detrend, '__call__'):
+ detrend_func = lambda seg: signaltools.detrend(seg, type=detrend)
+ elif axis != -1:
Ralf Gommers Owner
rgommers added a note

I don't quite get the logic of elif here.

Eric Moore Collaborator
ewmoore added a note

First the if checks if the user has passed a function. If they have, and we aren't working with the last axis, I wanted their function to receive an array with a shape that they could predict. So if necessary, I'm wrapping their function in some calls to np.rollaxis. I suppose I don't have to wrap it and could just note in the docstring that the shape will have rolled axes, but that seems to unnecessarily expose an implementation detail.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
scipy/signal/spectral.py
((195 lines not shown))
+ noverlap: int or None, optional
+ Number of points to overlap between segments. If None, `noverlap`
+ = nfft / 2. Defaults to None.
+
+ detrend : string or function, optional
+ Specifies how to detend each segment. If `detrend` is a string,
+ it is passed as the `type` argument to scipy.signal.detrend.
+ If it is a function, it takes a segment and returns a detrended
+ segment. Defaults to 'constant'
+
+ sides : { 'default', 'onesided', 'twosided' }, optional
+ Selects which sides of the periodogram to return. 'default'
+ computes a one-sided spectrum for real data and a two-sided
+ spectrum for complex data. 'onesided' or 'twosided' will force
+ a one- or two-sided spectrum. 'onesided' cannot be specified with
+ complex data.
Ralf Gommers Owner
rgommers added a note

Isn't this more complicated than it has to be? Alternative would be to use something like

return_onesided : bool, optional
    If True, return a one-sided spectrum for real data. If False, return a two-sided spectrum.
    Note that for complex data, the two-sided spectrum is always returned.
Eric Moore Collaborator
ewmoore added a note

That's a much better idea. I'll change it in both welch and periodogram.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Ralf Gommers
Owner

Overall looks good to me.

It looks like there's quite some overlap between the two functions (checks on input values, asarray call, etc.). Would it make sense to factor out those parts into a common function?

scipy/signal/spectral.py
((210 lines not shown))
+ complex data.
+
+ scaling : { 'density', 'spectrum' }, optional
+ Selects between computing the power spectral density ('density')
+ where Pxx has units of V**2/Hz if x is measured in V and computing
+ the power spectrum ('spectrum') where Pxx has units of V**2 if x is
+ measured in V. Defaults to 'density'
+
+ axis : int, optional
+ Axis along which the periodogram is computed; the default is over
+ the last axis (i.e. ``axis=-1``).
+
+ Returns
+ -------
+ f : ndarray
+ Array of sample frequencies.
Thomas Hisch
thisch added a note

wrong indentation

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Eric Moore
Collaborator

You're quite right that there is a lot of overlap between the two functions. Welch is naturally implemented in terms of periodogram, but it's easier to do it the other way around and implement periodogram in terms of welch, which I have done. I've also added padded FFT's to welch since periodogram allowed that before. I think I've also addressed all of the comments.

Ralf Gommers
Owner

Looks like you missed my comment on py3k above. Other than that it looks good. Will look in a little more detail at the last changes around Christmas and merge it.

Eric Moore
Collaborator

Right now if you pass a float32 array in you'll get a float64 array out. I can fix that but I don't think scipy is particularly careful about that right now anyway.

Ralf Gommers
Owner

float64 out should be fine.

Ralf Gommers
Owner

In tools/py3tool.py there's still this line missing to make things work:

os.path.join('signal', 'spectral.py'),
Ralf Gommers
Owner

Fixed the py3k issue and some minor things in https://github.com/rgommers/scipy/tree/spectral. Could you review and pull that?

Ralf Gommers
Owner

The behavior when putting in empty arrays or arrays with few data points is odd, due to the zero-padding. For empty input, I propose to return empty output. If the length of the input is shorter than one segment, maybe issue a warning?

Ralf Gommers
Owner

There's no test using a non-default version of noverlap. Also, tests covering the first and third elif clause of periodogram would be good to add. Other than that, test coverage is quite complete.

Eric Moore ewmoore BUG: signal: spectral: various bugs and missing tests
periodogram chokes if passed a list
    --> added call to np.asarray
periodogram throws an unhelpful division by zero error if passed an empty array
    --> added check and quick return if x.shape == (0,)
    --> added test_empty_input
periodogram needs test for nfft < x.shape[axis]
    --> added test_short_nfft
periodogram needs test for nfft == x.shape[axis]
    --> added test_nfft_is_xshape
welch gives big output when passing an empty array
    --> added check and quick return with empty arrays if x.shape == (0,)
    --> added test_empty_input
welch's output is unexpected if x.shape[axis] < nperseg
    --> check and warn that nperseg > x.shape[axis], but use nperseg == x.shape[axis]
    --> added test_short_data
welch has similar issue if using custom window and short data
    --> added a check that raises if window longer than data
    --> added test_long_window
aa74948
Eric Moore
Collaborator

Still needs a check for noverlap < nperseg and test for the non-default noverlap.

Eric Moore ewmoore BUG: signal: welch: add check for invalid noverlap
Also add test for invalid noverlap, and test for nondefault noverlap.
aa85b1b
Eric Moore
Collaborator

Ralph, I think I've addressed all of the issues you've raised. Thanks for all of the pointers. I think that where this is now is much improved.

Ralf Gommers rgommers MAINT: in periodogram/welch, improve empty array handling and test wa…
…rnings.

Also add a few more tests for complete test coverage.
ba58448
Ralf Gommers
Owner

Thanks Eric, this looks quite good now.

I added one more commit at https://github.com/rgommers/scipy/tree/spectral to address some minor issues. Could you have a look at that? Then it's ready to merge I think.

Eric Moore
Collaborator

Look good to me.

Eric Moore
Collaborator

Oh, I think 'with' is a python 2.5+ feature. I guess this doesn't matter since we are dropping support for 2.4 anyway.

Ralf Gommers
Owner

Yep, I used with on purpose. Will change that in more places soon, the try: WarningManager.__enter__()/ finally:WarningManager.__exit__() is getting old.

Ralf Gommers rgommers merged commit e08beea into from
Ralf Gommers
Owner

OK, merged. Thanks again, nice enhancement!

Ralf Gommers
Owner

I can't reproduce this with 1.6.1, but with numpy master I can reproduce one failure. The problem is with using assert_array_almost_equal_nulp:

In [19]: p[:,0,0]
Out[19]: 
array([  6.66666667e+00,   7.94623303e+01,   3.68524270e+00,
         2.05371537e-01,   1.96112367e-02,   0.00000000e+00])

In [20]: p[:,1,0]
Out[20]: 
array([  6.66666667e+00,   7.94623303e+01,   3.68524270e+00,
         2.05371537e-01,   1.96112367e-02,   5.25907270e-32])

In [21]: np.testing.utils.nulp_diff(p[:,0,0], p[:,1,0])
Out[21]: 
array([  5.40000000e+01,   0.00000000e+00,   0.00000000e+00,
         9.00000000e+00,   5.00000000e+00,   4.13910830e+18])

Changing to assert_allclose should fix that. That is annoying behavior by the way. I vaguely remember running into this before; never use nulps if you can have zeros in your array.

Pauli Virtanen
Owner
pv commented

@rgommers I think ulps is supposed to work that way, it's a feature, not a bug :)
For instance in special function value comparisons, this comes useful.

Ralf Gommers
Owner

I can't reproduce the other failure, could you check that https://github.com/rgommers/scipy/tree/fix-spectral-failures fixes it?

Eric Moore
Collaborator

Sorry for disappearing. I haven't actually ever tested against anything but python 2.7 and the latest numpy release (1.6.2, I guess.) It looks to me as though @rgommers fix-spectral-failures tree ought to fix things.

Ralf Gommers
Owner

OK I pushed the fix in 9530149, let me know if anything else has to be done.

Clemens ClemensFMN referenced this pull request from a commit
Commit has since been removed from the repository and is no longer available.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Commits on Dec 4, 2012
  1. Eric Moore
Commits on Dec 5, 2012
  1. Eric Moore

    BUG: signal: spectral: fix various bugs in periodogram and welch.

    ewmoore authored
    1. various off by 1 issues in dealing with the output of fftpack.rfft.
    2. the DC and nfft/2 points shouldn't be scaled by 2 for 'onesided'
    3. error in computing frequency axis in 'onesided'
    
    The output of welch exactly matches matplotlib.mlab.psd for both even
    and odd length  one-sided spectra and even length two-sided spectra of
    real data.  For odd length two-sided spectra of real data, welch matches
    the spectrum but matplotlib.mlab.psd return an incorrect frequcny axis.
  2. Eric Moore

    BUG: signal: spectral: Fix bug introduced in last commit

    ewmoore authored
    We want fftbins=True in signal.get_window.  matplotlib.mlab.psd
    effectively uses fftbin=False since this option doesn't exist in numpy.
    Swapped this in to do some testing.
Commits on Dec 7, 2012
  1. Eric Moore
  2. Eric Moore

    MAINT: signal: update binary attrib on _spectral.c

    ewmoore authored
    This marking was lost when the file was renamed.
  3. Eric Moore

    BUG: signal: welch: call detrend instead of doing it by hand.

    ewmoore authored
    To support detending via a polynomial, we could extend detrend.
Commits on Dec 8, 2012
  1. Eric Moore

    TST: signal: add test for welch

    ewmoore authored
    Also, fix bug in welch.
Commits on Dec 9, 2012
  1. Eric Moore
Commits on Dec 14, 2012
  1. Eric Moore
Commits on Dec 19, 2012
  1. Eric Moore

    FIX: signal: spectral: implement periodogram in terms of welch

    ewmoore authored
    Also finish addressing comments.
  2. Eric Moore
Commits on Dec 21, 2012
  1. Eric Moore
Commits on Dec 26, 2012
  1. Ralf Gommers
  2. Ralf Gommers
Commits on Jan 4, 2013
  1. Eric Moore

    BUG: signal: spectral: various bugs and missing tests

    ewmoore authored
    periodogram chokes if passed a list
        --> added call to np.asarray
    periodogram throws an unhelpful division by zero error if passed an empty array
        --> added check and quick return if x.shape == (0,)
        --> added test_empty_input
    periodogram needs test for nfft < x.shape[axis]
        --> added test_short_nfft
    periodogram needs test for nfft == x.shape[axis]
        --> added test_nfft_is_xshape
    welch gives big output when passing an empty array
        --> added check and quick return with empty arrays if x.shape == (0,)
        --> added test_empty_input
    welch's output is unexpected if x.shape[axis] < nperseg
        --> check and warn that nperseg > x.shape[axis], but use nperseg == x.shape[axis]
        --> added test_short_data
    welch has similar issue if using custom window and short data
        --> added a check that raises if window longer than data
        --> added test_long_window
Commits on Jan 5, 2013
  1. Eric Moore

    BUG: signal: welch: add check for invalid noverlap

    ewmoore authored
    Also add test for invalid noverlap, and test for nondefault noverlap.
  2. Ralf Gommers

    MAINT: in periodogram/welch, improve empty array handling and test wa…

    rgommers authored
    …rnings.
    
    Also add a few more tests for complete test coverage.
This page is out of date. Refresh to see the latest.
2  .gitattributes
View
@@ -10,7 +10,7 @@ scipy/interpolate/interpnd.c binary
scipy/io/matlab/mio_utils.c binary
scipy/io/matlab/mio5_utils.c binary
scipy/io/matlab/streams.c binary
-scipy/signal/spectral.c binary
+scipy/signal/_spectral.c binary
scipy/sparse/csgraph/_min_spanning_tree.c binary
scipy/sparse/csgraph/_shortest_path.c binary
scipy/sparse/csgraph/_tools.c binary
2  scipy/signal/SConscript
View
@@ -13,7 +13,7 @@ env.NumpyPythonExtension('sigtools',
'firfilter.c', \
'medianfilter.c'])
-env.NumpyPythonExtension('spectral', source='spectral.c')
+env.NumpyPythonExtension('_spectral', source='_spectral.c')
env.NumpyPythonExtension('spline',
source = ['splinemodule.c', 'S_bspline_util.c',
2  scipy/signal/__init__.py
View
@@ -215,6 +215,8 @@
.. autosummary::
:toctree: generated/
+ periodogram -- Computes a (modified) periodogram
+ welch -- Compute a periodogram using Welch's method
lombscargle -- Computes the Lomb-Scargle periodogram
"""
154 scipy/signal/spectral.c → scipy/signal/_spectral.c
View
@@ -1,4 +1,4 @@
-/* Generated by Cython 0.17.1 on Mon Oct 29 21:26:04 2012 */
+/* Generated by Cython 0.17.2 on Mon Dec 3 19:25:15 2012 */
#define PY_SSIZE_T_CLEAN
#include "Python.h"
@@ -53,12 +53,15 @@
(PyErr_Format(PyExc_TypeError, \
"expected index value, got %.200s", Py_TYPE(o)->tp_name), \
(PyObject*)0))
- #define PyIndex_Check(o) (PyNumber_Check(o) && !PyFloat_Check(o) && !PyComplex_Check(o))
+ #define __Pyx_PyIndex_Check(o) (PyNumber_Check(o) && !PyFloat_Check(o) && \
+ !PyComplex_Check(o))
+ #define PyIndex_Check __Pyx_PyIndex_Check
#define PyErr_WarnEx(category, message, stacklevel) PyErr_Warn(category, message)
#define __PYX_BUILD_PY_SSIZE_T "i"
#else
#define __PYX_BUILD_PY_SSIZE_T "n"
#define CYTHON_FORMAT_SSIZE_T "z"
+ #define __Pyx_PyIndex_Check PyIndex_Check
#endif
#if PY_VERSION_HEX < 0x02060000
#define Py_REFCNT(ob) (((PyObject*)(ob))->ob_refcnt)
@@ -246,8 +249,8 @@
#define _USE_MATH_DEFINES
#endif
#include <math.h>
-#define __PYX_HAVE__scipy__signal__spectral
-#define __PYX_HAVE_API__scipy__signal__spectral
+#define __PYX_HAVE__scipy__signal___spectral
+#define __PYX_HAVE_API__scipy__signal___spectral
#include "stdio.h"
#include "stdlib.h"
#include "numpy/arrayobject.h"
@@ -360,7 +363,7 @@ static const char *__pyx_filename;
static const char *__pyx_f[] = {
- "spectral.pyx",
+ "_spectral.pyx",
"numpy.pxd",
"type.pxd",
};
@@ -1000,16 +1003,16 @@ static CYTHON_INLINE char *__pyx_f_5numpy__util_dtypestring(PyArray_Descr *, cha
/* Module declarations from 'cython' */
-/* Module declarations from 'scipy.signal.spectral' */
+/* Module declarations from 'scipy.signal._spectral' */
static __Pyx_TypeInfo __Pyx_TypeInfo_nn___pyx_t_5numpy_float64_t = { "float64_t", NULL, sizeof(__pyx_t_5numpy_float64_t), { 0 }, 0, 'R', 0, 0 };
-#define __Pyx_MODULE_NAME "scipy.signal.spectral"
-int __pyx_module_is_main_scipy__signal__spectral = 0;
+#define __Pyx_MODULE_NAME "scipy.signal._spectral"
+int __pyx_module_is_main_scipy__signal___spectral = 0;
-/* Implementation of 'scipy.signal.spectral' */
+/* Implementation of 'scipy.signal._spectral' */
static PyObject *__pyx_builtin_ValueError;
static PyObject *__pyx_builtin_range;
static PyObject *__pyx_builtin_RuntimeError;
-static PyObject *__pyx_pf_5scipy_6signal_8spectral_lombscargle(CYTHON_UNUSED PyObject *__pyx_self, PyArrayObject *__pyx_v_x, PyArrayObject *__pyx_v_y, PyArrayObject *__pyx_v_freqs); /* proto */
+static PyObject *__pyx_pf_5scipy_6signal_9_spectral_lombscargle(CYTHON_UNUSED PyObject *__pyx_self, PyArrayObject *__pyx_v_x, PyArrayObject *__pyx_v_y, PyArrayObject *__pyx_v_freqs); /* proto */
static int __pyx_pf_5numpy_7ndarray___getbuffer__(PyArrayObject *__pyx_v_self, Py_buffer *__pyx_v_info, int __pyx_v_flags); /* proto */
static void __pyx_pf_5numpy_7ndarray_2__releasebuffer__(PyArrayObject *__pyx_v_self, Py_buffer *__pyx_v_info); /* proto */
static char __pyx_k_1[] = "Input arrays do not have the same size.";
@@ -1020,8 +1023,8 @@ static char __pyx_k_9[] = "unknown dtype code in numpy.pxd (%d)";
static char __pyx_k_10[] = "Format string allocated too short, see comment in numpy.pxd";
static char __pyx_k_13[] = "Format string allocated too short.";
static char __pyx_k_15[] = "Tools for spectral analysis of unequally sampled signals.";
-static char __pyx_k_18[] = "/home/rgommers/Code/scipy/scipy/signal/spectral.pyx";
-static char __pyx_k_19[] = "scipy.signal.spectral";
+static char __pyx_k_18[] = "/home/ewm/scipy/scipy/signal/_spectral.pyx";
+static char __pyx_k_19[] = "scipy.signal._spectral";
static char __pyx_k_20[] = "lombscargle (line 19)";
static char __pyx_k_21[] = "\n lombscargle(x, y, freqs)\n\n Computes the Lomb-Scargle periodogram.\n \n The Lomb-Scargle periodogram was developed by Lomb [1]_ and further\n extended by Scargle [2]_ to find, and test the significance of weak\n periodic signals with uneven temporal sampling.\n\n The computed periodogram is unnormalized, it takes the value\n ``(A**2) * N/4`` for a harmonic signal with amplitude A for sufficiently\n large N.\n\n Parameters\n ----------\n x : array_like\n Sample times.\n y : array_like\n Measurement values.\n freqs : array_like\n Angular frequencies for output periodogram.\n\n Returns\n -------\n pgram : array_like\n Lomb-Scargle periodogram.\n\n Raises\n ------\n ValueError\n If the input arrays `x` and `y` do not have the same shape.\n\n Notes\n -----\n This subroutine calculates the periodogram using a slightly\n modified algorithm due to Townsend [3]_ which allows the\n periodogram to be calculated using only a single pass through\n the input arrays for each frequency.\n\n The algorithm running time scales roughly as O(x * freqs) or O(N^2)\n for a large number of samples and frequencies.\n\n References\n ----------\n .. [1] N.R. Lomb \"Least-squares frequency analysis of unequally spaced\n data\", Astrophysics and Space Science, vol 39, pp. 447-462, 1976\n\n .. [2] J.D. Scargle \"Studies in astronomical time series analysis. II - \n Statistical aspects of spectral analysis of unevenly spaced data\",\n The Astrophysical Journal, vol 263, pp. 835-853, 1982\n\n .. [3] R.H.D. Townsend, \"Fast calculation of the Lomb-Scargle\n periodogram using graphics processing units.\", The Astrophysical\n Journal Supplement Series, vol 191, pp. 247-253, 2010\n\n Examples\n --------\n >>> import scipy.signal\n\n First define some input parameters for the signal:\n\n >>> A = 2.\n >>> w = 1.\n >>> phi = 0.5 * np.""pi\n >>> nin = 1000\n >>> nout = 100000\n >>> frac_points = 0.9 # Fraction of points to select\n \n Randomly select a fraction of an array with timesteps:\n\n >>> r = np.random.rand(nin)\n >>> x = np.linspace(0.01, 10*np.pi, nin)\n >>> x = x[r >= frac_points]\n >>> normval = x.shape[0] # For normalization of the periodogram\n \n Plot a sine wave for the selected times:\n\n >>> y = A * np.sin(w*x+phi)\n\n Define the array of frequencies for which to compute the periodogram:\n \n >>> f = np.linspace(0.01, 10, nout)\n \n Calculate Lomb-Scargle periodogram:\n\n >>> pgram = sp.signal.lombscargle(x, y, f)\n\n Now make a plot of the input data:\n\n >>> plt.subplot(2, 1, 1)\n <matplotlib.axes.AxesSubplot object at 0x102154f50>\n >>> plt.plot(x, y, 'b+')\n [<matplotlib.lines.Line2D object at 0x102154a10>]\n\n Then plot the normalized periodogram:\n\n >>> plt.subplot(2, 1, 2)\n <matplotlib.axes.AxesSubplot object at 0x104b0a990>\n >>> plt.plot(f, np.sqrt(4*(pgram/normval)))\n [<matplotlib.lines.Line2D object at 0x104b2f910>]\n >>> plt.show()\n \n ";
static char __pyx_k__B[] = "B";
@@ -1125,10 +1128,10 @@ static PyObject *__pyx_k_tuple_16;
static PyObject *__pyx_k_codeobj_17;
/* Python wrapper */
-static PyObject *__pyx_pw_5scipy_6signal_8spectral_1lombscargle(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/
-static char __pyx_doc_5scipy_6signal_8spectral_lombscargle[] = "\n lombscargle(x, y, freqs)\n\n Computes the Lomb-Scargle periodogram.\n \n The Lomb-Scargle periodogram was developed by Lomb [1]_ and further\n extended by Scargle [2]_ to find, and test the significance of weak\n periodic signals with uneven temporal sampling.\n\n The computed periodogram is unnormalized, it takes the value\n ``(A**2) * N/4`` for a harmonic signal with amplitude A for sufficiently\n large N.\n\n Parameters\n ----------\n x : array_like\n Sample times.\n y : array_like\n Measurement values.\n freqs : array_like\n Angular frequencies for output periodogram.\n\n Returns\n -------\n pgram : array_like\n Lomb-Scargle periodogram.\n\n Raises\n ------\n ValueError\n If the input arrays `x` and `y` do not have the same shape.\n\n Notes\n -----\n This subroutine calculates the periodogram using a slightly\n modified algorithm due to Townsend [3]_ which allows the\n periodogram to be calculated using only a single pass through\n the input arrays for each frequency.\n\n The algorithm running time scales roughly as O(x * freqs) or O(N^2)\n for a large number of samples and frequencies.\n\n References\n ----------\n .. [1] N.R. Lomb \"Least-squares frequency analysis of unequally spaced\n data\", Astrophysics and Space Science, vol 39, pp. 447-462, 1976\n\n .. [2] J.D. Scargle \"Studies in astronomical time series analysis. II - \n Statistical aspects of spectral analysis of unevenly spaced data\",\n The Astrophysical Journal, vol 263, pp. 835-853, 1982\n\n .. [3] R.H.D. Townsend, \"Fast calculation of the Lomb-Scargle\n periodogram using graphics processing units.\", The Astrophysical\n Journal Supplement Series, vol 191, pp. 247-253, 2010\n\n Examples\n --------\n >>> import scipy.signal\n\n First define some input parameters for the signal:\n\n >>> A = 2.\n >>> w = 1.\n >>> phi = 0.5 * np.""pi\n >>> nin = 1000\n >>> nout = 100000\n >>> frac_points = 0.9 # Fraction of points to select\n \n Randomly select a fraction of an array with timesteps:\n\n >>> r = np.random.rand(nin)\n >>> x = np.linspace(0.01, 10*np.pi, nin)\n >>> x = x[r >= frac_points]\n >>> normval = x.shape[0] # For normalization of the periodogram\n \n Plot a sine wave for the selected times:\n\n >>> y = A * np.sin(w*x+phi)\n\n Define the array of frequencies for which to compute the periodogram:\n \n >>> f = np.linspace(0.01, 10, nout)\n \n Calculate Lomb-Scargle periodogram:\n\n >>> pgram = sp.signal.lombscargle(x, y, f)\n\n Now make a plot of the input data:\n\n >>> plt.subplot(2, 1, 1)\n <matplotlib.axes.AxesSubplot object at 0x102154f50>\n >>> plt.plot(x, y, 'b+')\n [<matplotlib.lines.Line2D object at 0x102154a10>]\n\n Then plot the normalized periodogram:\n\n >>> plt.subplot(2, 1, 2)\n <matplotlib.axes.AxesSubplot object at 0x104b0a990>\n >>> plt.plot(f, np.sqrt(4*(pgram/normval)))\n [<matplotlib.lines.Line2D object at 0x104b2f910>]\n >>> plt.show()\n \n ";
-static PyMethodDef __pyx_mdef_5scipy_6signal_8spectral_1lombscargle = {__Pyx_NAMESTR("lombscargle"), (PyCFunction)__pyx_pw_5scipy_6signal_8spectral_1lombscargle, METH_VARARGS|METH_KEYWORDS, __Pyx_DOCSTR(__pyx_doc_5scipy_6signal_8spectral_lombscargle)};
-static PyObject *__pyx_pw_5scipy_6signal_8spectral_1lombscargle(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {
+static PyObject *__pyx_pw_5scipy_6signal_9_spectral_1lombscargle(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/
+static char __pyx_doc_5scipy_6signal_9_spectral_lombscargle[] = "\n lombscargle(x, y, freqs)\n\n Computes the Lomb-Scargle periodogram.\n \n The Lomb-Scargle periodogram was developed by Lomb [1]_ and further\n extended by Scargle [2]_ to find, and test the significance of weak\n periodic signals with uneven temporal sampling.\n\n The computed periodogram is unnormalized, it takes the value\n ``(A**2) * N/4`` for a harmonic signal with amplitude A for sufficiently\n large N.\n\n Parameters\n ----------\n x : array_like\n Sample times.\n y : array_like\n Measurement values.\n freqs : array_like\n Angular frequencies for output periodogram.\n\n Returns\n -------\n pgram : array_like\n Lomb-Scargle periodogram.\n\n Raises\n ------\n ValueError\n If the input arrays `x` and `y` do not have the same shape.\n\n Notes\n -----\n This subroutine calculates the periodogram using a slightly\n modified algorithm due to Townsend [3]_ which allows the\n periodogram to be calculated using only a single pass through\n the input arrays for each frequency.\n\n The algorithm running time scales roughly as O(x * freqs) or O(N^2)\n for a large number of samples and frequencies.\n\n References\n ----------\n .. [1] N.R. Lomb \"Least-squares frequency analysis of unequally spaced\n data\", Astrophysics and Space Science, vol 39, pp. 447-462, 1976\n\n .. [2] J.D. Scargle \"Studies in astronomical time series analysis. II - \n Statistical aspects of spectral analysis of unevenly spaced data\",\n The Astrophysical Journal, vol 263, pp. 835-853, 1982\n\n .. [3] R.H.D. Townsend, \"Fast calculation of the Lomb-Scargle\n periodogram using graphics processing units.\", The Astrophysical\n Journal Supplement Series, vol 191, pp. 247-253, 2010\n\n Examples\n --------\n >>> import scipy.signal\n\n First define some input parameters for the signal:\n\n >>> A = 2.\n >>> w = 1.\n >>> phi = 0.5 * np.""pi\n >>> nin = 1000\n >>> nout = 100000\n >>> frac_points = 0.9 # Fraction of points to select\n \n Randomly select a fraction of an array with timesteps:\n\n >>> r = np.random.rand(nin)\n >>> x = np.linspace(0.01, 10*np.pi, nin)\n >>> x = x[r >= frac_points]\n >>> normval = x.shape[0] # For normalization of the periodogram\n \n Plot a sine wave for the selected times:\n\n >>> y = A * np.sin(w*x+phi)\n\n Define the array of frequencies for which to compute the periodogram:\n \n >>> f = np.linspace(0.01, 10, nout)\n \n Calculate Lomb-Scargle periodogram:\n\n >>> pgram = sp.signal.lombscargle(x, y, f)\n\n Now make a plot of the input data:\n\n >>> plt.subplot(2, 1, 1)\n <matplotlib.axes.AxesSubplot object at 0x102154f50>\n >>> plt.plot(x, y, 'b+')\n [<matplotlib.lines.Line2D object at 0x102154a10>]\n\n Then plot the normalized periodogram:\n\n >>> plt.subplot(2, 1, 2)\n <matplotlib.axes.AxesSubplot object at 0x104b0a990>\n >>> plt.plot(f, np.sqrt(4*(pgram/normval)))\n [<matplotlib.lines.Line2D object at 0x104b2f910>]\n >>> plt.show()\n \n ";
+static PyMethodDef __pyx_mdef_5scipy_6signal_9_spectral_1lombscargle = {__Pyx_NAMESTR("lombscargle"), (PyCFunction)__pyx_pw_5scipy_6signal_9_spectral_1lombscargle, METH_VARARGS|METH_KEYWORDS, __Pyx_DOCSTR(__pyx_doc_5scipy_6signal_9_spectral_lombscargle)};
+static PyObject *__pyx_pw_5scipy_6signal_9_spectral_1lombscargle(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {
PyArrayObject *__pyx_v_x = 0;
PyArrayObject *__pyx_v_y = 0;
PyArrayObject *__pyx_v_freqs = 0;
@@ -1182,14 +1185,14 @@ static PyObject *__pyx_pw_5scipy_6signal_8spectral_1lombscargle(PyObject *__pyx_
__pyx_L5_argtuple_error:;
__Pyx_RaiseArgtupleInvalid("lombscargle", 1, 3, 3, PyTuple_GET_SIZE(__pyx_args)); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 19; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
__pyx_L3_error:;
- __Pyx_AddTraceback("scipy.signal.spectral.lombscargle", __pyx_clineno, __pyx_lineno, __pyx_filename);
+ __Pyx_AddTraceback("scipy.signal._spectral.lombscargle", __pyx_clineno, __pyx_lineno, __pyx_filename);
__Pyx_RefNannyFinishContext();
return NULL;
__pyx_L4_argument_unpacking_done:;
if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_x), __pyx_ptype_5numpy_ndarray, 1, "x", 0))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 19; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_y), __pyx_ptype_5numpy_ndarray, 1, "y", 0))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 20; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_freqs), __pyx_ptype_5numpy_ndarray, 1, "freqs", 0))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 21; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __pyx_r = __pyx_pf_5scipy_6signal_8spectral_lombscargle(__pyx_self, __pyx_v_x, __pyx_v_y, __pyx_v_freqs);
+ __pyx_r = __pyx_pf_5scipy_6signal_9_spectral_lombscargle(__pyx_self, __pyx_v_x, __pyx_v_y, __pyx_v_freqs);
goto __pyx_L0;
__pyx_L1_error:;
__pyx_r = NULL;
@@ -1198,7 +1201,7 @@ static PyObject *__pyx_pw_5scipy_6signal_8spectral_1lombscargle(PyObject *__pyx_
return __pyx_r;
}
-/* "scipy/signal/spectral.pyx":19
+/* "scipy/signal/_spectral.pyx":19
*
* @cython.boundscheck(False)
* def lombscargle(np.ndarray[np.float64_t, ndim=1] x, # <<<<<<<<<<<<<<
@@ -1206,7 +1209,7 @@ static PyObject *__pyx_pw_5scipy_6signal_8spectral_1lombscargle(PyObject *__pyx_
* np.ndarray[np.float64_t, ndim=1] freqs):
*/
-static PyObject *__pyx_pf_5scipy_6signal_8spectral_lombscargle(CYTHON_UNUSED PyObject *__pyx_self, PyArrayObject *__pyx_v_x, PyArrayObject *__pyx_v_y, PyArrayObject *__pyx_v_freqs) {
+static PyObject *__pyx_pf_5scipy_6signal_9_spectral_lombscargle(CYTHON_UNUSED PyObject *__pyx_self, PyArrayObject *__pyx_v_x, PyArrayObject *__pyx_v_y, PyArrayObject *__pyx_v_freqs) {
PyObject *__pyx_v_pgram = NULL;
Py_ssize_t __pyx_v_i;
Py_ssize_t __pyx_v_j;
@@ -1286,7 +1289,7 @@ static PyObject *__pyx_pf_5scipy_6signal_8spectral_lombscargle(CYTHON_UNUSED PyO
}
__pyx_pybuffernd_freqs.diminfo[0].strides = __pyx_pybuffernd_freqs.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_freqs.diminfo[0].shape = __pyx_pybuffernd_freqs.rcbuffer->pybuffer.shape[0];
- /* "scipy/signal/spectral.pyx":127
+ /* "scipy/signal/_spectral.pyx":127
*
* # Check input sizes
* if x.shape[0] != y.shape[0]: # <<<<<<<<<<<<<<
@@ -1296,7 +1299,7 @@ static PyObject *__pyx_pf_5scipy_6signal_8spectral_lombscargle(CYTHON_UNUSED PyO
__pyx_t_1 = ((__pyx_v_x->dimensions[0]) != (__pyx_v_y->dimensions[0]));
if (__pyx_t_1) {
- /* "scipy/signal/spectral.pyx":128
+ /* "scipy/signal/_spectral.pyx":128
* # Check input sizes
* if x.shape[0] != y.shape[0]:
* raise ValueError("Input arrays do not have the same size.") # <<<<<<<<<<<<<<
@@ -1312,7 +1315,7 @@ static PyObject *__pyx_pf_5scipy_6signal_8spectral_lombscargle(CYTHON_UNUSED PyO
}
__pyx_L3:;
- /* "scipy/signal/spectral.pyx":131
+ /* "scipy/signal/_spectral.pyx":131
*
* # Create empty array for output periodogram
* pgram = np.empty(freqs.shape[0], dtype=np.float64) # <<<<<<<<<<<<<<
@@ -1348,7 +1351,7 @@ static PyObject *__pyx_pf_5scipy_6signal_8spectral_lombscargle(CYTHON_UNUSED PyO
__pyx_v_pgram = __pyx_t_6;
__pyx_t_6 = 0;
- /* "scipy/signal/spectral.pyx":138
+ /* "scipy/signal/_spectral.pyx":138
* cdef double tau, c_tau, s_tau, c_tau2, s_tau2, cs_tau
*
* for i in range(freqs.shape[0]): # <<<<<<<<<<<<<<
@@ -1359,7 +1362,7 @@ static PyObject *__pyx_pf_5scipy_6signal_8spectral_lombscargle(CYTHON_UNUSED PyO
for (__pyx_t_8 = 0; __pyx_t_8 < __pyx_t_7; __pyx_t_8+=1) {
__pyx_v_i = __pyx_t_8;
- /* "scipy/signal/spectral.pyx":140
+ /* "scipy/signal/_spectral.pyx":140
* for i in range(freqs.shape[0]):
*
* xc = 0. # <<<<<<<<<<<<<<
@@ -1368,7 +1371,7 @@ static PyObject *__pyx_pf_5scipy_6signal_8spectral_lombscargle(CYTHON_UNUSED PyO
*/
__pyx_v_xc = 0.;
- /* "scipy/signal/spectral.pyx":141
+ /* "scipy/signal/_spectral.pyx":141
*
* xc = 0.
* xs = 0. # <<<<<<<<<<<<<<
@@ -1377,7 +1380,7 @@ static PyObject *__pyx_pf_5scipy_6signal_8spectral_lombscargle(CYTHON_UNUSED PyO
*/
__pyx_v_xs = 0.;
- /* "scipy/signal/spectral.pyx":142
+ /* "scipy/signal/_spectral.pyx":142
* xc = 0.
* xs = 0.
* cc = 0. # <<<<<<<<<<<<<<
@@ -1386,7 +1389,7 @@ static PyObject *__pyx_pf_5scipy_6signal_8spectral_lombscargle(CYTHON_UNUSED PyO
*/
__pyx_v_cc = 0.;
- /* "scipy/signal/spectral.pyx":143
+ /* "scipy/signal/_spectral.pyx":143
* xs = 0.
* cc = 0.
* ss = 0. # <<<<<<<<<<<<<<
@@ -1395,7 +1398,7 @@ static PyObject *__pyx_pf_5scipy_6signal_8spectral_lombscargle(CYTHON_UNUSED PyO
*/
__pyx_v_ss = 0.;
- /* "scipy/signal/spectral.pyx":144
+ /* "scipy/signal/_spectral.pyx":144
* cc = 0.
* ss = 0.
* cs = 0. # <<<<<<<<<<<<<<
@@ -1404,7 +1407,7 @@ static PyObject *__pyx_pf_5scipy_6signal_8spectral_lombscargle(CYTHON_UNUSED PyO
*/
__pyx_v_cs = 0.;
- /* "scipy/signal/spectral.pyx":146
+ /* "scipy/signal/_spectral.pyx":146
* cs = 0.
*
* for j in range(x.shape[0]): # <<<<<<<<<<<<<<
@@ -1415,7 +1418,7 @@ static PyObject *__pyx_pf_5scipy_6signal_8spectral_lombscargle(CYTHON_UNUSED PyO
for (__pyx_t_10 = 0; __pyx_t_10 < __pyx_t_9; __pyx_t_10+=1) {
__pyx_v_j = __pyx_t_10;
- /* "scipy/signal/spectral.pyx":148
+ /* "scipy/signal/_spectral.pyx":148
* for j in range(x.shape[0]):
*
* c = cos(freqs[i] * x[j]) # <<<<<<<<<<<<<<
@@ -1428,7 +1431,7 @@ static PyObject *__pyx_pf_5scipy_6signal_8spectral_lombscargle(CYTHON_UNUSED PyO
if (__pyx_t_12 < 0) __pyx_t_12 += __pyx_pybuffernd_x.diminfo[0].shape;
__pyx_v_c = cos(((*__Pyx_BufPtrStrided1d(__pyx_t_5numpy_float64_t *, __pyx_pybuffernd_freqs.rcbuffer->pybuffer.buf, __pyx_t_11, __pyx_pybuffernd_freqs.diminfo[0].strides)) * (*__Pyx_BufPtrStrided1d(__pyx_t_5numpy_float64_t *, __pyx_pybuffernd_x.rcbuffer->pybuffer.buf, __pyx_t_12, __pyx_pybuffernd_x.diminfo[0].strides))));
- /* "scipy/signal/spectral.pyx":149
+ /* "scipy/signal/_spectral.pyx":149
*
* c = cos(freqs[i] * x[j])
* s = sin(freqs[i] * x[j]) # <<<<<<<<<<<<<<
@@ -1441,7 +1444,7 @@ static PyObject *__pyx_pf_5scipy_6signal_8spectral_lombscargle(CYTHON_UNUSED PyO
if (__pyx_t_14 < 0) __pyx_t_14 += __pyx_pybuffernd_x.diminfo[0].shape;
__pyx_v_s = sin(((*__Pyx_BufPtrStrided1d(__pyx_t_5numpy_float64_t *, __pyx_pybuffernd_freqs.rcbuffer->pybuffer.buf, __pyx_t_13, __pyx_pybuffernd_freqs.diminfo[0].strides)) * (*__Pyx_BufPtrStrided1d(__pyx_t_5numpy_float64_t *, __pyx_pybuffernd_x.rcbuffer->pybuffer.buf, __pyx_t_14, __pyx_pybuffernd_x.diminfo[0].strides))));
- /* "scipy/signal/spectral.pyx":151
+ /* "scipy/signal/_spectral.pyx":151
* s = sin(freqs[i] * x[j])
*
* xc += y[j] * c # <<<<<<<<<<<<<<
@@ -1452,7 +1455,7 @@ static PyObject *__pyx_pf_5scipy_6signal_8spectral_lombscargle(CYTHON_UNUSED PyO
if (__pyx_t_15 < 0) __pyx_t_15 += __pyx_pybuffernd_y.diminfo[0].shape;
__pyx_v_xc = (__pyx_v_xc + ((*__Pyx_BufPtrStrided1d(__pyx_t_5numpy_float64_t *, __pyx_pybuffernd_y.rcbuffer->pybuffer.buf, __pyx_t_15, __pyx_pybuffernd_y.diminfo[0].strides)) * __pyx_v_c));
- /* "scipy/signal/spectral.pyx":152
+ /* "scipy/signal/_spectral.pyx":152
*
* xc += y[j] * c
* xs += y[j] * s # <<<<<<<<<<<<<<
@@ -1463,7 +1466,7 @@ static PyObject *__pyx_pf_5scipy_6signal_8spectral_lombscargle(CYTHON_UNUSED PyO
if (__pyx_t_16 < 0) __pyx_t_16 += __pyx_pybuffernd_y.diminfo[0].shape;
__pyx_v_xs = (__pyx_v_xs + ((*__Pyx_BufPtrStrided1d(__pyx_t_5numpy_float64_t *, __pyx_pybuffernd_y.rcbuffer->pybuffer.buf, __pyx_t_16, __pyx_pybuffernd_y.diminfo[0].strides)) * __pyx_v_s));
- /* "scipy/signal/spectral.pyx":153
+ /* "scipy/signal/_spectral.pyx":153
* xc += y[j] * c
* xs += y[j] * s
* cc += c * c # <<<<<<<<<<<<<<
@@ -1472,7 +1475,7 @@ static PyObject *__pyx_pf_5scipy_6signal_8spectral_lombscargle(CYTHON_UNUSED PyO
*/
__pyx_v_cc = (__pyx_v_cc + (__pyx_v_c * __pyx_v_c));
- /* "scipy/signal/spectral.pyx":154
+ /* "scipy/signal/_spectral.pyx":154
* xs += y[j] * s
* cc += c * c
* ss += s * s # <<<<<<<<<<<<<<
@@ -1481,7 +1484,7 @@ static PyObject *__pyx_pf_5scipy_6signal_8spectral_lombscargle(CYTHON_UNUSED PyO
*/
__pyx_v_ss = (__pyx_v_ss + (__pyx_v_s * __pyx_v_s));
- /* "scipy/signal/spectral.pyx":155
+ /* "scipy/signal/_spectral.pyx":155
* cc += c * c
* ss += s * s
* cs += c * s # <<<<<<<<<<<<<<
@@ -1491,7 +1494,7 @@ static PyObject *__pyx_pf_5scipy_6signal_8spectral_lombscargle(CYTHON_UNUSED PyO
__pyx_v_cs = (__pyx_v_cs + (__pyx_v_c * __pyx_v_s));
}
- /* "scipy/signal/spectral.pyx":157
+ /* "scipy/signal/_spectral.pyx":157
* cs += c * s
*
* tau = atan(2 * cs / (cc - ss)) / (2 * freqs[i]) # <<<<<<<<<<<<<<
@@ -1514,7 +1517,7 @@ static PyObject *__pyx_pf_5scipy_6signal_8spectral_lombscargle(CYTHON_UNUSED PyO
}
__pyx_v_tau = (__pyx_t_19 / __pyx_t_20);
- /* "scipy/signal/spectral.pyx":158
+ /* "scipy/signal/_spectral.pyx":158
*
* tau = atan(2 * cs / (cc - ss)) / (2 * freqs[i])
* c_tau = cos(freqs[i] * tau) # <<<<<<<<<<<<<<
@@ -1525,7 +1528,7 @@ static PyObject *__pyx_pf_5scipy_6signal_8spectral_lombscargle(CYTHON_UNUSED PyO
if (__pyx_t_21 < 0) __pyx_t_21 += __pyx_pybuffernd_freqs.diminfo[0].shape;
__pyx_v_c_tau = cos(((*__Pyx_BufPtrStrided1d(__pyx_t_5numpy_float64_t *, __pyx_pybuffernd_freqs.rcbuffer->pybuffer.buf, __pyx_t_21, __pyx_pybuffernd_freqs.diminfo[0].strides)) * __pyx_v_tau));
- /* "scipy/signal/spectral.pyx":159
+ /* "scipy/signal/_spectral.pyx":159
* tau = atan(2 * cs / (cc - ss)) / (2 * freqs[i])
* c_tau = cos(freqs[i] * tau)
* s_tau = sin(freqs[i] * tau) # <<<<<<<<<<<<<<
@@ -1536,7 +1539,7 @@ static PyObject *__pyx_pf_5scipy_6signal_8spectral_lombscargle(CYTHON_UNUSED PyO
if (__pyx_t_22 < 0) __pyx_t_22 += __pyx_pybuffernd_freqs.diminfo[0].shape;
__pyx_v_s_tau = sin(((*__Pyx_BufPtrStrided1d(__pyx_t_5numpy_float64_t *, __pyx_pybuffernd_freqs.rcbuffer->pybuffer.buf, __pyx_t_22, __pyx_pybuffernd_freqs.diminfo[0].strides)) * __pyx_v_tau));
- /* "scipy/signal/spectral.pyx":160
+ /* "scipy/signal/_spectral.pyx":160
* c_tau = cos(freqs[i] * tau)
* s_tau = sin(freqs[i] * tau)
* c_tau2 = c_tau * c_tau # <<<<<<<<<<<<<<
@@ -1545,7 +1548,7 @@ static PyObject *__pyx_pf_5scipy_6signal_8spectral_lombscargle(CYTHON_UNUSED PyO
*/
__pyx_v_c_tau2 = (__pyx_v_c_tau * __pyx_v_c_tau);
- /* "scipy/signal/spectral.pyx":161
+ /* "scipy/signal/_spectral.pyx":161
* s_tau = sin(freqs[i] * tau)
* c_tau2 = c_tau * c_tau
* s_tau2 = s_tau * s_tau # <<<<<<<<<<<<<<
@@ -1554,7 +1557,7 @@ static PyObject *__pyx_pf_5scipy_6signal_8spectral_lombscargle(CYTHON_UNUSED PyO
*/
__pyx_v_s_tau2 = (__pyx_v_s_tau * __pyx_v_s_tau);
- /* "scipy/signal/spectral.pyx":162
+ /* "scipy/signal/_spectral.pyx":162
* c_tau2 = c_tau * c_tau
* s_tau2 = s_tau * s_tau
* cs_tau = 2 * c_tau * s_tau # <<<<<<<<<<<<<<
@@ -1563,7 +1566,7 @@ static PyObject *__pyx_pf_5scipy_6signal_8spectral_lombscargle(CYTHON_UNUSED PyO
*/
__pyx_v_cs_tau = ((2.0 * __pyx_v_c_tau) * __pyx_v_s_tau);
- /* "scipy/signal/spectral.pyx":164
+ /* "scipy/signal/_spectral.pyx":164
* cs_tau = 2 * c_tau * s_tau
*
* pgram[i] = 0.5 * (((c_tau * xc + s_tau * xs)**2 / \ # <<<<<<<<<<<<<<
@@ -1572,7 +1575,7 @@ static PyObject *__pyx_pf_5scipy_6signal_8spectral_lombscargle(CYTHON_UNUSED PyO
*/
__pyx_t_19 = pow(((__pyx_v_c_tau * __pyx_v_xc) + (__pyx_v_s_tau * __pyx_v_xs)), 2.0);
- /* "scipy/signal/spectral.pyx":165
+ /* "scipy/signal/_spectral.pyx":165
*
* pgram[i] = 0.5 * (((c_tau * xc + s_tau * xs)**2 / \
* (c_tau2 * cc + cs_tau * cs + s_tau2 * ss)) + \ # <<<<<<<<<<<<<<
@@ -1585,7 +1588,7 @@ static PyObject *__pyx_pf_5scipy_6signal_8spectral_lombscargle(CYTHON_UNUSED PyO
{__pyx_filename = __pyx_f[0]; __pyx_lineno = 164; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
}
- /* "scipy/signal/spectral.pyx":166
+ /* "scipy/signal/_spectral.pyx":166
* pgram[i] = 0.5 * (((c_tau * xc + s_tau * xs)**2 / \
* (c_tau2 * cc + cs_tau * cs + s_tau2 * ss)) + \
* ((c_tau * xs - s_tau * xc)**2 / \ # <<<<<<<<<<<<<<
@@ -1594,7 +1597,7 @@ static PyObject *__pyx_pf_5scipy_6signal_8spectral_lombscargle(CYTHON_UNUSED PyO
*/
__pyx_t_17 = pow(((__pyx_v_c_tau * __pyx_v_xs) - (__pyx_v_s_tau * __pyx_v_xc)), 2.0);
- /* "scipy/signal/spectral.pyx":167
+ /* "scipy/signal/_spectral.pyx":167
* (c_tau2 * cc + cs_tau * cs + s_tau2 * ss)) + \
* ((c_tau * xs - s_tau * xc)**2 / \
* (c_tau2 * ss - cs_tau * cs + s_tau2 * cc))) # <<<<<<<<<<<<<<
@@ -1609,7 +1612,7 @@ static PyObject *__pyx_pf_5scipy_6signal_8spectral_lombscargle(CYTHON_UNUSED PyO
__pyx_t_6 = PyFloat_FromDouble((0.5 * ((__pyx_t_19 / __pyx_t_18) + (__pyx_t_17 / __pyx_t_23)))); if (unlikely(!__pyx_t_6)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 164; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_6);
- /* "scipy/signal/spectral.pyx":164
+ /* "scipy/signal/_spectral.pyx":164
* cs_tau = 2 * c_tau * s_tau
*
* pgram[i] = 0.5 * (((c_tau * xc + s_tau * xs)**2 / \ # <<<<<<<<<<<<<<
@@ -1620,7 +1623,7 @@ static PyObject *__pyx_pf_5scipy_6signal_8spectral_lombscargle(CYTHON_UNUSED PyO
__Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
}
- /* "scipy/signal/spectral.pyx":169
+ /* "scipy/signal/_spectral.pyx":169
* (c_tau2 * ss - cs_tau * cs + s_tau2 * cc)))
*
* return pgram # <<<<<<<<<<<<<<
@@ -1645,7 +1648,7 @@ static PyObject *__pyx_pf_5scipy_6signal_8spectral_lombscargle(CYTHON_UNUSED PyO
__Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_x.rcbuffer->pybuffer);
__Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_y.rcbuffer->pybuffer);
__Pyx_ErrRestore(__pyx_type, __pyx_value, __pyx_tb);}
- __Pyx_AddTraceback("scipy.signal.spectral.lombscargle", __pyx_clineno, __pyx_lineno, __pyx_filename);
+ __Pyx_AddTraceback("scipy.signal._spectral.lombscargle", __pyx_clineno, __pyx_lineno, __pyx_filename);
__pyx_r = NULL;
goto __pyx_L2;
__pyx_L0:;
@@ -2005,8 +2008,10 @@ static int __pyx_pf_5numpy_7ndarray___getbuffer__(PyArrayObject *__pyx_v_self, P
* cdef list stack
* cdef int offset
*/
- __Pyx_INCREF(((PyObject *)__pyx_v_self->descr));
- __pyx_v_descr = __pyx_v_self->descr;
+ __pyx_t_4 = ((PyObject *)__pyx_v_self->descr);
+ __Pyx_INCREF(__pyx_t_4);
+ __pyx_v_descr = ((PyArray_Descr *)__pyx_t_4);
+ __pyx_t_4 = 0;
/* "numpy.pxd":244
* cdef int offset
@@ -2081,7 +2086,8 @@ static int __pyx_pf_5numpy_7ndarray___getbuffer__(PyArrayObject *__pyx_v_self, P
* if ((descr.byteorder == c'>' and little_endian) or
* (descr.byteorder == c'<' and not little_endian)):
*/
- __pyx_v_t = __pyx_v_descr->type_num;
+ __pyx_t_5 = __pyx_v_descr->type_num;
+ __pyx_v_t = __pyx_t_5;
/* "numpy.pxd":255
* if not hasfields:
@@ -3636,7 +3642,7 @@ static PyMethodDef __pyx_methods[] = {
#if PY_MAJOR_VERSION >= 3
static struct PyModuleDef __pyx_moduledef = {
PyModuleDef_HEAD_INIT,
- __Pyx_NAMESTR("spectral"),
+ __Pyx_NAMESTR("_spectral"),
__Pyx_DOCSTR(__pyx_k_15), /* m_doc */
-1, /* m_size */
__pyx_methods /* m_methods */,
@@ -3705,7 +3711,7 @@ static int __Pyx_InitCachedConstants(void) {
__Pyx_RefNannyDeclarations
__Pyx_RefNannySetupContext("__Pyx_InitCachedConstants", 0);
- /* "scipy/signal/spectral.pyx":128
+ /* "scipy/signal/_spectral.pyx":128
* # Check input sizes
* if x.shape[0] != y.shape[0]:
* raise ValueError("Input arrays do not have the same size.") # <<<<<<<<<<<<<<
@@ -3803,7 +3809,7 @@ static int __Pyx_InitCachedConstants(void) {
__Pyx_GIVEREF(((PyObject *)__pyx_kp_u_13));
__Pyx_GIVEREF(((PyObject *)__pyx_k_tuple_14));
- /* "scipy/signal/spectral.pyx":19
+ /* "scipy/signal/_spectral.pyx":19
*
* @cython.boundscheck(False)
* def lombscargle(np.ndarray[np.float64_t, ndim=1] x, # <<<<<<<<<<<<<<
@@ -3887,11 +3893,11 @@ static int __Pyx_InitGlobals(void) {
}
#if PY_MAJOR_VERSION < 3
-PyMODINIT_FUNC initspectral(void); /*proto*/
-PyMODINIT_FUNC initspectral(void)
+PyMODINIT_FUNC init_spectral(void); /*proto*/
+PyMODINIT_FUNC init_spectral(void)
#else
-PyMODINIT_FUNC PyInit_spectral(void); /*proto*/
-PyMODINIT_FUNC PyInit_spectral(void)
+PyMODINIT_FUNC PyInit__spectral(void); /*proto*/
+PyMODINIT_FUNC PyInit__spectral(void)
#endif
{
PyObject *__pyx_t_1 = NULL;
@@ -3905,7 +3911,7 @@ PyMODINIT_FUNC PyInit_spectral(void)
Py_FatalError("failed to import 'refnanny' module");
}
#endif
- __Pyx_RefNannySetupContext("PyMODINIT_FUNC PyInit_spectral(void)", 0);
+ __Pyx_RefNannySetupContext("PyMODINIT_FUNC PyInit__spectral(void)", 0);
if ( __Pyx_check_binary_version() < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__pyx_empty_tuple = PyTuple_New(0); if (unlikely(!__pyx_empty_tuple)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__pyx_empty_bytes = PyBytes_FromStringAndSize("", 0); if (unlikely(!__pyx_empty_bytes)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
@@ -3927,11 +3933,19 @@ PyMODINIT_FUNC PyInit_spectral(void)
#endif
/*--- Module creation code ---*/
#if PY_MAJOR_VERSION < 3
- __pyx_m = Py_InitModule4(__Pyx_NAMESTR("spectral"), __pyx_methods, __Pyx_DOCSTR(__pyx_k_15), 0, PYTHON_API_VERSION); Py_XINCREF(__pyx_m);
+ __pyx_m = Py_InitModule4(__Pyx_NAMESTR("_spectral"), __pyx_methods, __Pyx_DOCSTR(__pyx_k_15), 0, PYTHON_API_VERSION); Py_XINCREF(__pyx_m);
#else
__pyx_m = PyModule_Create(&__pyx_moduledef);
#endif
if (unlikely(!__pyx_m)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ #if PY_MAJOR_VERSION >= 3
+ {
+ PyObject *modules = PyImport_GetModuleDict(); if (unlikely(!modules)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ if (!PyDict_GetItemString(modules, "scipy.signal._spectral")) {
+ if (unlikely(PyDict_SetItemString(modules, "scipy.signal._spectral", __pyx_m) < 0)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ }
+ }
+ #endif
__pyx_b = PyImport_AddModule(__Pyx_NAMESTR(__Pyx_BUILTIN_MODULE_NAME)); if (unlikely(!__pyx_b)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
#if CYTHON_COMPILING_IN_PYPY
Py_INCREF(__pyx_b);
@@ -3939,7 +3953,7 @@ PyMODINIT_FUNC PyInit_spectral(void)
if (__Pyx_SetAttrString(__pyx_m, "__builtins__", __pyx_b) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;};
/*--- Initialize various global constants etc. ---*/
if (unlikely(__Pyx_InitGlobals() < 0)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- if (__pyx_module_is_main_scipy__signal__spectral) {
+ if (__pyx_module_is_main_scipy__signal___spectral) {
if (__Pyx_SetAttrString(__pyx_m, "__name__", __pyx_n_s____main__) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;};
}
/*--- Builtin init code ---*/
@@ -3967,7 +3981,7 @@ PyMODINIT_FUNC PyInit_spectral(void)
/*--- Function import code ---*/
/*--- Execution code ---*/
- /* "scipy/signal/spectral.pyx":6
+ /* "scipy/signal/_spectral.pyx":6
* """Tools for spectral analysis of unequally sampled signals."""
*
* import numpy as np # <<<<<<<<<<<<<<
@@ -3979,7 +3993,7 @@ PyMODINIT_FUNC PyInit_spectral(void)
if (PyObject_SetAttr(__pyx_m, __pyx_n_s__np, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 6; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
- /* "scipy/signal/spectral.pyx":10
+ /* "scipy/signal/_spectral.pyx":10
* cimport cython
*
* __all__ = ['lombscargle'] # <<<<<<<<<<<<<<
@@ -3994,19 +4008,19 @@ PyMODINIT_FUNC PyInit_spectral(void)
if (PyObject_SetAttr(__pyx_m, __pyx_n_s____all__, ((PyObject *)__pyx_t_1)) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 10; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_DECREF(((PyObject *)__pyx_t_1)); __pyx_t_1 = 0;
- /* "scipy/signal/spectral.pyx":19
+ /* "scipy/signal/_spectral.pyx":19
*
* @cython.boundscheck(False)
* def lombscargle(np.ndarray[np.float64_t, ndim=1] x, # <<<<<<<<<<<<<<
* np.ndarray[np.float64_t, ndim=1] y,
* np.ndarray[np.float64_t, ndim=1] freqs):
*/
- __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_5scipy_6signal_8spectral_1lombscargle, NULL, __pyx_n_s_19); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 19; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_5scipy_6signal_9_spectral_1lombscargle, NULL, __pyx_n_s_19); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 19; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_1);
if (PyObject_SetAttr(__pyx_m, __pyx_n_s__lombscargle, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 19; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
- /* "scipy/signal/spectral.pyx":1
+ /* "scipy/signal/_spectral.pyx":1
* # Author: Pim Schellart # <<<<<<<<<<<<<<
* # 2010 - 2011
*
@@ -4028,10 +4042,10 @@ PyMODINIT_FUNC PyInit_spectral(void)
__pyx_L1_error:;
__Pyx_XDECREF(__pyx_t_1);
if (__pyx_m) {
- __Pyx_AddTraceback("init scipy.signal.spectral", __pyx_clineno, __pyx_lineno, __pyx_filename);
+ __Pyx_AddTraceback("init scipy.signal._spectral", __pyx_clineno, __pyx_lineno, __pyx_filename);
Py_DECREF(__pyx_m); __pyx_m = 0;
} else if (!PyErr_Occurred()) {
- PyErr_SetString(PyExc_ImportError, "init scipy.signal.spectral");
+ PyErr_SetString(PyExc_ImportError, "init scipy.signal._spectral");
}
__pyx_L0:;
__Pyx_RefNannyFinishContext();
0  scipy/signal/spectral.pyx → scipy/signal/_spectral.pyx
View
File renamed without changes
4 scipy/signal/bento.info
View
@@ -8,8 +8,8 @@ Library:
sigtoolsmodule.c,
firfilter.c,
medianfilter.c
- Extension: spectral
- Sources: spectral.c
+ Extension: _spectral
+ Sources: _spectral.c
Extension: spline
Sources:
splinemodule.c,
2  scipy/signal/setup.py
View
@@ -16,7 +16,7 @@ def configuration(parent_package='', top_path=None):
include_dirs=['.']
)
- config.add_extension('spectral', sources=['spectral.c'])
+ config.add_extension('_spectral', sources=['_spectral.c'])
config.add_extension('spline',
sources=['splinemodule.c', 'S_bspline_util.c', 'D_bspline_util.c',
368 scipy/signal/spectral.py
View
@@ -0,0 +1,368 @@
+"""Tools for spectral analysis.
+"""
+
+import numpy as np
+from scipy import fftpack
+import signaltools
+from windows import get_window
+from _spectral import lombscargle
+import warnings
+
+__all__ = ['periodogram', 'welch', 'lombscargle']
+
+
+def periodogram(x, fs=1.0, window=None, nfft=None, detrend='constant',
+ return_onesided=True, scaling='density', axis=-1):
+ """
+ Estimate power spectral density using a periodogram.
+
+ Parameters
+ ----------
+ x : array_like
+ Time series of measurement values
+ fs : float, optional
+ Sampling frequency of the `x` time series in units of Hz. Defaults
+ to 1.0.
+ window : str or tuple or array_like, optional
+ Desired window to use. See `get_window` for a list of windows and
+ required parameters. If `window` is an array it will be used
+ directly as the window. Defaults to None; equivalent to 'boxcar'.
+ nfft : int, optional
+ Length of the FFT used. If None the length of `x` will be used.
+ detrend : str or function, optional
+ Specifies how to detrend `x` prior to computing the spectrum. If
+ `detrend` is a string, it is passed as the ``type`` argument to
+ `detrend`. If it is a function, it should return a detrended array.
+ Defaults to 'constant'.
+ return_onesided : bool, optional
+ If True, return a one-sided spectrum for real data. If False return
+ a two-sided spectrum. Note that for complex data, a two-sided
+ spectrum is always returned.
+ scaling : { 'density', 'spectrum' }, optional
+ Selects between computing the power spectral density ('density')
+ where `Pxx` has units of V**2/Hz if `x` is measured in V and computing
+ the power spectrum ('spectrum') where `Pxx` has units of V**2 if `x` is
+ measured in V. Defaults to 'density'
+ axis : int, optional
+ Axis along which the periodogram is computed; the default is over
+ the last axis (i.e. ``axis=-1``).
+
+ Returns
+ -------
+ f : ndarray
+ Array of sample frequencies.
+ Pxx : ndarray
+ Power spectral density or power spectrum of `x`.
+
+ See Also
+ --------
+ welch: Estimate power spectral density using Welch's method
+ lombscargle: Lomb-Scargle periodogram for unevenly sampled data
+
+ Examples
+ --------
+ >>> from scipy import signal
Ralf Gommers Owner
rgommers added a note

Needs also import matplotlib.pyplot as plt.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
+ >>> import matplotlib.pyplot as plt
+
+ Generate a test signal, a 2 Vrms sine wave at 1234 Hz, corrupted by
+ 0.001 V**2/Hz of white noise sampled at 10 kHz.
+
+ >>> fs = 10e3
+ >>> N = 1e5
+ >>> amp = 2*np.sqrt(2)
+ >>> freq = 1234.0
+ >>> noise_power = 0.001 * fs / 2
+ >>> time = np.arange(N) / fs
+ >>> x = amp*np.sin(2*np.pi*freq*time)
+ >>> x += np.random.normal(scale=np.sqrt(noise_power), size=time.shape)
+
+ Compute and plot the power spectral density.
+
+ >>> f, Pxx_den = signal.periodogram(x, fs)
+ >>> plt.semilogy(f, Pxx_den)
+ >>> plt.xlabel('frequency [Hz]')
+ >>> plt.ylabel('PSD [V/sqrt(Hz)]')
Ralf Gommers Owner
rgommers added a note

Needs plt.show() here and in second example below.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
+ >>> plt.show()
+
+ If we average the last half of the spectral density, to exclude the
+ peak, we can recover the noise power on the signal.
+
+ >>> np.mean(Pxx_den[256:])
+ 0.0009924865443739191
+
+ Now compute and plot the power spectrum.
+
+ >>> f, Pxx_spec = signal.periodogram(x, fs, 'flattop', scaling='spectrum')
+ >>> plt.figure()
+ >>> plt.semilogy(f, np.sqrt(Pxx_spec))
+ >>> plt.xlabel('frequency [Hz]')
+ >>> plt.ylabel('Linear spectrum [V RMS]')
+ >>> plt.show()
+
+ The peak height in the power spectrum is an estimate of the RMS amplitude.
+
+ >>> np.sqrt(Pxx_spec.max())
+ 2.0077340678640727
+
+ """
+ x = np.asarray(x)
+
+ if x.size == 0:
+ return np.empty(x.shape), np.empty(x.shape)
+
+ if window is None:
+ window = 'boxcar'
+
+ if nfft is None:
+ nperseg = x.shape[axis]
+ elif nfft == x.shape[axis]:
+ nperseg = nfft
+ elif nfft > x.shape[axis]:
+ nperseg = x.shape[axis]
+ elif nfft < x.shape[axis]:
+ s = [np.s_[:]]*len(x.shape)
+ s[axis] = np.s_[:nfft]
+ x = x[s]
+ nperseg = nfft
+ nfft = None
+
+ return welch(x, fs, window, nperseg, 0, nfft, detrend, return_onesided,
+ scaling, axis)
+
+
+def welch(x, fs=1.0, window='hanning', nperseg=256, noverlap=None, nfft=None,
+ detrend='constant', return_onesided=True, scaling='density', axis=-1):
+ """
+ Estimate power spectral density using Welch's method.
+
+ Welch's method [1]_ computes an estimate of the power spectral density
+ by dividing the data into overlapping segments, computing a modified
+ periodogram for each segment and averaging the periodograms.
+
+ Parameters
+ ----------
+ x : array_like
+ Time series of measurement values
+ fs : float, optional
+ Sampling frequency of the `x` time series in units of Hz. Defaults
+ to 1.0.
+ window : str or tuple or array_like, optional
+ Desired window to use. See `get_window` for a list of windows and
+ required parameters. If `window` is array_like it will be used
+ directly as the window and its length will be used for nperseg.
+ Defaults to 'hanning'.
+ nperseg : int, optional
+ Length of each segment. Defaults to 256.
+ noverlap: int, optional
+ Number of points to overlap between segments. If None,
+ ``noverlap = nperseg / 2``. Defaults to None.
+ nfft : int, optional
+ Length of the FFT used, if a zero padded FFT is desired. If None,
+ the FFT length is `nperseg`. Defaults to None.
+ detrend : str or function, optional
+ Specifies how to detrend each segment. If `detrend` is a string,
+ it is passed as the ``type`` argument to `detrend`. If it is a
+ function, it takes a segment and returns a detrended segment.
+ Defaults to 'constant'.
+ return_onesided : bool, optional
+ If True, return a one-sided spectrum for real data. If False return
+ a two-sided spectrum. Note that for complex data, a two-sided
+ spectrum is always returned.
+ scaling : { 'density', 'spectrum' }, optional
+ Selects between computing the power spectral density ('density')
+ where Pxx has units of V**2/Hz if x is measured in V and computing
+ the power spectrum ('spectrum') where Pxx has units of V**2 if x is
+ measured in V. Defaults to 'density'.
+ axis : int, optional
+ Axis along which the periodogram is computed; the default is over
+ the last axis (i.e. ``axis=-1``).
+
+ Returns
+ -------
+ f : ndarray
+ Array of sample frequencies.
+ Pxx : ndarray
+ Power spectral density or power spectrum of x.
+
+ See Also
+ --------
+ periodogram: Simple, optionally modified periodogram
+ lombscargle: Lomb-Scargle periodogram for unevenly sampled data
+
+ Notes
+ -----
+ An appropriate amount of overlap will depend on the choice of window
+ and on your requirements. For the default 'hanning' window an
+ overlap of 50% is a reasonable trade off between accurately estimating
+ the signal power, while not over counting any of the data. Narrower
+ windows may require a larger overlap.
+
+ If `noverlap` is 0, this method is equivalent to Bartlett's method [2]_.
+
+ References
+ ----------
+ .. [1] P. Welch, "The use of the fast Fourier transform for the
+ estimation of power spectra: A method based on time averaging
+ over short, modified periodograms", IEEE Trans. Audio
+ Electroacoust. vol. 15, pp. 70-73, 1967.
+ .. [2] M.S. Bartlett, "Periodogram Analysis and Continuous Spectra",
+ Biometrika, vol. 37, pp. 1-16, 1950.
+
+ Examples
+ --------
+ >>> from scipy import signal
+ >>> import matplotlib.pyplot as plt
+
+ Generate a test signal, a 2 Vrms sine wave at 1234 Hz, corrupted by
+ 0.001 V**2/Hz of white noise sampled at 10 kHz.
+
+ >>> fs = 10e3
+ >>> N = 1e5
+ >>> amp = 2*np.sqrt(2)
+ >>> freq = 1234.0
+ >>> noise_power = 0.001 * fs / 2
+ >>> time = np.arange(N) / fs
+ >>> x = amp*np.sin(2*np.pi*freq*time)
+ >>> x += np.random.normal(scale=np.sqrt(noise_power), size=time.shape)
+
+ Compute and plot the power spectral density.
+
+ >>> f, Pxx_den = signal.welch(x, fs, 1024)
+ >>> plt.semilogy(f, Pxx_den)
+ >>> plt.xlabel('frequency [Hz]')
+ >>> plt.ylabel('PSD [V/sqrt(Hz)]')
+ >>> plt.show()
+
+ If we average the last half of the spectral density, to exclude the
+ peak, we can recover the noise power on the signal.
+
+ >>> np.mean(Pxx_den[256:])
+ 0.0009924865443739191
+
+ Now compute and plot the power spectrum.
+
+ >>> f, Pxx_spec = signal.welch(x, fs, 'flattop', 1024, scaling='spectrum')
+ >>> plt.figure()
+ >>> plt.semilogy(f, np.sqrt(Pxx_spec))
+ >>> plt.xlabel('frequency [Hz]')
+ >>> plt.ylabel('Linear spectrum [V RMS]')
+ >>> plt.show()
+
+ The peak height in the power spectrum is an estimate of the RMS amplitude.
+
+ >>> np.sqrt(Pxx_spec.max())
+ 2.0077340678640727
+
+ """
+ x = np.asarray(x)
+
+ if x.size == 0:
+ return np.empty(x.shape), np.empty(x.shape)
+
+ if axis != -1:
+ x = np.rollaxis(x, axis, len(x.shape))
+
+ if x.shape[-1] < nperseg:
+ warnings.warn('nperseg = %d, is greater than x.shape[%d] = %d, using '
+ 'nperseg = x.shape[%d]'
+ % (nperseg, axis, x.shape[axis], axis))
+ nperseg = x.shape[-1]
+
+ if isinstance(window, basestring) or type(window) is tuple:
+ win = get_window(window, nperseg)
+ else:
+ win = np.asarray(window)
+ if len(win.shape) != 1:
+ raise ValueError('window must be 1-D')
+ if win.shape[0] > x.shape[-1]:
+ raise ValueError('window is longer than x.')
+ nperseg = win.shape[0]
+
+ if scaling == 'density':
+ scale = 1.0 / (fs * (win*win).sum())
+ elif scaling == 'spectrum':
+ scale = 1.0 / win.sum()**2
+ else:
+ raise ValueError('Unknown scaling: %r' % scaling)
+
+ if noverlap is None:
+ noverlap = nperseg // 2
+ elif noverlap >= nperseg:
+ raise ValueError('noverlap must be less than nperseg.')
+
+ if nfft is None:
+ nfft = nperseg
+ elif nfft < nperseg:
+ raise ValueError('nfft must be greater than or equal to nperseg.')
+
+ if not hasattr(detrend, '__call__'):
+ detrend_func = lambda seg: signaltools.detrend(seg, type=detrend)
+ elif axis != -1:
+ # Wrap this function so that it receives a shape that it could
+ # reasonably expect to receive.
+ def detrend_func(seg):
+ seg = np.rollaxis(seg, -1, axis)
+ seg = detrend(seg)
+ return np.rollaxis(seg, axis, len(seg.shape))
+ else:
+ detrend_func = detrend
+
+ step = nperseg - noverlap
+ indices = np.arange(0, x.shape[-1]-nperseg+1, step)
+
+ if np.isrealobj(x) and return_onesided:
+ outshape = list(x.shape)
+ if nfft % 2 == 0: # even
+ outshape[-1] = nfft // 2 + 1
+ Pxx = np.empty(outshape, x.dtype)
+ for k, ind in enumerate(indices):
+ x_dt = detrend_func(x[..., ind:ind+nperseg])
+ xft = fftpack.rfft(x_dt*win, nfft)
+ # fftpack.rfft returns the positive frequency part of the fft
+ # as real values, packed r r i r i r i ...
+ # this indexing is to extract the matching real and imaginary
+ # parts, while also handling the pure real zero and nyquist
+ # frequencies.
+ if k == 0:
+ Pxx[..., (0,-1)] = xft[..., (0,-1)]**2
+ Pxx[..., 1:-1] = xft[..., 1:-1:2]**2 + xft[..., 2::2]**2
+ else:
+ Pxx *= k/(k+1.0)
+ Pxx[..., (0,-1)] += xft[..., (0,-1)]**2 / (k+1.0)
+ Pxx[..., 1:-1] += (xft[..., 1:-1:2]**2 + xft[..., 2::2]**2) \
+ / (k+1.0)
+ else: # odd
+ outshape[-1] = (nfft+1) // 2
+ Pxx = np.empty(outshape, x.dtype)
+ for k, ind in enumerate(indices):
+ x_dt = detrend_func(x[..., ind:ind+nperseg])
+ xft = fftpack.rfft(x_dt*win, nfft)
+ if k == 0:
+ Pxx[..., 0] = xft[..., 0]**2
+ Pxx[..., 1:] = xft[..., 1::2]**2 + xft[..., 2::2]**2
+ else:
+ Pxx *= k/(k+1.0)
+ Pxx[..., 0] += xft[..., 0]**2 / (k+1)
+ Pxx[..., 1:] += (xft[..., 1::2]**2 + xft[..., 2::2]**2) \
+ / (k+1.0)
+
+ Pxx[..., 1:-1] *= 2*scale
+ Pxx[..., (0,-1)] *= scale
+ f = np.arange(Pxx.shape[-1]) * (fs/nfft)
+ else:
+ for k, ind in enumerate(indices):
+ x_dt = detrend_func(x[..., ind:ind+nperseg])
+ xft = fftpack.fft(x_dt*win, nfft)
+ if k == 0:
+ Pxx = (xft * xft.conj()).real
+ else:
+ Pxx *= k/(k+1.0)
+ Pxx += (xft * xft.conj()).real / (k+1.0)
+ Pxx *= scale
+ f = fftpack.fftfreq(nfft, 1.0/fs)
+
+ if axis != -1:
+ Pxx = np.rollaxis(Pxx, -1, axis)
+
+ return f, Pxx
+
278 scipy/signal/tests/test_spectral.py
View
@@ -1,11 +1,279 @@
-"""Unit tests for Lomb Scargle routines.
-"""
-
+import warnings
import numpy as np
from numpy.testing import assert_raises, assert_approx_equal, \
- assert_, run_module_suite
+ assert_, run_module_suite, TestCase,\
+ assert_allclose, assert_array_equal,\
+ assert_array_almost_equal_nulp
+from scipy import signal, fftpack
+from scipy.signal import periodogram, welch, lombscargle
+
+
+class TestPeriodogram(TestCase):
+ def test_real_onesided_even(self):
+ x = np.zeros(16)
+ x[0] = 1
+ f, p = periodogram(x)
+ assert_allclose(f, np.linspace(0, 0.5, 9))
+ q = np.ones(9)
+ q[0] = 0
+ q[-1] /= 2.0
+ q /= 8
+ assert_allclose(p, q)
+
+ def test_real_onesided_odd(self):
+ x = np.zeros(15)
+ x[0] = 1
+ f, p = periodogram(x)
+ assert_allclose(f, np.arange(8.0)/15.0)
+ q = np.ones(8)
+ q[0] = 0
+ q[-1] /= 2.0
+ q *= 2.0/15.0
+ assert_allclose(p, q, atol=1e-15)
+
+ def test_real_twosided(self):
+ x = np.zeros(16)
+ x[0] = 1
+ f, p = periodogram(x, return_onesided=False)
+ assert_allclose(f, fftpack.fftfreq(16, 1.0))
+ q = np.ones(16)/16.0
+ q[0] = 0
+ assert_allclose(p, q)
+
+ def test_real_spectrum(self):
+ x = np.zeros(16)
+ x[0] = 1
+ f, p = periodogram(x, scaling='spectrum')
+ g, q = periodogram(x, scaling='density')
+ assert_allclose(f, np.linspace(0, 0.5, 9))
+ assert_allclose(p, q/16.0)
+
+ def test_complex(self):
+ x = np.zeros(16, np.complex128)
+ x[0] = 1.0 + 2.0j
+ f, p = periodogram(x)
+ assert_allclose(f, fftpack.fftfreq(16, 1.0))
+ q = 5.0*np.ones(16)/16.0
+ q[0] = 0
+ assert_allclose(p, q)
+
+ def test_unk_scaling(self):
+ assert_raises(ValueError, periodogram, np.zeros(4, np.complex128),
+ scaling='foo')
+
+ def test_nd_axis_m1(self):
+ x = np.zeros(20, dtype=np.float64)
+ x = x.reshape((2,1,10))
+ x[:,:,0] = 1.0
+ f, p = periodogram(x)
+ assert_array_equal(p.shape, (2, 1, 6))
+ assert_array_almost_equal_nulp(p[0,0,:], p[1,0,:], 60)
+ f0, p0 = periodogram(x[0,0,:])
+ assert_array_almost_equal_nulp(p0[np.newaxis,:], p[1,:], 60)
+
+ def test_nd_axis_0(self):
+ x = np.zeros(20, dtype=np.float64)
+ x = x.reshape((10,2,1))
+ x[0,:,:] = 1.0
+ f, p = periodogram(x, axis=0)
+ assert_array_equal(p.shape, (6,2,1))
+ assert_array_almost_equal_nulp(p[:,0,0], p[:,1,0], 60)
+ f0, p0 = periodogram(x[:,0,0])
+ assert_array_almost_equal_nulp(p0, p[:,1,0])
+
+ def test_window_external(self):
+ x = np.zeros(16)
+ x[0] = 1
+ f, p = periodogram(x, 10, 'hanning')
+ win = signal.get_window('hanning', 16)
+ fe, pe = periodogram(x, 10, win)
+ assert_array_almost_equal_nulp(p, pe)
+ assert_array_almost_equal_nulp(f, fe)
+
+ def test_padded_fft(self):
+ x = np.zeros(16)
+ x[0] = 1
+ f, p = periodogram(x)
+ fp, pp = periodogram(x, nfft=32)
+ assert_allclose(f, fp[::2])
+ assert_allclose(p, pp[::2])
+ assert_array_equal(pp.shape, (17,))
+
+ def test_empty_input(self):
+ f, p = periodogram([])
+ assert_array_equal(f.shape, (0,))
+ assert_array_equal(p.shape, (0,))
+ for shape in [(0,), (3,0), (0,5,2)]:
+ f, p = periodogram(np.empty(shape))
+ assert_array_equal(f.shape, shape)
+ assert_array_equal(p.shape, shape)
+
+ def test_short_nfft(self):
+ x = np.zeros(18)
+ x[0] = 1
+ f, p = periodogram(x, nfft=16)
+ assert_allclose(f, np.linspace(0, 0.5, 9))
+ q = np.ones(9)
+ q[0] = 0
+ q[-1] /= 2.0
+ q /= 8
+ assert_allclose(p, q)
+
+ def test_nfft_is_xshape(self):
+ x = np.zeros(16)
+ x[0] = 1
+ f, p = periodogram(x, nfft=16)
+ assert_allclose(f, np.linspace(0, 0.5, 9))
+ q = np.ones(9)
+ q[0] = 0
+ q[-1] /= 2.0
+ q /= 8
+ assert_allclose(p, q)
+
+
+class TestWelch(TestCase):
+ def test_real_onesided_even(self):
+ x = np.zeros(16)
+ x[0] = 1
+ x[8] = 1
+ f, p = welch(x, nperseg=8)
+ assert_allclose(f, np.linspace(0, 0.5, 5))
+ assert_allclose(p, np.array([ 0.08333333, 0.15277778, 0.22222222,
+ 0.22222222, 0.11111111]))
+
+ def test_real_onesided_odd(self):
+ x = np.zeros(16)
+ x[0] = 1
+ x[8] = 1
+ f, p = welch(x, nperseg=9)
+ assert_allclose(f, np.arange(5.0)/9.0)
+ assert_allclose(p, np.array([ 0.15958226, 0.24193954, 0.24145223,
+ 0.24100919, 0.12188675]))
+
+ def test_real_twosided(self):
+ x = np.zeros(16)
+ x[0] = 1
+ x[8] = 1
+ f, p = welch(x, nperseg=8, return_onesided=False)
+ assert_allclose(f, fftpack.fftfreq(8, 1.0))
+ assert_allclose(p, np.array([ 0.08333333, 0.07638889, 0.11111111,
+ 0.11111111, 0.11111111, 0.11111111, 0.11111111, 0.07638889]))
+
+ def test_real_spectrum(self):
+ x = np.zeros(16)
+ x[0] = 1
+ x[8] = 1
+ f, p = welch(x, nperseg=8, scaling='spectrum')
+ assert_allclose(f, np.linspace(0, 0.5, 5))
+ assert_allclose(p, np.array([ 0.015625, 0.028645833333333332,
+ 0.041666666666666664, 0.041666666666666664, 0.020833333333333332]))
+
+ def test_complex(self):
+ x = np.zeros(16, np.complex128)
+ x[0] = 1.0 + 2.0j
+ x[8] = 1.0 + 2.0j
+ f, p = welch(x, nperseg=8)
+ assert_allclose(f, fftpack.fftfreq(8, 1.0))
+ assert_allclose(p, np.array([ 0.41666667, 0.38194444, 0.55555556,
+ 0.55555556, 0.55555556, 0.55555556, 0.55555556, 0.38194444]))
+
+ def test_unk_scaling(self):
+ assert_raises(ValueError, welch, np.zeros(4, np.complex128),
+ scaling='foo', nperseg=4)
+
+ def test_detrend_linear(self):
+ x = np.arange(10, dtype=np.float64)+0.04
+ f, p = welch(x, nperseg=10, detrend='linear')
+ assert_allclose(p, np.zeros_like(p), atol=1e-15)
+
+ def test_detrend_external(self):
+ x = np.arange(10, dtype=np.float64)+0.04
+ f, p = welch(x, nperseg=10,
+ detrend=lambda seg: signal.detrend(seg, type='l'))
+ assert_allclose(p, np.zeros_like(p), atol=1e-15)
+
+ def test_detrend_external_nd_m1(self):
+ x = np.arange(40, dtype=np.float64)+0.04
+ x = x.reshape((2,2,10))
+ f, p = welch(x, nperseg=10,
+ detrend=lambda seg: signal.detrend(seg, type='l'))
+ assert_allclose(p, np.zeros_like(p), atol=1e-15)
+
+ def test_detrend_external_nd_0(self):
+ x = np.arange(20, dtype=np.float64)+0.04
+ x = x.reshape((2,1,10))
+ x = np.rollaxis(x, 2, 0)
+ f, p = welch(x, nperseg=10, axis=0,
+ detrend=lambda seg: signal.detrend(seg, axis=0, type='l'))
+ assert_allclose(p, np.zeros_like(p), atol=1e-15)
+
+ def test_nd_axis_m1(self):
+ x = np.arange(20, dtype=np.float64)+0.04
+ x = x.reshape((2,1,10))
+ f, p = welch(x, nperseg=10)
+ assert_array_equal(p.shape, (2, 1, 6))
+ assert_array_almost_equal_nulp(p[0,0,:], p[1,0,:], 60)
+ f0, p0 = welch(x[0,0,:], nperseg=10)
+ assert_array_almost_equal_nulp(p0[np.newaxis,:], p[1,:], 60)
+
+ def test_nd_axis_0(self):
+ x = np.arange(20, dtype=np.float64)+0.04
+ x = x.reshape((10,2,1))
+ f, p = welch(x, nperseg=10, axis=0)
+ assert_array_equal(p.shape, (6,2,1))
+ assert_array_almost_equal_nulp(p[:,0,0], p[:,1,0], 60)
+ f0, p0 = welch(x[:,0,0], nperseg=10)
+ assert_array_almost_equal_nulp(p0, p[:,1,0])
+
+ def test_window_external(self):
+ x = np.zeros(16)
+ x[0] = 1
+ x[8] = 1
+ f, p = welch(x, 10, 'hanning', 8)
+ win = signal.get_window('hanning', 8)
+ fe, pe = welch(x, 10, win, 8)
+ assert_array_almost_equal_nulp(p, pe)
+ assert_array_almost_equal_nulp(f, fe)
+
+ def test_empty_input(self):
+ f, p = welch([])
+ assert_array_equal(f.shape, (0,))
+ assert_array_equal(p.shape, (0,))
+ for shape in [(0,), (3,0), (0,5,2)]:
+ f, p = welch(np.empty(shape))
+ assert_array_equal(f.shape, shape)
+ assert_array_equal(p.shape, shape)
+
+ def test_short_data(self):
+ x = np.zeros(8)
+ x[0] = 1
+ with warnings.catch_warnings():
+ warnings.simplefilter('ignore', UserWarning)
+ f, p = welch(x)
+
+ f1, p1 = welch(x, nperseg=8)
+ assert_allclose(f, f1)
+ assert_allclose(p, p1)
+
+ def test_window_long_or_nd(self):
+ with warnings.catch_warnings():
+ warnings.simplefilter('ignore', UserWarning)
+ assert_raises(ValueError, welch, np.zeros(4), 1, np.array([1,1,1,1,1]))
+ assert_raises(ValueError, welch, np.zeros(4), 1,
+ np.arange(6).reshape((2,3)))
+
+ def test_nondefault_noverlap(self):
+ x = np.zeros(64)
+ x[::8] = 1
+ f, p = welch(x, nperseg=16, noverlap=4)
+ q = np.array([0, 1./12., 1./3., 1./5., 1./3., 1./5., 1./3., 1./5., 1./6.])
+ assert_allclose(p, q, atol=1e-12)
+
+ def test_bad_noverlap(self):
+ assert_raises(ValueError, welch, np.zeros(4), 1, 'hanning', 2, 7)
-from scipy.signal.spectral import lombscargle
+ def test_nfft_too_short(self):
+ assert_raises(ValueError, welch, np.ones(12), nfft=3, nperseg=4)
class TestLombscargle:
3  tools/py3tool.py
View
@@ -179,6 +179,7 @@ def custom_mangling(filename):
os.path.join('signal', '__init__.py'),
os.path.join('signal', 'bsplines.py'),
os.path.join('signal', 'signaltools.py'),
+ os.path.join('signal', 'spectral.py'),
os.path.join('signal', 'fir_filter_design.py'),
os.path.join('special', '__init__.py'),
os.path.join('special', 'add_newdocs.py'),
@@ -214,7 +215,7 @@ def custom_mangling(filename):
'_ufuncs', '_ufuncs_cxx',
'_minpack', '_zeros', '_lbfgsb', '_cobyla', '_slsqp',
'_nnls',
- 'sigtools', 'spline', 'spectral',
+ 'sigtools', 'spline', '_spectral',
'_fitpack', 'dfitpack', '_interpolate',
'_odepack', '_quadpack', 'vode', '_dop', 'lsoda',
'vonmises_cython', '_rank',
Something went wrong with that request. Please try again.