From 49189c3f4890e99032210442c405709aeefd7b03 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fran=C3=A7ois=20Boulogne?= Date: Tue, 23 Jul 2024 10:56:29 +0200 Subject: [PATCH] init --- .gitignore | 132 +++++++++++++++++++++++++++++++++++++++ README.md | 2 + mysignal/__init__.py | 1 + mysignal/mysignal.py | 143 +++++++++++++++++++++++++++++++++++++++++++ requirements.txt | 3 + setup.py | 72 ++++++++++++++++++++++ 6 files changed, 353 insertions(+) create mode 100644 .gitignore create mode 100644 README.md create mode 100644 mysignal/__init__.py create mode 100644 mysignal/mysignal.py create mode 100644 requirements.txt create mode 100644 setup.py diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..a95bc6e --- /dev/null +++ b/.gitignore @@ -0,0 +1,132 @@ +# Byte-compiled / optimized / DLL files +__pycache__/ +*.py[cod] +*$py.class + +# C extensions +*.so + +# Distribution / packaging +.Python +build/ +develop-eggs/ +dist/ +downloads/ +eggs/ +.eggs/ +lib/ +lib64/ +parts/ +sdist/ +var/ +wheels/ +*.egg-info/ +.installed.cfg +*.egg + +# PyInstaller +# Usually these files are written by a python script from a template +# before PyInstaller builds the exe, so as to inject date/other infos into it. +*.manifest +*.spec + +# Installer logs +pip-log.txt +pip-delete-this-directory.txt + +# Unit test / coverage reports +htmlcov/ +.tox/ +.nox/ +.coverage +.coverage.* +.cache +nosetests.xml +coverage.xml +*.cover +*.py,cover +.hypothesis/ +.pytest_cache/ +.coverage.* + +# Translations +*.mo +*.pot + +# Django stuff: +*.log +local_settings.py +db.sqlite3 + +# Flask stuff: +instance/ +.webassets-cache + +# Scrapy stuff: +.scrapy + +# Sphinx documentation +docs/_build/ +docs/autosummary/ + +# Jupyter Notebook +.ipynb_checkpoints +*.ipynb + +# IPython +profile_default/ +ipython_config.py + +# pyenv +# For a library or package, you might want to ignore these files since the code is +# intended to run in multiple environments; see https://github.com/pyenv/pyenv/issues/528 +.python-version + +# PEP 582; used by e.g. github.com/David-OConnor/pyflow and https://pdm.fming.dev +__pypackages__/ + +# Celery stuff +celerybeat-schedule +celerybeat.pid + +# SageMath parsed files +*.sage.py + +# Environments +.env +.venv +env/ +venv/ +ENV/ +env.bak/ +venv.bak/ +env.d/ +venv.d/ + +# Spyder project settings +.spyderproject +.spyproject + +# Rope project settings +.ropeproject + +# mkdocs documentation +/site + +# mypy +.mypy_cache/ +.dmypy.json +dmypy.json + +# Pyre type checker +.pyre/ + +# pytype static type analyzer +.pytype/ + +# Cython debug symbols +cython_debug/ + +# VSCode settings +.vscode/ + diff --git a/README.md b/README.md new file mode 100644 index 0000000..d98dc01 --- /dev/null +++ b/README.md @@ -0,0 +1,2 @@ +mysignal +======= diff --git a/mysignal/__init__.py b/mysignal/__init__.py new file mode 100644 index 0000000..b794fd4 --- /dev/null +++ b/mysignal/__init__.py @@ -0,0 +1 @@ +__version__ = '0.1.0' diff --git a/mysignal/mysignal.py b/mysignal/mysignal.py new file mode 100644 index 0000000..299a8db --- /dev/null +++ b/mysignal/mysignal.py @@ -0,0 +1,143 @@ +import functools +import numpy as np +import pandas as pd +from scipy.fft import fft, fftfreq +from scipy.signal import find_peaks +from scipy.signal import correlate +import matplotlib.pyplot as plt + + +def normalize_signal(s): + """ + Set the mean at 0, and max at 1. + """ + s -= s.mean() + s /= s.max() + return s + + +def get_global_mode(modes): + """ + For a series of modes, compute the period of the resulting signal. + + Ex (1, 2, 3) -> 6 + """ + return functools.reduce(np.lcm, modes) + + +def sync_phase_signal_in(df_, modes, baseperiod): + df = df_.copy() + dt = np.mean(np.diff(df['time'])) + N = len(df['time']) + t = np.linspace(0., N*dt, N, endpoint=False) + global_mode = get_global_mode(modes) + + signal_ref = 0 + for mode in modes: + signal_ref += np.sin(2 * np.pi * t / (mode * baseperiod)) + correlation = correlate(df['RHset'], signal_ref, mode='full', method='direct') + corr_peaks = find_peaks(correlation, distance=.95 * baseperiod * global_mode)[0] + x_corr = np.arange(0, len(correlation), 1) - (N-4) + # Take the first positive peak + corr_idx = corr_peaks[corr_peaks>np.argwhere(x_corr>0)[0][0]][0] + idx = x_corr[corr_idx] + + df = df[idx:] + df['time'] -= df['time'].iloc[0] + df = df.reset_index() + del df['index'] + return df + + +def resample_df(df_, modes, baseperiod, id_t='time', id_in='RH', id_out='dm_m'): + + t_origin = df_[id_t].to_numpy() + in_origin = df_[id_in].to_numpy() + out_origin = df_[id_out].to_numpy() + + + dt = np.mean(np.diff(t_origin)) + + + signal_period = get_global_mode(modes) * baseperiod + alpha = np.floor(len(t_origin) * dt / signal_period) + print(alpha) + + N = int(np.round(alpha * signal_period / dt)) + print(f'N_max: {len(t_origin)}, N: {N}') + # reconstructed time (strictly evenly spaced) + t = np.linspace(0., N*dt, N, endpoint=False) + signal_in = in_origin[:N] + signal_out = out_origin[:N] + return t, dt, signal_in, signal_out + + +def find_prominent_frequencies(signal, sampling_rate, num_frequencies=1): + # Compute the FFT of the signal + fft_result = fft(signal) + + # Compute the power spectrum + power_spectrum = np.abs(fft_result) ** 2 + + # Generate the frequency axis + freqs = np.fft.fftfreq(len(signal), d=1/sampling_rate) + + # Only consider the positive frequencies + positive_freqs = freqs[:len(freqs) // 2] + positive_power_spectrum = power_spectrum[:len(power_spectrum) // 2] + + # Find the indices of the peaks in the power spectrum + peaks, _ = find_peaks(positive_power_spectrum) + + # Get the indices of the highest peaks + highest_peaks = np.argsort(positive_power_spectrum[peaks])[-num_frequencies:] + + # Get the corresponding frequencies of the highest peaks + prominent_frequencies = positive_freqs[peaks][highest_peaks] + + return prominent_frequencies + + +def get_IQ_coeff(signal, sampling_rate, frequencies, debug=False): + N = len(signal) + + # if not hasattr(frequencies, '__iter__'): + # frequencies = (frequencies,) + + # FFT + Yf = fft(signal)[:N//2] + xf = fftfreq(N, 1/sampling_rate)[:N//2] + + # Get coefficients + Bcos_fft = [] + Bsin_fft = [] + + if debug: + fig, ax = plt.subplots(nrows=2, figsize=(8, 6)) + + for f in frequencies: + idx = np.argmin(np.abs(xf - f)) + print(f'distance xf f: {np.abs(xf[idx] - f)/f}') + #print('position freq:') + #print(idx, len(xf), idx / len(xf)) + Y_f = Yf[idx] + Bcos_fft.append(2.0 * Y_f.real / N) + Bsin_fft.append(-2.0 * Y_f.imag / N) + + if debug: + ax[0].plot(xf[idx-10:idx+10], Yf[idx-10:idx+10].real, '+-'); + ax[0].plot(xf[idx], Yf[idx].real, 'o', label=f'f={f:.1e}'); + ax[1].plot(xf[idx-10:idx+10], Yf[idx-10:idx+10].imag, 'x-'); + ax[1].plot(xf[idx], Yf[idx].imag, 'o', label=f'f={f:.1e}'); + + if debug: + plt.tight_layout() + ax[0].set_title('cos') + ax[0].legend() + ax[1].set_title('sin') + ax[1].legend() + + # Display + for i, f in enumerate(frequencies): + print(f"Freq {f:.2e} Hz, Period {1/f:.1f} s: A_I = {Bsin_fft[i]:.2e}, A_Q = {Bcos_fft[i]:.2e}, tan delta = {Bcos_fft[i]/Bsin_fft[i]:.2f}") + return np.array(Bcos_fft), np.array(Bsin_fft) diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..475fbd5 --- /dev/null +++ b/requirements.txt @@ -0,0 +1,3 @@ +numpy +scipy +pandas diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..b9b32fd --- /dev/null +++ b/setup.py @@ -0,0 +1,72 @@ + +import os +import re +from setuptools import setup, find_packages + +def read_version(name): + with open(os.path.join(os.path.dirname(__file__), name, '__init__.py')) as f: + for line in f: + match = re.match(r"^__version__ = ['\"]([^'\"]*)['\"]", line) + if match: + return match.group(1) + raise RuntimeError("Unable to find version string.") + +def parse_requirements(filename): + with open(filename, 'r') as file: + return file.read().splitlines() + +# Get the long description from the README file +here = os.path.abspath(os.path.dirname(__file__)) +with open(os.path.join(here, 'README.md'), encoding='utf-8') as f: + long_description = f.read() + + +setup( + name='mysignal', # Required + version = read_version('mysignal'), + description='A short description of the project', # Optional + long_description=long_description, # Optional + long_description_content_type='text/markdown', # Optional (see note above) + url='https://github.com/yourusername/yourproject', # Optional + author='Your Name', # Optional + author_email='your.email@example.com', # Optional + + # Automatically find packages in the current directory + packages=find_packages(), # Required + + # Include additional files specified in MANIFEST.in + include_package_data=True, # Optional + + # Define your dependencies in a separate file + install_requires=parse_requirements('requirements.txt'), # Optional + + # Classifiers help users find your project by categorizing it + classifiers=[ + 'Development Status :: 3 - Alpha', + 'Intended Audience :: Developers', + 'Topic :: Software Development :: Build Tools', + 'License :: OSI Approved :: MIT License', + 'Programming Language :: Python :: 3', + 'Programming Language :: Python :: 3.6', + 'Programming Language :: Python :: 3.7', + 'Programming Language :: Python :: 3.8', + 'Programming Language :: Python :: 3.9', + ], + + # Specify the Python versions you support + python_requires='>=3.6, <4', + + # Specify additional groups of dependencies (e.g., for testing) + # extras_require={ + # 'dev': ['check-manifest'], + # 'test': ['coverage'], + # }, + + # You can also specify entry points for command-line scripts + # entry_points={ + # 'console_scripts': [ + # 'sample=sample:main', + # ], + # }, +) +