claudiodsf / stockwell Goto Github PK

View Code? Open in Web Editor NEW

78.0 8.0 24.0 302 KB

Stockwell transform for Python

License: GNU General Public License v3.0

Shell 0.36% Python 90.77% C 8.57% Batchfile 0.30%

signal processing transform time-frequency-analysis

stockwell's Introduction

Stockwell

Python package for time-frequency analysis through Stockwell transform.

Based on original code from NIMH MEG Core Facility.

Installation

Using Anaconda

If you use Anaconda, the latest release of Stockwell is available via conda-forge.

To install, simply run:

conda install -c conda-forge stockwell

Using pip and PyPI

The latest release of Stockwell is available on the Python Package Index.

You can install it easily through pip:

pip install stockwell

Installation from source

If no precompiled package is available for you architecture on PyPI, or if you want to work on the source code, you will need to compile this package from source.

To obtain the source code, download the latest release from the releases page, or clone the GitHub project.

C compiler

Part of Stockwell is written in C, so you will need a C compiler.

On Linux (Debian or Ubuntu), install the build-essential package:

sudo apt install build-essential

On macOS, install the XCode Command Line Tools:

xcode-select --install

On Windows, install the Microsoft C++ Build Tools.

FFTW

To compile Stockwell, you will need to have FFTW installed.

If you use Anaconda (Linux, macOS, Windows):

conda install fftw

If you use Homebrew (macOS)

brew install fftw

If you use apt (Debian or Ubuntu)

sudo apt install libfftw3-dev

Install the Python package from source

Finally, install this Python package using pip:

pip install .

Or, alternatively, in "editable" mode:

pip install -e .

Usage

Example usage:

import numpy as np
from scipy.signal import chirp
import matplotlib.pyplot as plt
from stockwell import st

t = np.linspace(0, 10, 5001)
w = chirp(t, f0=12.5, f1=2.5, t1=10, method='linear')

fmin = 0  # Hz
fmax = 25  # Hz
df = 1./(t[-1]-t[0])  # sampling step in frequency domain (Hz)
fmin_samples = int(fmin/df)
fmax_samples = int(fmax/df)
stock = st.st(w, fmin_samples, fmax_samples)
extent = (t[0], t[-1], fmin, fmax)

fig, ax = plt.subplots(2, 1, sharex=True)
ax[0].plot(t, w)
ax[0].set(ylabel='amplitude')
ax[1].imshow(np.abs(stock), origin='lower', extent=extent)
ax[1].axis('tight')
ax[1].set(xlabel='time (s)', ylabel='frequency (Hz)')
plt.show()

You should get the following output:

You can also compute the inverse Stockwell transform, ex:

inv_stock = st.ist(stock, fmin_samples, fmax_samples)
fig, ax = plt.subplots(2, 1, sharex=True)
ax[0].plot(t, w, label='original signal')
ax[0].plot(t, inv_stock, label='inverse Stockwell')
ax[0].set(ylabel='amplitude')
ax[0].legend(loc='upper right')
ax[1].plot(t, w - inv_stock)
ax[1].set_xlim(0, 10)
ax[1].set(xlabel='time (s)', ylabel='amplitude difference')
plt.show()

References

Stockwell, R.G., Mansinha, L. & Lowe, R.P., 1996. Localization of the complex spectrum: the S transform, IEEE Trans. Signal Process., 44(4), 998–1001, doi:10.1109/78.492555

S transform on Wikipedia.

stockwell's People

Contributors

Stargazers

Watchers

stockwell's Issues

Installing problems in Python 3.9.12 Numpy 1.21.5

Hi,
I cannot install Stockwell in Linux with Python 3.9.12 Numpy 1.21.5.
To reproduce the error, download the code, extract and then

 xxx@xxx: pip install .
Processing ~/Downloads/stockwell-master
  DEPRECATION: A future pip version will change local packages to be built in-place without first copying to a temporary directory. We recommend you use --use-feature=in-tree-build to test your packages with this new behavior before it becomes the default.
   pip 21.3 will remove support for this functionality. You can find discussion regarding this at https://github.com/pypa/pip/issues/7555.
  Installing build dependencies ... done
  Getting requirements to build wheel ... done
    Preparing wheel metadata ... done
Requirement already satisfied: numpy>=1.18 in ~/miniconda3/lib/python3.9/site-packages (from stockwell==0+unknown) (1.21.5)
Building wheels for collected packages: stockwell
  Building wheel for stockwell (PEP 517) ... done
  Created wheel for stockwell: filename=stockwell-0+unknown-cp39-cp39-linux_x86_64.whl size=30531 sha256=3668038edf9336ccb613e1dbfd7071e10b09a673d1800140ca5f2b3de2e25e00
  Stored in directory: ~/.cache/pip/wheels/90/b6/62/b90786979ab2e9d77c3e53be01e275394e6f616accf04d0535
Successfully built stockwell
Installing collected packages: stockwell
Successfully installed stockwell-0+unknown

