commit 4dd36f8dc20dbd8146ab3901a0b637f8d10b657d Author: Benoît Sierro Date: Wed Apr 6 13:41:07 2022 +0200 initial commit diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..7e12e49 --- /dev/null +++ b/.gitignore @@ -0,0 +1,163 @@ +.vscode +.DS_Store + +# 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/ +share/python-wheels/ +*.egg-info/ +.installed.cfg +*.egg +MANIFEST + +# 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/ +cover/ + +# Translations +*.mo +*.pot + +# Django stuff: +*.log +local_settings.py +db.sqlite3 +db.sqlite3-journal + +# Flask stuff: +instance/ +.webassets-cache + +# Scrapy stuff: +.scrapy + +# Sphinx documentation +docs/_build/ + +# PyBuilder +.pybuilder/ +target/ + +# Jupyter Notebook +.ipynb_checkpoints + +# 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; otherwise, check them in: +# .python-version + +# pipenv +# According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control. +# However, in case of collaboration, if having platform-specific dependencies or dependencies +# having no cross-platform support, pipenv may install dependencies that don't work, or not +# install all needed dependencies. +#Pipfile.lock + +# poetry +# Similar to Pipfile.lock, it is generally recommended to include poetry.lock in version control. +# This is especially recommended for binary packages to ensure reproducibility, and is more +# commonly ignored for libraries. +# https://python-poetry.org/docs/basic-usage/#commit-your-poetrylock-file-to-version-control +#poetry.lock + +# pdm +# Similar to Pipfile.lock, it is generally recommended to include pdm.lock in version control. +#pdm.lock +# pdm stores project-wide configurations in .pdm.toml, but it is recommended to not include it +# in version control. +# https://pdm.fming.dev/#use-with-ide +.pdm.toml + +# PEP 582; used by e.g. github.com/David-OConnor/pyflow and github.com/pdm-project/pdm +__pypackages__/ + +# Celery stuff +celerybeat-schedule +celerybeat.pid + +# SageMath parsed files +*.sage.py + +# Environments +.env +.venv +env/ +venv/ +ENV/ +env.bak/ +venv.bak/ + +# 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/ + +# PyCharm +# JetBrains specific template is maintained in a separate JetBrains.gitignore that can +# be found at https://github.com/github/gitignore/blob/main/Global/JetBrains.gitignore +# and can be added to the global gitignore or merged into this file. For a more nuclear +# option (not recommended) you can uncomment the following to ignore the entire idea folder. +#.idea/ \ No newline at end of file diff --git a/data/Argon_750nm.CSV b/data/Argon_750nm.CSV new file mode 100755 index 0000000..7fdef79 --- /dev/null +++ b/data/Argon_750nm.CSV @@ -0,0 +1,530 @@ +70CSV +// AQ6370 OPTICAL SPECTRUM ANALYZER // +25 +"CTRWL", 750.180000 +"SPAN", 2.000000 +"START WL", 749.180000 +"STOP WL", 751.180000 +"WLFREQ", 0 +"REFL", -69.900 +"LSCL",-3.0 +"BASEL",0.00000000000 +"RESLN",0.020 +"AVG", 10 +"SMPLAUTO", 1 +"SMPL", 501 +"SMPLINTVL",0.0040 +"HIGH3" +"MEAS" +"LSUNT",0 +"NMSKH","OFF" +"RESCOR",0 +"RESPARM",14283 +"FREQPARM",14283 + + + + + +"[TRACE DATA]" + 749.1800,1.886E-009 + 749.1840,-1.814E-010 + 749.1880,4.393E-009 + 749.1920,5.302E-009 + 749.1960,1.564E-009 + 749.2000,2.578E-009 + 749.2040,3.173E-009 + 749.2080,2.543E-009 + 749.2120,4.090E-009 + 749.2160,2.254E-009 + 749.2200,1.613E-009 + 749.2240,2.656E-009 + 749.2280,4.871E-009 + 749.2320,2.881E-009 + 749.2360,2.523E-009 + 749.2400,1.763E-009 + 749.2440,5.209E-009 + 749.2480,3.911E-009 + 749.2520,3.312E-009 + 749.2560,2.970E-009 + 749.2600,3.186E-009 + 749.2640,3.533E-009 + 749.2680,6.071E-009 + 749.2720,6.027E-009 + 749.2760,2.802E-009 + 749.2800,4.396E-009 + 749.2840,4.363E-009 + 749.2880,4.154E-009 + 749.2920,5.001E-009 + 749.2960,6.928E-009 + 749.3000,3.123E-009 + 749.3040,4.637E-009 + 749.3080,3.466E-009 + 749.3120,4.200E-009 + 749.3160,3.010E-009 + 749.3200,4.786E-009 + 749.3240,5.492E-009 + 749.3280,3.710E-009 + 749.3320,5.558E-009 + 749.3360,3.259E-009 + 749.3400,4.192E-009 + 749.3440,1.736E-009 + 749.3480,2.128E-009 + 749.3520,5.626E-009 + 749.3560,1.442E-009 + 749.3600,2.901E-009 + 749.3640,4.523E-009 + 749.3680,4.928E-009 + 749.3720,5.645E-009 + 749.3760,4.817E-009 + 749.3800,4.527E-009 + 749.3840,4.506E-009 + 749.3880,3.123E-009 + 749.3920,3.111E-009 + 749.3960,2.992E-009 + 749.4000,5.482E-009 + 749.4040,4.208E-009 + 749.4080,6.201E-009 + 749.4120,5.391E-009 + 749.4160,4.743E-009 + 749.4200,5.388E-009 + 749.4240,5.527E-009 + 749.4280,7.186E-009 + 749.4320,4.669E-009 + 749.4360,5.038E-009 + 749.4400,4.892E-009 + 749.4440,2.480E-009 + 749.4480,2.981E-009 + 749.4520,2.782E-009 + 749.4560,5.581E-009 + 749.4600,6.543E-009 + 749.4640,6.392E-009 + 749.4680,7.106E-009 + 749.4720,3.108E-009 + 749.4760,5.250E-009 + 749.4800,3.147E-009 + 749.4840,6.431E-009 + 749.4880,4.606E-009 + 749.4920,4.045E-009 + 749.4960,6.164E-009 + 749.5000,3.298E-009 + 749.5040,4.967E-009 + 749.5080,5.383E-009 + 749.5120,4.850E-009 + 749.5160,5.807E-009 + 749.5200,5.329E-009 + 749.5240,7.289E-009 + 749.5280,5.368E-009 + 749.5320,6.169E-009 + 749.5360,1.040E-008 + 749.5400,4.093E-009 + 749.5440,5.866E-009 + 749.5480,3.356E-009 + 749.5520,7.623E-009 + 749.5560,4.415E-009 + 749.5600,3.819E-009 + 749.5640,3.544E-009 + 749.5680,7.250E-009 + 749.5720,4.572E-009 + 749.5760,3.755E-009 + 749.5800,5.222E-009 + 749.5840,5.988E-009 + 749.5880,5.257E-009 + 749.5920,6.610E-009 + 749.5960,3.157E-009 + 749.6000,3.219E-009 + 749.6040,4.556E-009 + 749.6080,6.296E-009 + 749.6120,5.584E-009 + 749.6160,4.192E-009 + 749.6200,7.069E-009 + 749.6240,5.588E-009 + 749.6280,4.157E-009 + 749.6320,3.856E-009 + 749.6360,5.057E-009 + 749.6400,5.239E-009 + 749.6440,5.843E-009 + 749.6480,2.581E-009 + 749.6520,6.027E-009 + 749.6560,6.581E-009 + 749.6600,5.121E-009 + 749.6640,6.653E-009 + 749.6680,4.271E-009 + 749.6720,5.540E-009 + 749.6760,6.515E-009 + 749.6800,7.306E-009 + 749.6840,4.190E-009 + 749.6880,3.709E-009 + 749.6920,6.302E-009 + 749.6960,6.388E-009 + 749.7000,5.680E-009 + 749.7040,5.440E-009 + 749.7080,5.455E-009 + 749.7120,7.778E-009 + 749.7160,4.530E-009 + 749.7200,6.271E-009 + 749.7240,8.175E-009 + 749.7280,6.554E-009 + 749.7320,5.844E-009 + 749.7360,6.875E-009 + 749.7400,7.670E-009 + 749.7440,6.135E-009 + 749.7480,6.296E-009 + 749.7520,6.238E-009 + 749.7560,8.811E-009 + 749.7600,1.027E-008 + 749.7640,3.792E-009 + 749.7680,6.693E-009 + 749.7720,5.717E-009 + 749.7760,5.449E-009 + 749.7800,8.478E-009 + 749.7840,7.188E-009 + 749.7880,9.100E-009 + 749.7920,7.598E-009 + 749.7960,6.140E-009 + 749.8000,7.155E-009 + 749.8040,7.712E-009 + 749.8080,7.128E-009 + 749.8120,4.515E-009 + 749.8160,6.684E-009 + 749.8200,7.515E-009 + 749.8240,5.411E-009 + 749.8280,6.048E-009 + 749.8320,9.447E-009 + 749.8360,9.474E-009 + 749.8400,6.826E-009 + 749.8440,6.150E-009 + 749.8480,9.486E-009 + 749.8520,8.923E-009 + 749.8560,8.326E-009 + 749.8600,9.442E-009 + 749.8640,4.681E-009 + 749.8680,7.626E-009 + 749.8720,7.288E-009 + 749.8760,5.657E-009 + 749.8800,9.260E-009 + 749.8840,7.951E-009 + 749.8880,8.935E-009 + 749.8920,6.514E-009 + 749.8960,8.464E-009 + 749.9000,9.091E-009 + 749.9040,1.053E-008 + 749.9080,1.103E-008 + 749.9120,9.890E-009 + 749.9160,6.340E-009 + 749.9200,8.855E-009 + 749.9240,1.044E-008 + 749.9280,9.490E-009 + 749.9320,1.076E-008 + 749.9360,9.019E-009 + 749.9400,1.063E-008 + 749.9440,9.949E-009 + 749.9480,9.264E-009 + 749.9520,1.116E-008 + 749.9560,1.238E-008 + 749.9600,1.025E-008 + 749.9640,1.183E-008 + 749.9680,1.042E-008 + 749.9720,9.502E-009 + 749.9760,7.598E-009 + 749.9800,1.096E-008 + 749.9840,1.079E-008 + 749.9880,1.368E-008 + 749.9920,1.069E-008 + 749.9960,1.192E-008 + 750.0000,1.211E-008 + 750.0040,1.476E-008 + 750.0080,1.534E-008 + 750.0120,1.506E-008 + 750.0160,1.322E-008 + 750.0200,1.406E-008 + 750.0240,1.508E-008 + 750.0280,1.574E-008 + 750.0320,1.629E-008 + 750.0360,1.545E-008 + 750.0400,1.585E-008 + 750.0440,1.731E-008 + 750.0480,1.625E-008 + 750.0520,1.876E-008 + 750.0560,1.791E-008 + 750.0600,1.896E-008 + 750.0640,1.980E-008 + 750.0680,2.010E-008 + 750.0720,2.023E-008 + 750.0760,2.128E-008 + 750.0800,2.544E-008 + 750.0840,2.544E-008 + 750.0880,2.731E-008 + 750.0920,2.901E-008 + 750.0960,2.839E-008 + 750.1000,3.096E-008 + 750.1040,3.367E-008 + 750.1080,3.454E-008 + 750.1120,4.052E-008 + 750.1160,4.003E-008 + 750.1200,4.248E-008 + 750.1240,4.626E-008 + 750.1280,5.003E-008 + 750.1320,5.483E-008 + 750.1360,5.708E-008 + 750.1400,6.262E-008 + 750.1440,6.452E-008 + 750.1480,6.785E-008 + 750.1520,6.922E-008 + 750.1560,7.622E-008 + 750.1600,7.845E-008 + 750.1640,8.219E-008 + 750.1680,8.486E-008 + 750.1720,8.939E-008 + 750.1760,9.025E-008 + 750.1800,9.551E-008 + 750.1840,9.807E-008 + 750.1880,1.009E-007 + 750.1920,1.024E-007 + 750.1960,1.022E-007 + 750.2000,1.002E-007 + 750.2040,1.001E-007 + 750.2080,1.027E-007 + 750.2120,1.021E-007 + 750.2160,1.022E-007 + 750.2200,9.834E-008 + 750.2240,9.332E-008 + 750.2280,9.466E-008 + 750.2320,8.955E-008 + 750.2360,8.618E-008 + 750.2400,8.028E-008 + 750.2440,7.696E-008 + 750.2480,7.295E-008 + 750.2520,6.936E-008 + 750.2560,6.697E-008 + 750.2600,6.196E-008 + 750.2640,5.951E-008 + 750.2680,5.537E-008 + 750.2720,5.195E-008 + 750.2760,4.552E-008 + 750.2800,4.439E-008 + 750.2840,4.067E-008 + 750.2880,3.739E-008 + 750.2920,3.261E-008 + 750.2960,3.446E-008 + 750.3000,3.258E-008 + 750.3040,2.899E-008 + 750.3080,2.870E-008 + 750.3120,2.625E-008 + 750.3160,2.627E-008 + 750.3200,2.699E-008 + 750.3240,2.223E-008 + 750.3280,2.266E-008 + 750.3320,2.202E-008 + 750.3360,2.112E-008 + 750.3400,1.899E-008 + 750.3440,1.847E-008 + 750.3480,1.961E-008 + 750.3520,1.872E-008 + 750.3560,1.725E-008 + 750.3600,1.668E-008 + 750.3640,1.558E-008 + 750.3680,1.712E-008 + 750.3720,1.662E-008 + 750.3760,1.451E-008 + 750.3800,1.382E-008 + 750.3840,1.517E-008 + 750.3880,1.397E-008 + 750.3920,1.554E-008 + 750.3960,1.302E-008 + 750.4000,1.672E-008 + 750.4040,1.285E-008 + 750.4080,1.271E-008 + 750.4120,1.242E-008 + 750.4160,1.341E-008 + 750.4200,1.195E-008 + 750.4240,1.081E-008 + 750.4280,1.136E-008 + 750.4320,8.758E-009 + 750.4360,1.142E-008 + 750.4400,1.043E-008 + 750.4440,8.205E-009 + 750.4480,1.255E-008 + 750.4520,8.982E-009 + 750.4560,9.364E-009 + 750.4600,8.464E-009 + 750.4640,1.071E-008 + 750.4680,1.159E-008 + 750.4720,6.774E-009 + 750.4760,1.051E-008 + 750.4800,9.242E-009 + 750.4840,9.927E-009 + 750.4880,6.740E-009 + 750.4920,9.382E-009 + 750.4960,8.447E-009 + 750.5000,1.111E-008 + 750.5040,9.997E-009 + 750.5080,8.018E-009 + 750.5120,1.122E-008 + 750.5160,9.959E-009 + 750.5200,8.895E-009 + 750.5240,1.157E-008 + 750.5280,8.677E-009 + 750.5320,9.303E-009 + 750.5360,7.186E-009 + 750.5400,9.815E-009 + 750.5440,7.531E-009 + 750.5480,1.158E-008 + 750.5520,8.225E-009 + 750.5560,7.345E-009 + 750.5600,6.519E-009 + 750.5640,7.842E-009 + 750.5680,8.510E-009 + 750.5720,9.191E-009 + 750.5760,1.218E-008 + 750.5800,6.444E-009 + 750.5840,9.170E-009 + 750.5880,9.678E-009 + 750.5920,7.806E-009 + 750.5960,7.396E-009 + 750.6000,6.904E-009 + 750.6040,9.053E-009 + 750.6080,1.029E-008 + 750.6120,5.697E-009 + 750.6160,1.016E-008 + 750.6200,1.036E-008 + 750.6240,8.116E-009 + 750.6280,6.741E-009 + 750.6320,8.353E-009 + 750.6360,7.724E-009 + 750.6400,9.610E-009 + 750.6440,8.374E-009 + 750.6480,1.078E-008 + 750.6520,9.833E-009 + 750.6560,7.249E-009 + 750.6600,6.146E-009 + 750.6640,6.232E-009 + 750.6680,6.907E-009 + 750.6720,6.143E-009 + 750.6760,6.774E-009 + 750.6800,6.977E-009 + 750.6840,6.737E-009 + 750.6880,8.807E-009 + 750.6920,8.040E-009 + 750.6960,7.598E-009 + 750.7000,7.425E-009 + 750.7040,5.808E-009 + 750.7080,8.750E-009 + 750.7120,9.276E-009 + 750.7160,8.646E-009 + 750.7200,1.036E-008 + 750.7240,8.932E-009 + 750.7280,8.940E-009 + 750.7320,1.151E-008 + 750.7360,8.142E-009 + 750.7400,8.845E-009 + 750.7440,7.983E-009 + 750.7480,1.041E-008 + 750.7520,6.024E-009 + 750.7560,9.695E-009 + 750.7600,8.548E-009 + 750.7640,8.050E-009 + 750.7680,9.533E-009 + 750.7720,7.497E-009 + 750.7760,7.554E-009 + 750.7800,1.029E-008 + 750.7840,7.741E-009 + 750.7880,7.026E-009 + 750.7920,9.383E-009 + 750.7960,6.949E-009 + 750.8000,7.401E-009 + 750.8040,8.791E-009 + 750.8080,8.840E-009 + 750.8120,5.809E-009 + 750.8160,7.133E-009 + 750.8200,6.826E-009 + 750.8240,5.498E-009 + 750.8280,6.650E-009 + 750.8320,6.834E-009 + 750.8360,6.541E-009 + 750.8400,9.232E-009 + 750.8440,5.792E-009 + 750.8480,9.510E-009 + 750.8520,7.648E-009 + 750.8560,1.193E-008 + 750.8600,1.032E-008 + 750.8640,7.720E-009 + 750.8680,6.292E-009 + 750.8720,1.025E-008 + 750.8760,7.721E-009 + 750.8800,5.053E-009 + 750.8840,7.349E-009 + 750.8880,7.235E-009 + 750.8920,7.434E-009 + 750.8960,6.194E-009 + 750.9000,5.896E-009 + 750.9040,6.173E-009 + 750.9080,8.190E-009 + 750.9120,7.514E-009 + 750.9160,5.786E-009 + 750.9200,6.912E-009 + 750.9240,7.847E-009 + 750.9280,8.771E-009 + 750.9320,8.483E-009 + 750.9360,7.402E-009 + 750.9400,9.425E-009 + 750.9440,7.235E-009 + 750.9480,7.314E-009 + 750.9520,9.199E-009 + 750.9560,7.407E-009 + 750.9600,9.833E-009 + 750.9640,8.369E-009 + 750.9680,8.188E-009 + 750.9720,5.897E-009 + 750.9760,1.121E-008 + 750.9800,7.247E-009 + 750.9840,1.030E-008 + 750.9880,9.059E-009 + 750.9920,6.305E-009 + 750.9960,7.803E-009 + 751.0000,6.260E-009 + 751.0040,6.632E-009 + 751.0080,8.439E-009 + 751.0120,6.287E-009 + 751.0160,9.148E-009 + 751.0200,8.275E-009 + 751.0240,8.422E-009 + 751.0280,7.888E-009 + 751.0320,8.916E-009 + 751.0360,7.581E-009 + 751.0400,9.430E-009 + 751.0440,9.312E-009 + 751.0480,9.281E-009 + 751.0520,9.729E-009 + 751.0560,1.015E-008 + 751.0600,9.065E-009 + 751.0640,1.002E-008 + 751.0680,7.607E-009 + 751.0720,8.622E-009 + 751.0760,1.189E-008 + 751.0800,9.802E-009 + 751.0840,8.685E-009 + 751.0880,6.030E-009 + 751.0920,9.090E-009 + 751.0960,1.158E-008 + 751.1000,5.451E-009 + 751.1040,9.767E-009 + 751.1080,6.368E-009 + 751.1120,7.536E-009 + 751.1160,1.124E-008 + 751.1200,1.017E-008 + 751.1240,9.669E-009 + 751.1280,9.600E-009 + 751.1320,8.671E-009 + 751.1360,1.134E-008 + 751.1400,9.225E-009 + 751.1440,1.207E-008 + 751.1480,8.720E-009 + 751.1520,5.573E-009 + 751.1560,9.581E-009 + 751.1600,9.566E-009 + 751.1640,1.122E-008 + 751.1680,1.144E-008 + 751.1720,1.069E-008 + 751.1760,1.181E-008 + 751.1800,1.227E-008 diff --git a/data/data.png b/data/data.png new file mode 100644 index 0000000..e056d95 Binary files /dev/null and b/data/data.png differ diff --git a/deconvolution.py b/deconvolution.py new file mode 100644 index 0000000..72f87c1 --- /dev/null +++ b/deconvolution.py @@ -0,0 +1,400 @@ +from dataclasses import dataclass, field +from pathlib import Path +from typing import Callable + +import click +import labmodule as lab +import numpy as np +from customfunc.app import PlotApp +from numpy.fft import fft, fftshift, ifft +from scipy.optimize import curve_fit + +LN2 = np.log(2) + + +@click.group() +def main(): + pass + + +@dataclass +class FitResult: + curve: np.ndarray = field(repr=False) + covariance: np.ndarray + amp: float + width: float + x0: float = 0 + params: tuple[float, ...] = field(default_factory=tuple) + + +@dataclass +class VoigtFitResult: + name: str + curve: np.ndarray = field(repr=False) + covariance: np.ndarray + gaussian_width: float + lorentzian_width: float + + @property + def descr(self) -> str: + return ( + f"{self.name} - Gaussian width : {self.gaussian_width:.2f}, " + f"Lorentzian width : {self.lorentzian_width:.2f}" + ) + + +def setup_axes(n: int, lims=(-10, 10)) -> tuple[np.ndarray, float, np.ndarray, float]: + """Creates a time-like and a frequency-like array and returns their respective spacing""" + x, dx = np.linspace(*lims, n, retstep=True) + f = np.fft.fftfreq(n, dx) + return x, dx, f, f[1] - f[0] + + +def freq_axis(x: np.ndarray) -> tuple[np.ndarray, np.ndarray, float]: + """ + given an array with equally spaced values, returns a corresponding frequency array, + an array of indices to automatically sort ffts and the x spacing + """ + dx = np.diff(x).mean() + f = np.fft.fftfreq(len(x), dx) + ind = np.argsort(f) + return f[ind], ind, dx + + +def lorentzian(x: np.ndarray, amp: float, width: float, x0: float) -> np.ndarray: + """Lorentzian function of max `amp`, FWHM `width` and position `x0`""" + gamma = (width * 0.5) ** 2 + return amp * gamma / ((x - x0) ** 2 + gamma) + + +def unit_lorentzian(x: np.ndarray, width: float, x0: float) -> np.ndarray: + """Normalized Lorentzian function of FWHM `width` and position `x0`""" + return lorentzian(x, 1 / (np.pi * width * 0.5), width, x0) + + +def peak_func(x: np.ndarray, amp: float, width: float, x0: float) -> np.ndarray: + """2-sided decaying exp of max `amp`, FWHM `width` and position `x0`""" + a = 2 * LN2 / width + return amp * np.exp(-np.abs(x - x0) * a) + + +def fft_lorentzian(f: np.ndarray, dx: float, amp: float, width: float) -> np.ndarray: + """ + Same as `peak_func` but parametrized as the fourier transform + of a Lorentzian function of max `amp` and FWHM `width` + """ + gamma_sqrt = width / 2 + return amp * np.pi * gamma_sqrt * np.exp(-2 * np.pi * gamma_sqrt * np.abs(f)) / dx + + +def gaussian(x: np.ndarray, amp: float, width: float, x0: float) -> np.ndarray: + """Gaussian function of max `amp`, FWHM `width` and position `x0`""" + sig2 = width ** 2 / (8 * LN2) + return amp * np.exp(-((x - x0) ** 2) / (2 * sig2)) + + +def elevated_gaussian(x: np.ndarray, amp: float, width: float, x0: float, y0: float) -> np.ndarray: + """Gaussian function of max `amp`, FWHM `width`, position `x0` and baseline `y0`""" + sig2 = width ** 2 / (8 * LN2) + return amp * np.exp(-((x - x0) ** 2) / (2 * sig2)) + y0 + + +def unit_gaussian(x: np.ndarray, width: float, x0: float) -> np.ndarray: + """Normalized Gaussian function of FWHM `width` and position `x0`""" + return gaussian(x, np.sqrt(4 * LN2 / np.pi) / width, width, x0) + + +def fft_gaussian(f: np.ndarray, dx: float, amp: float, width: float) -> np.ndarray: + """ + Same as `gaussian` but parametrized as the fourier transform + of another Gaussian function of max `amp` and FWHM `width` + """ # exp_param = + freq_amp = amp * width / 2 * np.sqrt(np.pi / LN2) / dx + freq_width = 4 * LN2 / (width * np.pi) + return gaussian(f, freq_amp, freq_width, 0) + + +def voigt(x: np.ndarray, gaussian_width: float, lorentzian_width: float) -> np.ndarray: + """Voigt profile of max 1 with specified Gaussian and Lorentzian FWHM""" + conv = convolve(unit_gaussian(x, gaussian_width, 0), unit_lorentzian(x, lorentzian_width, 0)) + if conv.max() > 0: + return conv / conv.max() + return conv + + +def log_voigt(x: np.ndarray, gaussian_width: float, lorentzian_width: float) -> np.ndarray: + """log of `voigt`""" + return np.log(voigt(x, gaussian_width, lorentzian_width)) + + +def fft_voigt( + f: np.ndarray, dx: float, gaussian_width: float, lorentzian_width: float +) -> np.ndarray: + """ + returns the product of a `peak_func` and a `gaussian` + parameterized as the Fourier transform of a Voigt profile + given by `gaussian_width` and `lorentzian_width` + """ + profile = fft_gaussian(f, dx, 1, gaussian_width) * fft_lorentzian(f, dx, 1, lorentzian_width) + return profile / profile.max() + + +def single_func_fit( + x: np.ndarray, + y: np.ndarray, + fct: Callable[[np.ndarray, float, float, float], np.ndarray], + force_center=False, + use_weights=False, + extra_params: int = 0, +) -> FitResult: + """fits a single bell-like function (Gaussian or Lorentian, not Voigt) + + Parameters + ---------- + x : np.ndarray + x input + y : np.ndarray + corresponding y data + fct : Callable[[np.ndarray, float, float, float], np.ndarray] + function to fit + force_center : bool, optional + assume the bell is centered at x=0, by default False + use_weights : bool, optional + give more importance to higher y values, by default False + + Returns + ------- + FitResult + resuts of the fit + """ + if force_center: + + def _fct(_x: np.ndarray, amp: float, width: float): + return fct(_x, amp, width, 0) + + fct = _fct + am = y.argmax() + params, cov = curve_fit( + fct, x, y, sigma=1 / y if use_weights else None, p0=[y[am], 1, x[am]] + [0] * extra_params + ) + curve = fct(x, *params) + return FitResult(curve, cov, params[0], params[1], params[2], tuple(params)) + + +def linear_direct_voigt_fit(x: np.ndarray, y: np.ndarray, use_weights=False) -> VoigtFitResult: + params, cov = curve_fit(voigt, x, y, sigma=1 / y if use_weights else None) + curve = voigt(x, *params) + return VoigtFitResult("Linear direct fit", curve, cov, *params) + + +def log_direct_voigt_fit(x: np.ndarray, y: np.ndarray) -> VoigtFitResult: + ind = y > 0 + y_fit = np.log(y[ind]) + x_fit = x[ind] + params, cov = curve_fit(log_voigt, x_fit, y_fit, bounds=[1e-12, np.inf]) + curve = np.exp(log_voigt(x, *params)) + return VoigtFitResult("Log direct fit", curve, cov, *params) + + +def linear_fft_voigt_fit(f: np.ndarray, A: np.ndarray, x: np.ndarray, dx: float) -> VoigtFitResult: + params, cov = curve_fit( + lambda _f, _g, _l: fft_voigt(_f, dx, _g, _l), f, A, bounds=[1e-12, np.inf] + ) + curve = voigt(x, *params) + return VoigtFitResult("Linear FFT fit", curve, cov, *params) + + +def convolve(gau: np.ndarray, lor: np.ndarray) -> np.ndarray: + return np.convolve(gau, lor, mode="same") + + +def add_noise(arr: np.ndarray, abs_fac: float, rel_fac: float): + return ( + arr + abs_fac * np.random.rand(len(arr)) + rel_fac * arr * (np.random.rand(len(arr)) - 0.5) + ) + + +@main.command(help="fit a lorentzian with a gaussian") +def demo1(): + x = np.linspace(-10, 10, 501) + with PlotApp( + FWHM=np.geomspace(0.1, 10), + relative_noise=np.linspace(0, 2), + absolute_noise=[0, *np.geomspace(0.001, 1)], + ) as app: + app.set_antialiasing(True) + app.params["FWHM"].value = 3.7 + + @app.update + def draw(FWHM, relative_noise, absolute_noise): + lo = lorentzian(x, 1, FWHM, 0) + lo = add_noise(lo, absolute_noise, relative_noise * 0.04) + fit = single_func_fit(x, lo, gaussian) + app[0].set_line_data("Lorentzian", x, lo, width=2) + app[0].set_line_data( + "gaus fit", x, fit.curve, width=2, label=f"Gaussian fit : FWHM = {fit.width}" + ) + + +@main.command(name="default", help="illustrate the influence of noises and fit a voigt profile") +def demo2(): + with PlotApp( + n=np.arange(256, 2046), + lorentz_fwhm=np.geomspace(0.1, 1), + gauss_fwhm=np.linspace(1, 10), + relative_noise=np.linspace(0, 2), + absolute_noise=[0, *np.geomspace(1e-6, 1, 201)], + ) as app: + app.set_antialiasing(True) + ax = app["Main Plot"] + + app.params["gauss_fwhm"].value = 5 + app.params["n"].value = 1024 + app.params_layout.setSpacing(0) + + @app.update + def draw(n, lorentz_fwhm, gauss_fwhm, absolute_noise, relative_noise): + x, dx, f, df = setup_axes(n, (-20, 20)) + gaus = gaussian(x, 1, gauss_fwhm, 0) + lore = lorentzian(x, 1, lorentz_fwhm, 0) + ax.set_line_data("Gaussian", x, gaus / gaus.max(), width=2) + ax.set_line_data("Lorentz", x, lore / lore.max(), width=2) + + conv = convolve(gaus, lore) + conv /= conv.max() + conv_noise = add_noise(conv, absolute_noise, relative_noise * 0.1) + gaus_noise = add_noise(gaus, absolute_noise, relative_noise * 0.1) + ax.set_line_data("Voigt", x, conv, width=2) + ax.set_line_data("Noisy Voigt", x, conv_noise, width=2) + ax.set_line_data("Noisy Gaussian", x, gaus_noise, width=2) + + lin_fit = linear_direct_voigt_fit(x, conv_noise) + log_fit = log_direct_voigt_fit(x, conv_noise) + transfo = np.fft.fft(np.fft.fftshift(conv_noise)).real + fft_fit = linear_fft_voigt_fit(f, transfo / transfo.max(), x, dx) + for fit in lin_fit, log_fit, fft_fit: + display_fit(x, fit) + + def display_fit(x: np.ndarray, fit: VoigtFitResult): + ax.set_line_data(fit.name, x, fit.curve) + ax.set_line_data( + fit.name + "gauss", + x, + gaussian(x, 1, fit.gaussian_width, 0), + label=f"{fit.name} - Gaussian {fit.gaussian_width:.2f}", + width=1.5, + ) + ax.set_line_data( + fit.name + "loren", + x, + lorentzian(x, 1, fit.lorentzian_width, 0), + label=f"{fit.name} - Lorentzian {fit.lorentzian_width:.2f}", + width=1.5, + ) + + +@main.command(help="Explore fft fitting") +def demo3(): + x, dx, f, df = setup_axes(1001, (-20, 20)) + with PlotApp( + lorentz_fwhm=np.geomspace(0.1, 1), + gauss_fwhm=np.linspace(1, 10), + relative_noise=np.linspace(0, 2), + absolute_noise=[0, *np.geomspace(1e-6, 1, 201)], + ) as app: + app.set_antialiasing(True) + + ax = app[0] + ind = np.argsort(f) + + @app.update + def draw(lorentz_fwhm, gauss_fwhm, absolute_noise, relative_noise): + gaus_clean = gaussian(x, 1, gauss_fwhm, 0) + lore_clean = lorentzian(x, 1, lorentz_fwhm, 0) + voig_clean = convolve(gaus_clean, lore_clean) + + gaus = add_noise(gaus_clean, absolute_noise, relative_noise * 0.1) + lore = add_noise(lore_clean, absolute_noise, relative_noise * 0.1) + voig = add_noise(voig_clean, absolute_noise, relative_noise * 0.1) + + gaus_f = np.fft.fft(np.fft.fftshift(gaus)).real + lore_f = np.fft.fft(np.fft.fftshift(lore)).real + voig_f = np.fft.fft(np.fft.fftshift(voig)).real + + ax.set_line_data("Transformed Gaussian", f[ind], gaus_f[ind]) + ax.set_line_data("Transformed Lorentz", f[ind], lore_f[ind]) + ax.set_line_data("Transformed Voigt", f[ind], voig_f[ind]) + + gaus_ff = fft_gaussian(f, dx, 1, gauss_fwhm) + lore_ff = fft_lorentzian(f, dx, 1, lorentz_fwhm) + voig_ff = fft_voigt(f, dx, gauss_fwhm, lorentz_fwhm) * voig_clean.sum() + + ax.set_line_data("Expected Gaussian", f[ind], gaus_ff[ind], width=1) + ax.set_line_data("Expected Lorentz", f[ind], lore_ff[ind], width=1) + ax.set_line_data("Expected Voigt", f[ind], voig_ff[ind], width=1) + + +@click.argument("file", type=click.Path(True, dir_okay=False)) +@main.command(name="fit") +def fit_measurement(file): + wl, intens = lab.osa.load_osa_spectrum(file) + wl *= 1e9 + gauss_fit = single_func_fit(wl, intens, elevated_gaussian, extra_params=1) + wl0 = gauss_fit.x0 + baseline = gauss_fit.params[-1] + intens -= baseline + gauss_fit.curve -= baseline + intens[intens < 0] = 0 + wl -= wl0 + with PlotApp(n=np.arange(len(wl) // 2)[1:]) as app: + app.set_antialiasing(True) + + ax1 = app["Measurement"] + ax2 = app["Fourier domain"] + ax1.set_line_data("Measured spectrum", wl, intens, width=2, color="red") + ax1.set_line_data("Initial Guassian fit", wl, gauss_fit.curve, width=2, color="blue") + ax1.plot_widget.plotItem.setLogMode(y=True) + + @app.update + def draw(n): + s = slice(n, -n) if n else slice(None) + + wl_fit = wl[s] + intens_fit = intens[s] - intens[s].min() + f, ind, dwl = freq_axis(wl_fit) + + ax1.set_line_data("Fitted spectrum", wl_fit, intens_fit) + + trans = np.abs(np.fft.fft(intens_fit))[ind] + m = trans.max() + ax2.set_line_data("Transformed", f, trans) + + fft_fit = linear_fft_voigt_fit(f, trans / m, wl_fit, dwl) + lin_fit = linear_direct_voigt_fit(wl_fit, intens_fit) + log_fit = log_direct_voigt_fit(wl_fit, intens_fit) + for fit in fft_fit, lin_fit, log_fit: + display_fit(wl_fit, dwl, f, fit, m) + + def display_fit(x: np.ndarray, dx: float, f: np.ndarray, fit: VoigtFitResult, amp: float): + ax2.set_line_data( + fit.name, f, amp * fft_voigt(f, dx, fit.gaussian_width, fit.lorentzian_width) + ) + ax1.set_line_data(fit.name, x, fit.curve) + ax1.set_line_data( + fit.name + "gauss", + x, + gaussian(x, 1, fit.gaussian_width, 0), + label=f"{fit.name} - Gaussian {fit.gaussian_width:.2f}", + width=1.5, + ) + ax1.set_line_data( + fit.name + "lorentz", + x, + lorentzian(x, 1, fit.lorentzian_width, 0), + label=f"{fit.name} - Lorentzian {fit.lorentzian_width:.2f}", + width=1.5, + ) + + +if __name__ == "__main__": + main() diff --git a/test_scaling.py b/test_scaling.py new file mode 100644 index 0000000..b0cc2c7 --- /dev/null +++ b/test_scaling.py @@ -0,0 +1,46 @@ +"""test what the scaling factore of numpy's fft is""" +import numpy as np +from customfunc.app import PlotApp + + +def setup_axes(n: int, lims=(-10, 10)) -> tuple[np.ndarray, float, np.ndarray, float]: + x, dx = np.linspace(*lims, n, retstep=True) + f = np.fft.fftfreq(n, dx) + return x, dx, f, f[1] - f[0] + + +def abs2(z): + return z.real ** 2 + z.imag ** 2 + + +def g(x, w): + return np.exp(-w * x ** 2) + + +def main(): + + with PlotApp(exp_param=np.linspace(1, 50), n=np.arange(32, 1024)) as app: + + app.set_antialiasing(True) + fax = app["Frequency"] + + @app.update + def draw(exp_param, n): + x, dx, f, df = setup_axes(n) + ind = np.argsort(f) + time_amp = 1 + + freq_amp = np.sqrt(np.pi / exp_param) * time_amp / dx + freq_w = np.pi ** 2 / exp_param + + signal = g(x, exp_param) * time_amp + + transformed = np.fft.fft(np.fft.fftshift(signal)) + fax.set_line_data("transformed", f[ind], transformed[ind].real) + + modeled = g(f, freq_w) * freq_amp + fax.set_line_data("expected", f[ind], modeled[ind]) + + +if __name__ == "__main__": + main()