xxx@xxx:~/Downloads/stockwell-master$ python
Python 3.9.12 (main, Apr  5 2022, 06:56:58) 
[GCC 7.5.0] :: Anaconda, Inc. on linux
Type "help", "copyright", "credits" or "license" for more information.    
>>> from stockwell import st
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "~/Downloads/stockwell-master/stockwell/__init__.py", line 2, in <module>
    from .lib import st  #NOQA
ImportError: cannot import name 'st' from 'stockwell.lib' (~/Downloads/stockwell-master/stockwell/lib/__init__.py)

Thank you!
Regards,
Ignacio

PS: I have the same problem using pip install stockwell

does the gamma and window paprameters are included in the stoackwell transfrom ?

gamma = 1.3
window = 'gauss' # or 'kazemi'
st.st(data, window, gamma, 100, 1800)

Does it included in the stockwell transfrom .

when i implement something like this stock = st.st(w, window, gamma, fmin_samples, fmax_samples)
I am getting error as
TypeError: function takes at most 3 arguments (5 given)

Collaborate with a similar package

Hello, @claudiodsf ! I was recently able to open-source a package similar to this one that I've been using for particle motion analysis. I call it Particleman. The root of the package is also the Stockwell transform, and I've added a module implementing the normalized inner product of Meza-Fajardo et al., 2015 for the extraction of surface waves.

Your package has some great features, like a command-line script for reading/processing/plotting data, and I'm interested in collaborating. Is this something you're interested in doing? Thank you!

-Jon

Hilbert phase

Would you consider making the small change in the C code found in this github repository? The author states that he "...made a small change to avoid shifting the Hilbert transform in frequency when multiplying by the Gaussian window. Avoiding this shift means the normal Hilbert instaneous phase comes out in the spectrum." This is useful for doing phase-based analysis with the ST.

Power spectrum from Fourier and Stockwell transforms

Due to the normalized Gaussian weight, Stockwell integral in time (dx) must produce the Fourier spectrum (frequency domain) exactly:

F(f) = integral { S(x,f) . dx }

Therefore the power spectrum of Fourier transform must be the same as power spectrum of Stockwell transform (following your example):

dx = t[1] - t[0]
fk = np.fft.rfft(w) / (len(w))# normalized FFT
pk = np.abs(fk)**2
freq = np.fft.rfftfreq( len(w) , dx )
stock_p = (np.abs( (stock*dx).sum(axis=1)))**2

these two power spectrum are almost the same but meet each other only using a scale factor of 400!

scale = 400 # 400? 1?
fig,ax = plt.subplots(figsize=(8,8))
plt.plot(freq[0:251],stock_p/scale,linewidth=3, label='Stockwell transform')
plt.plot(freq[0:251],pk[0:251],'*',linewidth=3,color='k', label='conventional Fourier transform')
ax.set(xlabel='frequency (Hz)', ylabel='Power spectrum P(f)')
leg = plt.legend( prop={'size':10})

The scaling factor differs from case to case (e.g. Zimin et al, 2003 test needs the scale factor of ca. 640).

Do I forget something?!

Another question is concerning the periodicity. Does this library consider the signal as a periodic signal? what about the edges in non-periodic signals? Do we need something like "cone of influence"?!

Thanks in advance for your kind answer.

Python 3 can't import st

I'm probably forgetting something obvious, but when I start Python 3 in the main directory, it can't import st from stockwell:

>>> from stockwell import st
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/bubbarich/Documents/GitHub/stockwell/stockwell/__init__.py", line 1, in <module>
    from .lib import st
ImportError: cannot import name 'st'

I'm not familiar enough with importing from externally-defined namespaces to know how it is supposed to work. I tried several variations, and building and installing the project again.

Installing stockwell in Windows

Please provide instructions to install stockwell in Windows

Stockwell and Obspy

Hello, I have a problem when I try to use obspy stockwell with obspy. The problem exactly is:

RuntimeError Traceback (most recent call last)
RuntimeError: module compiled against API version 0xf but this version of numpy is 0xe

ImportError Traceback (most recent call last)
/tmp/ipykernel_101293/699505916.py in
3 from obspy import read
4 from scipy.signal import hilbert
----> 5 from stockwell import st

~/anaconda3/envs/obspy/lib/python3.9/site-packages/stockwell/init.py in
1 """Time-frequency analysis through Stockwell transform."""
----> 2 from .lib import st #NOQA
3 from .lib import sine #NOQA
4 from ._version import get_versions
5 version = get_versions()['version']

ImportError: numpy.core.multiarray failed to import

This problem, I think is a problem with numpy version. I try to solve this installing version 1.23.0, but obspy doesn't work.

Do you have some options?

Thanks in advance

Joan

Hi buddy，is there an official documentation?

Recently I've been studying the S-transform，prepare to use the S-transform in the my paper,I used it in Matlab，but I want to do it in Python，could you provide API documentation for Python?Thx