diff --git a/.drone.yml b/.drone.yml new file mode 100644 index 0000000..d8ce678 --- /dev/null +++ b/.drone.yml @@ -0,0 +1,110 @@ +kind: pipeline +type: docker +name: archlinux + +clone: + depth: 50 + +steps: + - name: archlinux_build + image: archlinux_build:latest + pull: if-not-exists + volumes: + - name: archlinux_ccache + path: /root/.ccache + commands: + # The following command is not required, we included this in the docker + # image of archlinux_build + # - pacman -S --noconfirm ccache openblas fftw pulseaudio pybind11 + - git submodule update --init --recursive + - cmake . + # More than two makes ascee2 irresponsive for now + - make -j2 + + - name: archlinux_test + image: archlinux_build:latest + pull: if-not-exists + commands: + # The following command is not required, we included this in the docker + # image of archlinux_build + # - pacman -S --noconfirm openblas python-pytest fftw pulseaudio python-pip python-scipy python-h5py + - pip install -r requirements.txt + - pip install . + - pytest + + # - name: release-arch + # commands: + # - + +volumes: + - name: archlinux_ccache + host: + path: /tmp/archlinux_ccache + +--- +kind: pipeline +type: docker +name: ubuntu + +clone: + depth: 3 + +volumes: +- name: archlinux_ccache + path: /root/.ccache + +steps: + - name: ubuntu_build + image: ubuntu_build:latest + pull: if-not-exists + volumes: + - name: ubuntu_ccache + path: /root/.ccache + commands: + # The following commands are not required, we included this in the docker + # image of ubuntu_build + #- apt update + #- apt install -y git cmake python3-pybind11 libopenblas-dev python3-pip python3-scipy libusb-1.0-0-dev libpulse-dev python3-h5py fftw-dev + - git submodule update --init --recursive + - cmake . + # More than two makes ascee2 irresponsive for now + - make -j2 + + - name: ubuntu_test + image: ubuntu_build:latest + pull: if-not-exists + commands: + # The following commands are not required, we included this in the docker + # image of ubuntu_build + #- apt update + #- apt install -y python3-pytest fftw pulseaudio python3-pip python3-scipy python3-h5py + - pip install -r requirements.txt + - pip install . + - pytest-3 + +volumes: + - name: ubuntu_ccache + host: + path: /tmp/ubuntu_ccache + +--- +kind: pipeline +type: docker +name: documentation_build + +clone: + depth: 3 + + +steps: + - name: build_docker_master + image: plugins/docker + settings: + repo: ascee/lasp_ascee_nl + tags: latest + username: + from_secret: docker_username + password: + from_secret: docker_password + when: + branch: master diff --git a/.gitignore b/.gitignore index c6e64f4..6caaa0c 100644 --- a/.gitignore +++ b/.gitignore @@ -20,3 +20,4 @@ doc .spyproject .cache _skbuild +acme_log.log diff --git a/Dockerfile b/Dockerfile new file mode 100644 index 0000000..2dcced5 --- /dev/null +++ b/Dockerfile @@ -0,0 +1,16 @@ +FROM archlinux +MAINTAINER J.A. de Jong - j.a.dejong@ascee.nl +RUN pacman --noconfirm -Sy +RUN pacman --noconfirm -S git doxygen graphviz lighttpd python-pip +RUN pip install doxypypy +WORKDIR /root +RUN git clone https://code.ascee.nl/ascee/lasp +WORKDIR /root/lasp +RUN doxygen +RUN rm -rf /srv//http +RUN mv doc/html /srv/http +CMD /usr/bin/lighttpd -D -f /etc/lighttpd/lighttpd.conf + +EXPOSE 80 + + diff --git a/Doxyfile b/Doxyfile index bcd6a6f..b9f844f 100644 --- a/Doxyfile +++ b/Doxyfile @@ -1,4 +1,4 @@ -# Doxyfile 1.8.17 +# Doxyfile 1.9.1 # This file describes the settings to be used by the documentation system # doxygen (www.doxygen.org) for a project. @@ -51,7 +51,7 @@ PROJECT_BRIEF = "Library for Acoustic Signal Processing" # pixels and the maximum width should not exceed 200 pixels. Doxygen will copy # the logo to the output directory. -PROJECT_LOGO = /home/anne/wip/mycode/lasp/img/LASP_200px.png +PROJECT_LOGO = img/LASP_200px.png # The OUTPUT_DIRECTORY tag is used to specify the (relative or absolute) path # into which the generated documentation will be written. If a relative path is @@ -93,6 +93,14 @@ ALLOW_UNICODE_NAMES = NO OUTPUT_LANGUAGE = English +# The OUTPUT_TEXT_DIRECTION tag is used to specify the direction in which all +# documentation generated by doxygen is written. Doxygen will use this +# information to generate all generated output in the proper direction. +# Possible values are: None, LTR, RTL and Context. +# The default value is: None. + +OUTPUT_TEXT_DIRECTION = None + # If the BRIEF_MEMBER_DESC tag is set to YES, doxygen will include brief member # descriptions after the members that are listed in the file and class # documentation (similar to Javadoc). Set to NO to disable this. @@ -225,7 +233,7 @@ MULTILINE_CPP_IS_BRIEF = NO # documentation blocks is shown as doxygen documentation. # The default value is: YES. -PYTHON_DOCSTRING = YES +PYTHON_DOCSTRING = NO # If the INHERIT_DOCS tag is set to YES then an undocumented member inherits the # documentation from any documented member that it re-implements. @@ -250,16 +258,16 @@ TAB_SIZE = 4 # the documentation. An alias has the form: # name=value # For example adding -# "sideeffect=@par Side Effects:^^" +# "sideeffect=@par Side Effects:\n" # will allow you to put the command \sideeffect (or @sideeffect) in the # documentation, which will result in a user-defined paragraph with heading -# "Side Effects:". Note that you cannot put \n's in the value part of an alias -# to insert newlines (in the resulting output). You can put ^^ in the value part -# of an alias to insert a newline as if a physical newline was in the original -# file. When you need a literal { or } or , in the value part of an alias you -# have to escape them by means of a backslash (\), this can lead to conflicts -# with the commands \{ and \} for these it is advised to use the version @{ and -# @} or use a double escape (\\{ and \\}) +# "Side Effects:". You can put \n's in the value part of an alias to insert +# newlines (in the resulting output). You can put ^^ in the value part of an +# alias to insert a newline as if a physical newline was in the original file. +# When you need a literal { or } or , in the value part of an alias you have to +# escape them by means of a backslash (\), this can lead to conflicts with the +# commands \{ and \} for these it is advised to use the version @{ and @} or use +# a double escape (\\{ and \\}) ALIASES = @@ -304,8 +312,8 @@ OPTIMIZE_OUTPUT_SLICE = NO # extension. Doxygen has a built-in mapping, but you can override or extend it # using this tag. The format is ext=language, where ext is a file extension, and # language is one of the parsers supported by doxygen: IDL, Java, JavaScript, -# Csharp (C#), C, C++, Lex, D, PHP, md (Markdown), Objective-C, Python, Slice, -# VHDL, Fortran (fixed format Fortran: FortranFixed, free formatted Fortran: +# Csharp (C#), C, C++, D, PHP, md (Markdown), Objective-C, Python, Slice, VHDL, +# Fortran (fixed format Fortran: FortranFixed, free formatted Fortran: # FortranFree, unknown formatted Fortran: Fortran. In the later case the parser # tries to guess whether the code is fixed or free formatted code, this is the # default for Fortran type files). For instance to make doxygen treat .inc files @@ -458,7 +466,7 @@ LOOKUP_CACHE_SIZE = 0 # than 0 to get more control over the balance between CPU load and processing # speed. At this moment only the input processing can be done using multiple # threads. Since this is still an experimental feature the default is set to 1, -# which effectively disables parallel processing. Please report any issues you +# which efficively disables parallel processing. Please report any issues you # encounter. Generating dot graphs in parallel is controlled by the # DOT_NUM_THREADS setting. # Minimum value: 0, maximum value: 32, default value: 1. @@ -602,12 +610,6 @@ HIDE_SCOPE_NAMES = NO HIDE_COMPOUND_REFERENCE= NO -# If the SHOW_HEADERFILE tag is set to YES then the documentation for a class -# will show which file needs to be included to use the class. -# The default value is: YES. - -SHOW_HEADERFILE = YES - # If the SHOW_INCLUDE_FILES tag is set to YES then doxygen will put a list of # the files that are included by a file in the documentation of that file. # The default value is: YES. @@ -765,8 +767,7 @@ FILE_VERSION_FILTER = # output files in an output format independent way. To create the layout file # that represents doxygen's defaults, run doxygen with the -l option. You can # optionally specify a file name after the option, if omitted DoxygenLayout.xml -# will be used as the name of the layout file. See also section "Changing the -# layout of pages" for information. +# will be used as the name of the layout file. # # Note that if you run doxygen from a directory containing a file called # DoxygenLayout.xml, doxygen will parse it automatically even if the LAYOUT_FILE @@ -812,26 +813,18 @@ WARNINGS = YES WARN_IF_UNDOCUMENTED = YES # If the WARN_IF_DOC_ERROR tag is set to YES, doxygen will generate warnings for -# potential errors in the documentation, such as documenting some parameters in -# a documented function twice, or documenting parameters that don't exist or -# using markup commands wrongly. +# potential errors in the documentation, such as not documenting some parameters +# in a documented function, or documenting parameters that don't exist or using +# markup commands wrongly. # The default value is: YES. WARN_IF_DOC_ERROR = YES -# If WARN_IF_INCOMPLETE_DOC is set to YES, doxygen will warn about incomplete -# function parameter documentation. If set to NO, doxygen will accept that some -# parameters have no documentation without warning. -# The default value is: YES. - -WARN_IF_INCOMPLETE_DOC = YES - # This WARN_NO_PARAMDOC option can be enabled to get warnings for functions that # are documented, but have no documentation for their parameters or return -# value. If set to NO, doxygen will only warn about wrong parameter -# documentation, but not about the absence of documentation. If EXTRACT_ALL is -# set to YES then this flag will automatically be disabled. See also -# WARN_IF_INCOMPLETE_DOC +# value. If set to NO, doxygen will only warn about wrong or incomplete +# parameter documentation, but not about the absence of documentation. If +# EXTRACT_ALL is set to YES then this flag will automatically be disabled. # The default value is: NO. WARN_NO_PARAMDOC = NO @@ -857,10 +850,7 @@ WARN_FORMAT = "$file:$line: $text" # The WARN_LOGFILE tag can be used to specify a file to which warning and error # messages should be written. If left blank the output is written to standard -# error (stderr). In case the file specified cannot be opened for writing the -# warning and error messages are written to standard error. When as file - is -# specified the warning and error messages are written to standard output -# (stdout). +# error (stderr). WARN_LOGFILE = @@ -898,55 +888,16 @@ INPUT_ENCODING = UTF-8 # # If left blank the following patterns are tested:*.c, *.cc, *.cxx, *.cpp, # *.c++, *.java, *.ii, *.ixx, *.ipp, *.i++, *.inl, *.idl, *.ddl, *.odl, *.h, -# *.hh, *.hxx, *.hpp, *.h++, *.l, *.cs, *.d, *.php, *.php4, *.php5, *.phtml, -# *.inc, *.m, *.markdown, *.md, *.mm, *.dox (to be provided as doxygen C -# comment), *.py, *.pyw, *.f90, *.f95, *.f03, *.f08, *.f18, *.f, *.for, *.vhd, -# *.vhdl, *.ucf, *.qsf and *.ice. +# *.hh, *.hxx, *.hpp, *.h++, *.cs, *.d, *.php, *.php4, *.php5, *.phtml, *.inc, +# *.m, *.markdown, *.md, *.mm, *.dox (to be provided as doxygen C comment), +# *.py, *.pyw, *.f90, *.f95, *.f03, *.f08, *.f18, *.f, *.for, *.vhd, *.vhdl, +# *.ucf, *.qsf and *.ice. FILE_PATTERNS = *.c \ - *.cc \ - *.cxx \ *.cpp \ - *.c++ \ - *.java \ - *.ii \ - *.ixx \ - *.ipp \ - *.i++ \ - *.inl \ - *.idl \ - *.ddl \ - *.odl \ *.h \ - *.hh \ - *.hxx \ - *.hpp \ - *.h++ \ - *.cs \ - *.d \ - *.php \ - *.php4 \ - *.php5 \ - *.phtml \ - *.inc \ - *.m \ - *.markdown \ *.md \ - *.mm \ - *.dox \ - *.py \ - *.pyw \ - *.f90 \ - *.f95 \ - *.f03 \ - *.f08 \ - *.f \ - *.for \ - *.tcl \ - *.vhd \ - *.vhdl \ - *.ucf \ - *.qsf + *.py=scripts/py_filter # The RECURSIVE tag can be used to specify whether or not subdirectories should # be searched for input files as well. @@ -983,7 +934,7 @@ EXCLUDE_PATTERNS = # (namespaces, classes, functions, etc.) that should be excluded from the # output. The symbol name can be a fully qualified name, a word, or if the # wildcard * is used, a substring. Examples: ANamespace, AClass, -# ANamespace::AClass, ANamespace::*Test +# AClass::ANamespace, ANamespace::*Test # # Note that the wildcards are matched against the file with absolute path, so to # exclude all test directories use the pattern */test/* @@ -1158,6 +1109,44 @@ USE_HTAGS = NO VERBATIM_HEADERS = YES +# If the CLANG_ASSISTED_PARSING tag is set to YES then doxygen will use the +# clang parser (see: +# http://clang.llvm.org/) for more accurate parsing at the cost of reduced +# performance. This can be particularly helpful with template rich C++ code for +# which doxygen's built-in parser lacks the necessary type information. +# Note: The availability of this option depends on whether or not doxygen was +# generated with the -Duse_libclang=ON option for CMake. +# The default value is: NO. + +CLANG_ASSISTED_PARSING = NO + +# If clang assisted parsing is enabled and the CLANG_ADD_INC_PATHS tag is set to +# YES then doxygen will add the directory of each input to the include path. +# The default value is: YES. + +CLANG_ADD_INC_PATHS = YES + +# If clang assisted parsing is enabled you can provide the compiler with command +# line options that you would normally use when invoking the compiler. Note that +# the include paths will already be set by doxygen for the files and directories +# specified with INPUT and INCLUDE_PATH. +# This tag requires that the tag CLANG_ASSISTED_PARSING is set to YES. + +CLANG_OPTIONS = + +# If clang assisted parsing is enabled you can provide the clang parser with the +# path to the directory containing a file called compile_commands.json. This +# file is the compilation database (see: +# http://clang.llvm.org/docs/HowToSetupToolingForLLVM.html) containing the +# options used when the source files were built. This is equivalent to +# specifying the -p option to a clang tool, such as clang-check. These options +# will then be passed to the parser. Any options specified with CLANG_OPTIONS +# will be added as well. +# Note: The availability of this option depends on whether or not doxygen was +# generated with the -Duse_libclang=ON option for CMake. + +CLANG_DATABASE_PATH = + #--------------------------------------------------------------------------- # Configuration options related to the alphabetical class index #--------------------------------------------------------------------------- @@ -1268,7 +1257,7 @@ HTML_EXTRA_FILES = # The HTML_COLORSTYLE_HUE tag controls the color of the HTML output. Doxygen # will adjust the colors in the style sheet and background images according to -# this color. Hue is specified as an angle on a color-wheel, see +# this color. Hue is specified as an angle on a colorwheel, see # https://en.wikipedia.org/wiki/Hue for more information. For instance the value # 0 represents red, 60 is yellow, 120 is green, 180 is cyan, 240 is blue, 300 # purple, and 360 is red again. @@ -1278,7 +1267,7 @@ HTML_EXTRA_FILES = HTML_COLORSTYLE_HUE = 220 # The HTML_COLORSTYLE_SAT tag controls the purity (or saturation) of the colors -# in the HTML output. For a value of 0 the output will use gray-scales only. A +# in the HTML output. For a value of 0 the output will use grayscales only. A # value of 255 will produce the most vivid colors. # Minimum value: 0, maximum value: 255, default value: 100. # This tag requires that the tag GENERATE_HTML is set to YES. @@ -1360,13 +1349,6 @@ GENERATE_DOCSET = NO DOCSET_FEEDNAME = "Doxygen generated docs" -# This tag determines the URL of the docset feed. A documentation feed provides -# an umbrella under which multiple documentation sets from a single provider -# (such as a company or product suite) can be grouped. -# This tag requires that the tag GENERATE_DOCSET is set to YES. - -DOCSET_FEEDURL = - # This tag specifies a string that should uniquely identify the documentation # set bundle. This should be a reverse domain-name style string, e.g. # com.mycompany.MyDocSet. Doxygen will append .docset to the name. @@ -1392,12 +1374,8 @@ DOCSET_PUBLISHER_NAME = Publisher # If the GENERATE_HTMLHELP tag is set to YES then doxygen generates three # additional HTML index files: index.hhp, index.hhc, and index.hhk. The # index.hhp is a project file that can be read by Microsoft's HTML Help Workshop -# on Windows. In the beginning of 2021 Microsoft took the original page, with -# a.o. the download links, offline the HTML help workshop was already many years -# in maintenance mode). You can download the HTML help workshop from the web -# archives at Installation executable (see: -# http://web.archive.org/web/20160201063255/http://download.microsoft.com/downlo -# ad/0/A/9/0A939EF6-E31C-430F-A3DF-DFAE7960D564/htmlhelp.exe). +# (see: +# https://www.microsoft.com/en-us/download/details.aspx?id=21138) on Windows. # # The HTML Help Workshop contains a compiler that can convert all HTML output # generated by doxygen into a single compiled HTML file (.chm). Compiled HTML @@ -1556,28 +1534,16 @@ DISABLE_INDEX = NO # to work a browser that supports JavaScript, DHTML, CSS and frames is required # (i.e. any modern browser). Windows users are probably better off using the # HTML help feature. Via custom style sheets (see HTML_EXTRA_STYLESHEET) one can -# further fine tune the look of the index (see "Fine-tuning the output"). As an -# example, the default style sheet generated by doxygen has an example that -# shows how to put an image at the root of the tree instead of the PROJECT_NAME. -# Since the tree basically has the same information as the tab index, you could -# consider setting DISABLE_INDEX to YES when enabling this option. +# further fine-tune the look of the index. As an example, the default style +# sheet generated by doxygen has an example that shows how to put an image at +# the root of the tree instead of the PROJECT_NAME. Since the tree basically has +# the same information as the tab index, you could consider setting +# DISABLE_INDEX to YES when enabling this option. # The default value is: NO. # This tag requires that the tag GENERATE_HTML is set to YES. GENERATE_TREEVIEW = YES -# When both GENERATE_TREEVIEW and DISABLE_INDEX are set to YES, then the -# FULL_SIDEBAR option determines if the side bar is limited to only the treeview -# area (value NO) or if it should extend to the full height of the window (value -# YES). Setting this to YES gives a layout similar to -# https://docs.readthedocs.io with more room for contents, but less room for the -# project logo, title, and description. If either GENERATE_TREEVIEW or -# DISABLE_INDEX is set to NO, this option has no effect. -# The default value is: NO. -# This tag requires that the tag GENERATE_HTML is set to YES. - -FULL_SIDEBAR = NO - # The ENUM_VALUES_PER_LINE tag can be used to set the number of enum values that # doxygen will group on one line in the generated HTML documentation. # @@ -1602,13 +1568,6 @@ TREEVIEW_WIDTH = 250 EXT_LINKS_IN_WINDOW = NO -# If the OBFUSCATE_EMAILS tag is set to YES, doxygen will obfuscate email -# addresses. -# The default value is: YES. -# This tag requires that the tag GENERATE_HTML is set to YES. - -OBFUSCATE_EMAILS = YES - # If the HTML_FORMULA_FORMAT option is set to svg, doxygen will use the pdf2svg # tool (see https://github.com/dawbarton/pdf2svg) or inkscape (see # https://inkscape.org) to generate formulas as SVG images instead of PNGs for @@ -1657,29 +1616,11 @@ FORMULA_MACROFILE = USE_MATHJAX = NO -# With MATHJAX_VERSION it is possible to specify the MathJax version to be used. -# Note that the different versions of MathJax have different requirements with -# regards to the different settings, so it is possible that also other MathJax -# settings have to be changed when switching between the different MathJax -# versions. -# Possible values are: MathJax_2 and MathJax_3. -# The default value is: MathJax_2. -# This tag requires that the tag USE_MATHJAX is set to YES. - -MATHJAX_VERSION = MathJax_2 - # When MathJax is enabled you can set the default output format to be used for -# the MathJax output. For more details about the output format see MathJax -# version 2 (see: -# http://docs.mathjax.org/en/v2.7-latest/output.html) and MathJax version 3 -# (see: -# http://docs.mathjax.org/en/latest/web/components/output.html). +# the MathJax output. See the MathJax site (see: +# http://docs.mathjax.org/en/v2.7-latest/output.html) for more details. # Possible values are: HTML-CSS (which is slower, but has the best -# compatibility. This is the name for Mathjax version 2, for MathJax version 3 -# this will be translated into chtml), NativeMML (i.e. MathML. Only supported -# for NathJax 2. For MathJax version 3 chtml will be used instead.), chtml (This -# is the name for Mathjax version 3, for MathJax version 2 this will be -# translated into HTML-CSS) and SVG. +# compatibility), NativeMML (i.e. MathML) and SVG. # The default value is: HTML-CSS. # This tag requires that the tag USE_MATHJAX is set to YES. @@ -1692,21 +1633,15 @@ MATHJAX_FORMAT = HTML-CSS # MATHJAX_RELPATH should be ../mathjax. The default value points to the MathJax # Content Delivery Network so you can quickly see the result without installing # MathJax. However, it is strongly recommended to install a local copy of -# MathJax from https://www.mathjax.org before deployment. The default value is: -# - in case of MathJax version 2: https://cdn.jsdelivr.net/npm/mathjax@2 -# - in case of MathJax version 3: https://cdn.jsdelivr.net/npm/mathjax@3 +# MathJax from https://www.mathjax.org before deployment. +# The default value is: https://cdn.jsdelivr.net/npm/mathjax@2. # This tag requires that the tag USE_MATHJAX is set to YES. MATHJAX_RELPATH = https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.2/ # The MATHJAX_EXTENSIONS tag can be used to specify one or more MathJax # extension names that should be enabled during MathJax rendering. For example -# for MathJax version 2 (see -# https://docs.mathjax.org/en/v2.7-latest/tex.html#tex-and-latex-extensions): # MATHJAX_EXTENSIONS = TeX/AMSmath TeX/AMSsymbols -# For example for MathJax version 3 (see -# http://docs.mathjax.org/en/latest/input/tex/extensions/index.html): -# MATHJAX_EXTENSIONS = ams # This tag requires that the tag USE_MATHJAX is set to YES. MATHJAX_EXTENSIONS = @@ -1886,31 +1821,29 @@ PAPER_TYPE = a4 EXTRA_PACKAGES = -# The LATEX_HEADER tag can be used to specify a user-defined LaTeX header for -# the generated LaTeX document. The header should contain everything until the -# first chapter. If it is left blank doxygen will generate a standard header. It -# is highly recommended to start with a default header using -# doxygen -w latex new_header.tex new_footer.tex new_stylesheet.sty -# and then modify the file new_header.tex. See also section "Doxygen usage" for -# information on how to generate the default header that doxygen normally uses. +# The LATEX_HEADER tag can be used to specify a personal LaTeX header for the +# generated LaTeX document. The header should contain everything until the first +# chapter. If it is left blank doxygen will generate a standard header. See +# section "Doxygen usage" for information on how to let doxygen write the +# default header to a separate file. # -# Note: Only use a user-defined header if you know what you are doing! -# Note: The header is subject to change so you typically have to regenerate the -# default header when upgrading to a newer version of doxygen. The following -# commands have a special meaning inside the header (and footer): For a -# description of the possible markers and block names see the documentation. +# Note: Only use a user-defined header if you know what you are doing! The +# following commands have a special meaning inside the header: $title, +# $datetime, $date, $doxygenversion, $projectname, $projectnumber, +# $projectbrief, $projectlogo. Doxygen will replace $title with the empty +# string, for the replacement values of the other commands the user is referred +# to HTML_HEADER. # This tag requires that the tag GENERATE_LATEX is set to YES. LATEX_HEADER = -# The LATEX_FOOTER tag can be used to specify a user-defined LaTeX footer for -# the generated LaTeX document. The footer should contain everything after the -# last chapter. If it is left blank doxygen will generate a standard footer. See +# The LATEX_FOOTER tag can be used to specify a personal LaTeX footer for the +# generated LaTeX document. The footer should contain everything after the last +# chapter. If it is left blank doxygen will generate a standard footer. See # LATEX_HEADER for more information on how to generate a default footer and what -# special commands can be used inside the footer. See also section "Doxygen -# usage" for information on how to generate the default footer that doxygen -# normally uses. Note: Only use a user-defined footer if you know what you are -# doing! +# special commands can be used inside the footer. +# +# Note: Only use a user-defined footer if you know what you are doing! # This tag requires that the tag GENERATE_LATEX is set to YES. LATEX_FOOTER = @@ -1955,7 +1888,8 @@ USE_PDFLATEX = YES # If the LATEX_BATCHMODE tag is set to YES, doxygen will add the \batchmode # command to the generated LaTeX files. This will instruct LaTeX to keep running -# if errors occur, instead of asking the user for help. +# if errors occur, instead of asking the user for help. This option is also used +# when generating formulas in HTML. # The default value is: NO. # This tag requires that the tag GENERATE_LATEX is set to YES. @@ -1968,6 +1902,16 @@ LATEX_BATCHMODE = NO LATEX_HIDE_INDICES = NO +# If the LATEX_SOURCE_CODE tag is set to YES then doxygen will include source +# code with syntax highlighting in the LaTeX output. +# +# Note that which sources are shown also depends on other settings such as +# SOURCE_BROWSER. +# The default value is: NO. +# This tag requires that the tag GENERATE_LATEX is set to YES. + +LATEX_SOURCE_CODE = NO + # The LATEX_BIB_STYLE tag can be used to specify the style to use for the # bibliography, e.g. plainnat, or ieeetr. See # https://en.wikipedia.org/wiki/BibTeX and \cite for more info. @@ -2048,6 +1992,16 @@ RTF_STYLESHEET_FILE = RTF_EXTENSIONS_FILE = +# If the RTF_SOURCE_CODE tag is set to YES then doxygen will include source code +# with syntax highlighting in the RTF output. +# +# Note that which sources are shown also depends on other settings such as +# SOURCE_BROWSER. +# The default value is: NO. +# This tag requires that the tag GENERATE_RTF is set to YES. + +RTF_SOURCE_CODE = NO + #--------------------------------------------------------------------------- # Configuration options related to the man page output #--------------------------------------------------------------------------- @@ -2144,6 +2098,15 @@ GENERATE_DOCBOOK = NO DOCBOOK_OUTPUT = docbook +# If the DOCBOOK_PROGRAMLISTING tag is set to YES, doxygen will include the +# program listings (including syntax highlighting and cross-referencing +# information) to the DOCBOOK output. Note that enabling this will significantly +# increase the size of the DOCBOOK output. +# The default value is: NO. +# This tag requires that the tag GENERATE_DOCBOOK is set to YES. + +DOCBOOK_PROGRAMLISTING = NO + #--------------------------------------------------------------------------- # Configuration options for the AutoGen Definitions output #--------------------------------------------------------------------------- @@ -2156,6 +2119,10 @@ DOCBOOK_OUTPUT = docbook GENERATE_AUTOGEN_DEF = NO +#--------------------------------------------------------------------------- +# Configuration options related to Sqlite3 output +#--------------------------------------------------------------------------- + #--------------------------------------------------------------------------- # Configuration options related to the Perl module output #--------------------------------------------------------------------------- @@ -2322,6 +2289,15 @@ EXTERNAL_PAGES = YES # Configuration options related to the dot tool #--------------------------------------------------------------------------- +# If the CLASS_DIAGRAMS tag is set to YES, doxygen will generate a class diagram +# (in HTML and LaTeX) for classes with base or super classes. Setting the tag to +# NO turns the diagrams off. Note that this option also works with HAVE_DOT +# disabled, but it is recommended to install and use dot, since it yields more +# powerful graphs. +# The default value is: YES. + +CLASS_DIAGRAMS = YES + # You can include diagrams made with dia in doxygen documentation. Doxygen will # then run dia to produce the diagram and insert it in the documentation. The # DIA_PATH tag allows you to specify the directory where the dia binary resides. @@ -2340,7 +2316,7 @@ HIDE_UNDOC_RELATIONS = YES # http://www.graphviz.org/), a graph visualization toolkit from AT&T and Lucent # Bell Labs. The other options in this section have no effect if this option is # set to NO -# The default value is: NO. +# The default value is: YES. HAVE_DOT = YES @@ -2378,14 +2354,11 @@ DOT_FONTSIZE = 10 DOT_FONTPATH = -# If the CLASS_GRAPH tag is set to YES (or GRAPH) then doxygen will generate a -# graph for each documented class showing the direct and indirect inheritance -# relations. In case HAVE_DOT is set as well dot will be used to draw the graph, -# otherwise the built-in generator will be used. If the CLASS_GRAPH tag is set -# to TEXT the direct and indirect inheritance relations will be shown as texts / -# links. -# Possible values are: NO, YES, TEXT and GRAPH. +# If the CLASS_GRAPH tag is set to YES then doxygen will generate a graph for +# each documented class showing the direct and indirect inheritance relations. +# Setting this tag to YES will force the CLASS_DIAGRAMS tag to NO. # The default value is: YES. +# This tag requires that the tag HAVE_DOT is set to YES. CLASS_GRAPH = YES @@ -2514,13 +2487,6 @@ GRAPHICAL_HIERARCHY = YES DIRECTORY_GRAPH = YES -# The DIR_GRAPH_MAX_DEPTH tag can be used to limit the maximum number of levels -# of child directories generated in directory dependency graphs by dot. -# Minimum value: 1, maximum value: 25, default value: 1. -# This tag requires that the tag DIRECTORY_GRAPH is set to YES. - -DIR_GRAPH_MAX_DEPTH = 1 - # The DOT_IMAGE_FORMAT tag can be used to set the image format of the images # generated by dot. For an explanation of the image formats see the section # output formats in the documentation of the dot tool (Graphviz (see: @@ -2528,7 +2494,9 @@ DIR_GRAPH_MAX_DEPTH = 1 # Note: If you choose svg you need to set HTML_FILE_EXTENSION to xhtml in order # to make the SVG files visible in IE 9+ (other browsers do not have this # requirement). -# Possible values are: png, jpg, gif, svg, png:gd, png:gd:gd, png:cairo, +# Possible values are: png, png:cairo, png:cairo:cairo, png:cairo:gd, png:gd, +# png:gd:gd, jpg, jpg:cairo, jpg:cairo:gd, jpg:gd, jpg:gd:gd, gif, gif:cairo, +# gif:cairo:gd, gif:gd, gif:gd:gd, svg, png:gd, png:gd:gd, png:cairo, # png:cairo:gd, png:cairo:cairo, png:cairo:gdiplus, png:gdiplus and # png:gdiplus:gdiplus. # The default value is: png. @@ -2574,10 +2542,10 @@ MSCFILE_DIRS = DIAFILE_DIRS = # When using plantuml, the PLANTUML_JAR_PATH tag should be used to specify the -# path where java can find the plantuml.jar file or to the filename of jar file -# to be used. If left blank, it is assumed PlantUML is not used or called during -# a preprocessing step. Doxygen will generate a warning when it encounters a -# \startuml command in this case and will not generate output for the diagram. +# path where java can find the plantuml.jar file. If left blank, it is assumed +# PlantUML is not used or called during a preprocessing step. Doxygen will +# generate a warning when it encounters a \startuml command in this case and +# will not generate output for the diagram. PLANTUML_JAR_PATH = @@ -2639,8 +2607,6 @@ DOT_MULTI_TARGETS = NO # If the GENERATE_LEGEND tag is set to YES doxygen will generate a legend page # explaining the meaning of the various boxes and arrows in the dot generated # graphs. -# Note: This tag requires that UML_LOOK isn't set, i.e. the doxygen internal -# graphical representation for inheritance and collaboration diagrams is used. # The default value is: YES. # This tag requires that the tag HAVE_DOT is set to YES. @@ -2649,8 +2615,8 @@ GENERATE_LEGEND = YES # If the DOT_CLEANUP tag is set to YES, doxygen will remove the intermediate # files that are used to generate the various graphs. # -# Note: This setting is not only used for dot files but also for msc temporary -# files. +# Note: This setting is not only used for dot files but also for msc and +# plantuml temporary files. # The default value is: YES. DOT_CLEANUP = YES diff --git a/README.md b/README.md index 98bfcc0..738d6dd 100644 --- a/README.md +++ b/README.md @@ -1,108 +1,112 @@ # Library for Acoustic Signal Processing -Welcome to LASP: Library for Acoustic Signal Processing. LASP is a C++ library -with a Python interface which is supposed to process (multi-) microphone -acoustic data in real time on a PC and output results. +- Master branch: [![Build Status](https://drone.ascee.nl/api/badges/ASCEE/lasp/status.svg?ref=refs/heads/master)](https://drone.ascee.nl/ASCEE/lasp) +- Develop branch: [![Build Status](https://drone.ascee.nl/api/badges/ASCEE/lasp/status.svg?ref=refs/heads/develop)](https://drone.ascee.nl/ASCEE/lasp) -At the point in time of this writing, we are yet unsure whether the Raspberry -PI will have enough computational power to this end, but may be by the time it -is finished, we have a new faster generation :). + +Welcome to LASP: Library for Acoustic Signal Processing. LASP is a C++ library +with a Python interface which is supposed to acquire and process (multi) sensor data in real time on a PC and output results. Current features that are implemented: -- Compile-time determination of the floating-point accuracy (32/64 bit) -- Fast convolution FIR filter implementation -- Sample rate decimation by an integer factor of 4. -- Octave filterbank FIR filters designed to comply with IEC 61260 - (1995). + +- Communication with data acquisition (DAQ) devices, of which: + - Internal sound cards via the [RtAudio](http://www.music.mcgill.ca/~gary/rtaudio) backend. Many thanks to Gary P. Scavone et al. + - [Measurement Computing](https://www.mccdaq.com) [DT9838A](https://www.mccdaq.com/Products/Sound-Vibration-DAQ/DT9837) signal analyzer. +- Configuration of DAQ devices: AC coupling, IEPE, sensitivity physical + quantities. +- Recording of signals from these DAQ devices, and storing in a HDF5 file. +- Filter designers to create A/C sound pressure weighting +- Biquad filter designers for low pass, high pass, peaking and notch filters +- A Peak Programme Meter (PPM) to monitor signal levels from DAQ and to watch + for signal clipping. +- A signal generator to create sine waves, sweeps and noise (white / pink). +- Equalizers to equalize the output prior to sending. - Averaged power spectra and power spectral density determination using Welch' method. Taper functions of Hann, Hamming, Bartlett and Blackman are provided. -- A thread-safe job queue including routines to create worker threads. -- Several linear algebra routines (wrappers around BLAS and LAPACK). -- A nice debug tracer implementation -- Third octave filter bank FIR filters designed to comply with IEC 61260 +- (One third) octave filter bank filters designed to comply with IEC 61260 (1995). - Slow and fast time updates of (A/C/Z) weighted sound pressure levels +- Full Sound Level Meter implementation +- Real time Sound Level meter, Power / Transfer function estimator +- Spectra data smoothing algorithms +- Sensor calibration for microphones Future features (wish-list) -- Conventional and delay-and-sum beam-forming algorithms -For now, the source code is well-documented but it requires some -additional documentation (the math behind it). This will be published -in a sister repository (https://code.ascee.nl/ascee/lasp-doc). +- Conventional and delay-and-sum beam-forming algorithms +- Impedance tube measurement processing + +For now, the source code is well-documented on [lasp.ascee.nl](https://lasp.ascee.nl) but it requires some +additional documentation (the math behind it). This is maintained +in a sister repository [lasp-doc](https://code.ascee.nl/ascee/lasp-doc). The +most recent If you have any question(s), please feel free to contact us: info@ascee.nl. +# Installation -# Building from source +## Dependencies -Two commands that install all requirements (for Ubuntu / Linux Mint) - -- `pip install scipy numpy build scikit-build appdirs` -- `sudo apt install libusb-dev libpulse-dev libboost-dev` - -## Runtime dependencies (Linux) - -- FFTW (For really fast FFT's). If compiled with Ffftpack, this library is not - required. -- libUlDAQ, for the Measurement Computing DT9837A USB DAQ box - - GNU Autotools, for compiling libUlDAQ -- RtAudio, for Audio DAQ backends -- libusb -- BLAS (OpenBLAS, other). - -## Editable install - -In the root directory of the repository, run: - -- `pip3 isntall --user -e .` -- `cmake .` -- `make -j` +- `$ sudo apt install python3-pybind11 libopenblas-dev python3-pip python3-scipy libusb-1.0-0-dev libpulse-dev cmake-curses-gui python3-h5py` +- `$ pip3 install --user -r requirements.txt` - -## Build dependencies - -Optional dependencies, which can be turned ON/OFF using CMake: - -- Build tools: compiler [http://cmake.org](CMake), the Python packages: - - Scipy - - Numpy - - appdirs - These can all be installed using: - -- The following Python packages need also be available: - - `Scipy` (which includes Numpy). Install with `sudo apt install - python3-scipy`, or `pacman -S scipy`. - - `appdirs`, which can be grabbed from [https://pypi.org](Pypi) - - - -## Compilation of LASP - -### Archlinux - -Compiling the code on Archlinux requires the following packages to be available: - -- openblas-lapack (AUR) -- Python>=3.7 -- Numpy (Python-numpy) -- Cython - -### Ubuntu / Linux Mint - -*Only tested with Linux Mint 18.04*, we require the following packages for -compilation: - -- build-essential -- cython -- python3-numpy -- libopenblas - -If building RtAudio with the ALSA backend: +If building RtAudio with the ALSA backend, you will also require the following packages: - libclalsadrv-dev -- libopenblas-base -- libopenblas-dev -- libusb-1.0-dev +If building RtAudio with the Jack Audio Connection Kit (JACK) backend, you will also require the following packages: + +- libjack-jackd2-dev + +## Download & build + +- `$ git clone --recursive https://code.ascee.nl/ASCEE/lasp.git` +- `$ cd lasp` + +For a release build: + +- `$ cmake .` + +or optionally for a custom build: + +- `$ ccmake .` + +Configure and run: + +- `$ make -j` + +### Build documentation + +In directory: + +`$ sudo apt install doxygen graphviz` +`$ pip install doxypypy` + +While still in lasp dir: + +`$ doxygen` + +This will build the documentation. It can be read by: + +`$ doc/html/index.html` + +Or via docker: + +`$ docker build -t lasp_ascee_nl:latest .` + +## Install + +For an editable install (while developing): + +- `$ pip3 install --prefix=$HOME/.local -e .` + +To install locally, for a fixed version: + +- `$ pip3 install --prefix=$HOME/.local` + +## Usage + +- See examples directories for IPython notebooks. +- Please refer to the [documentation](https://lasp.ascee.nl/) for features. diff --git a/test/test_biquadbank.py b/examples/example_biquadbank.py similarity index 88% rename from test/test_biquadbank.py rename to examples/example_biquadbank.py index be0a1da..1f5315a 100644 --- a/test/test_biquadbank.py +++ b/examples/example_biquadbank.py @@ -2,15 +2,10 @@ import numpy as np from lasp import SeriesBiquad, AvPowerSpectra from lasp.filter import SPLFilterDesigner - import matplotlib.pyplot as plt from scipy.signal import sosfreqz -# plt.close('all') +plt.close('all') -# def test_cppslm2(): -# """ -# Generate a sine wave, now A-weighted -# """ fs = 48000 omg = 2*np.pi*1000 diff --git a/pyproject.toml b/pyproject.toml deleted file mode 100644 index 6ac7893..0000000 --- a/pyproject.toml +++ /dev/null @@ -1,17 +0,0 @@ -[project] # Project metadata -name = "lasp" -readme = "README.md" -requires-python = ">=3.8" -license = { "file" = "LICENSE" } -authors = [{ "name" = "J.A. de Jong et al.", "email" = "info@ascee.nl" }] - -classifiers = [ - "Topic :: Scientific/Engineering", - "Programming Language :: Python :: 3.8", - "Operating System :: POSIX :: Linux", - "Operating System :: Microsoft :: Windows", -] - -# urls = { "Documentation" = "https://" } -dynamic = ["version", "description"] - diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..c7a9274 --- /dev/null +++ b/requirements.txt @@ -0,0 +1,3 @@ +appdirs +dataclasses_json +matplotlib diff --git a/scripts/py_filter b/scripts/py_filter new file mode 100755 index 0000000..2d6b1f0 --- /dev/null +++ b/scripts/py_filter @@ -0,0 +1,4 @@ +#!/bin/bash +# Script used to filter Python files such that it can be used together with +# Doxygen +doxypypy -a -c $1 diff --git a/setup.py b/setup.py index b46819d..66ceecd 100644 --- a/setup.py +++ b/setup.py @@ -1,16 +1,22 @@ -import glob +import glob, os import platform from setuptools import setup if 'Linux' in platform.platform(): - extension = list(glob.glob('src/lasp/lasp_cpp.cpython*')) - if len(extension) == 0: + ext_name_glob = 'lasp_cpp.cpython*' + extensions = list(glob.glob('src/lasp/' + ext_name_glob)) + + # Split of path from file. + ext_names = [os.path.split(a)[1] for a in extensions] + + print(extensions) + if len(extensions) == 0: raise RuntimeError('Please first run CMake to build extension') - elif len(extension) > 1: + elif len(extensions) > 1: raise RuntimeError('Too many extension files found') - pkgdata = extension + pkgdata = ext_names else: raise RuntimeError('Not yet Windows-proof') @@ -24,9 +30,6 @@ classifiers = [ keywords = ["DSP", "DAQ", "Signal processing"] -with open('README.md', 'r') as f: - readme = f.read() - setup( name="lasp", @@ -40,12 +43,11 @@ setup( classifiers=classifiers, keywords=keywords, license="MIT", - readme=readme, dependencies=["numpy", "scipy", "appdirs", "h5py", "appdirs", - "dataclasses_json"], + "dataclasses_json"], + package_dir={"": "src"}, packages=['lasp', 'lasp.filter', 'lasp.tools'], - data_files = pkgdata, include_package_data=True, - package_dir={'': 'src'}, + package_data={'lasp': pkgdata}, python_requires='>=3.8', ) diff --git a/src/lasp/__init__.py b/src/lasp/__init__.py index d48bd0d..b6f3e90 100644 --- a/src/lasp/__init__.py +++ b/src/lasp/__init__.py @@ -1,4 +1,3 @@ -# Comments are what is imported, state of 6-8-2021 """ LASP: Library for Acoustic Signal Processing @@ -14,6 +13,7 @@ from .lasp_octavefilter import * from .lasp_slm import * # SLM, Dummy from .lasp_record import * # RecordStatus, Recording from .lasp_daqconfigs import * +from .lasp_measurementset import * # from .lasp_siggen import * # SignalType, NoiseType, SiggenMessage, SiggenData, Siggen # from .lasp_weighcal import * # WeighCal # from .tools import * # SmoothingType, smoothSpectralData, SmoothingWidth diff --git a/src/lasp/device/CMakeLists.txt b/src/lasp/device/CMakeLists.txt index 052d997..e576d23 100644 --- a/src/lasp/device/CMakeLists.txt +++ b/src/lasp/device/CMakeLists.txt @@ -8,6 +8,7 @@ add_library(lasp_device_lib OBJECT lasp_rtaudiodaq.cpp lasp_streammgr.cpp lasp_uldaq.cpp + lasp_uldaq_impl.cpp ) # Callback requires certain arguments that are not used by code. This disables diff --git a/src/lasp/device/lasp_daq.h b/src/lasp/device/lasp_daq.h index a921a09..dbc7392 100644 --- a/src/lasp/device/lasp_daq.h +++ b/src/lasp/device/lasp_daq.h @@ -44,7 +44,6 @@ public: systemError, threadError, logicError, - apiSpecificError }; /** @@ -59,6 +58,7 @@ public: {StreamError::threadError, "Thread error"}, {StreamError::logicError, "Logic error (probably a bug)"}, }; + bool isRunning = false; /** * @brief Check if stream has error @@ -78,6 +78,25 @@ public: * @return as described above. */ bool runningOK() const { return isRunning && !error(); } + + }; // End of class StreamStatus + + using rte = std::runtime_error; + /** + * @brief Used for internal throwing of exceptions. + */ + class StreamException : public rte { + using StreamError = StreamStatus::StreamError; + + public: + StreamStatus::StreamError e; + StreamException(const StreamStatus::StreamError e) + : rte(StreamStatus::errorMessages.at(e)), e(e) {} + StreamException(const StreamStatus::StreamError e, + const std::string &additional_info) + : rte(StreamStatus::errorMessages.at(e) + ": " + additional_info), + e(e) {} + operator StreamError() { return e; } }; /** @@ -161,7 +180,7 @@ public: * * * @return Maximum offset from 0 before clipping. */ - dvec inputRangeForEnabledChannels(const bool include_monitor=true) const; + dvec inputRangeForEnabledChannels(const bool include_monitor = true) const; /** * @brief Returns datatype (enum) corresponding to the datatype of the @@ -176,7 +195,7 @@ public: * * @return A DataTypeDescriptor */ - const DataTypeDescriptor& dtypeDescr() const; + const DataTypeDescriptor &dtypeDescr() const; /** * @brief The number of frames that is send in a block of DaqData. diff --git a/src/lasp/device/lasp_daqdata.cpp b/src/lasp/device/lasp_daqdata.cpp index bd4c56d..f191ef0 100644 --- a/src/lasp/device/lasp_daqdata.cpp +++ b/src/lasp/device/lasp_daqdata.cpp @@ -1,8 +1,8 @@ /* #define DEBUGTRACE_ENABLED */ -#include "debugtrace.hpp" -#include #include "lasp_daqdata.h" +#include "debugtrace.hpp" #include "lasp_mathtypes.h" +#include #include #include @@ -34,19 +34,20 @@ DaqData::DaqData(const us nframes, const us nchannels, DaqData::DaqData(const DaqData &o) : DaqData(o.nframes, o.nchannels, o.dtype) { DEBUGTRACE_ENTER; - /* std::copy(o._data, &o._data[sw * nchannels * nframes], _data); */ memcpy(_data, o._data, sw * nchannels * nframes); } -/* DaqData::DaqData(DaqData &&o) */ -/* : nframes(o.nframes), nchannels(o.nchannels), dtype(o.dtype), */ -/* dtype_descr(std::move(o.dtype_descr)), sw(o.sw) { */ +DaqData::DaqData(DaqData &&o) + : nframes(o.nframes), nchannels(o.nchannels), dtype(o.dtype), + dtype_descr(std::move(o.dtype_descr)), sw(o.sw) { -/* DEBUGTRACE_ENTER; */ -/* _data = o._data; */ -/* /// Nullptrs do not get deleted */ -/* o._data = nullptr; */ -/* } */ + /// Steal from other one + + DEBUGTRACE_ENTER; + _data = o._data; + /// Nullptrs do not get deleted + o._data = nullptr; +} DaqData::~DaqData() { DEBUGTRACE_ENTER; @@ -58,19 +59,22 @@ void DaqData::copyInFromRaw(const std::vector &ptrs) { DEBUGTRACE_ENTER; us ch = 0; assert(ptrs.size() == nchannels); - for (auto& ptr : ptrs) { - assert(ch < nchannels); - memcpy(&_data[sw * ch * nframes], ptr, sw * nframes); - /* std::copy(ptr, ptr + sw * nframes, &_data[sw * ch * nframes]); */ + for (auto &ptr : ptrs) { + copyInFromRaw(ch, ptr); ch++; } } +void DaqData::copyInFromRaw(const us channel, const byte_t *ptr) { + DEBUGTRACE_ENTER; + assert(ptr); + memcpy(&_data[sw * channel * nframes], ptr, sw * nframes); +} void DaqData::copyToRaw(const us channel, byte_t *ptr) { /* std::copy(raw_ptr(0, channel), raw_ptr(nframes, channel), ptr); */ - assert(channel @@ -78,7 +82,7 @@ d DaqData::toFloat(const us frame, const us channel) const { /* DEBUGTRACE_ENTER; */ if constexpr (std::is_integral::value) { return static_cast(value(frame, channel)) / - std::numeric_limits::max(); + std::numeric_limits::max(); } else { return static_cast(value(frame, channel)); } @@ -87,7 +91,7 @@ d DaqData::toFloat(const us frame, const us channel) const { template vd DaqData::toFloat(const us channel) const { DEBUGTRACE_ENTER; #if LASP_DEBUG == 1 - check_type(); + check_type(); #endif vd res(nframes); for (us i = 0; i < nframes; i++) { @@ -99,7 +103,7 @@ template dmat DaqData::toFloat() const { DEBUGTRACE_ENTER; #if LASP_DEBUG == 1 - check_type(); + check_type(); #endif dmat res(nframes, nchannels); @@ -141,7 +145,7 @@ d DaqData::toFloat(const us frame, const us channel) const { vd DaqData::toFloat(const us channel_no) const { DEBUGTRACE_ENTER; using DataType = DataTypeDescriptor::DataType; - cerr << (int) dtype << endl; + /* cerr << (int)dtype << endl; */ switch (dtype) { case (DataType::dtype_int8): { return toFloat(channel_no); @@ -166,7 +170,6 @@ vd DaqData::toFloat(const us channel_no) const { return vd(); } - dmat DaqData::toFloat() const { DEBUGTRACE_ENTER; @@ -204,17 +207,92 @@ dmat DaqData::toFloat() const { return dmat(); } +template +void DaqData::fromFloat(const us frame, const us channel, const d val) { + DEBUGTRACE_ENTER; +#if LASP_DEBUG == 1 + check_type(); +#endif + if constexpr (std::is_integral::value) { + value(frame, channel) = + static_cast(val * std::numeric_limits::max()); + } else { + value(frame, channel) = static_cast(val); + } +} + +void DaqData::fromFloat(const us frame, const us channel, const d val) { + + using DataType = DataTypeDescriptor::DataType; + + switch (dtype) { + case (DataType::dtype_int8): + return fromFloat(frame, channel, val); + break; + case (DataType::dtype_int16): + return fromFloat(frame, channel, val); + break; + case (DataType::dtype_int32): + return fromFloat(frame, channel, val); + break; + case (DataType::dtype_fl32): + return fromFloat(frame, channel, val); + break; + case (DataType::dtype_fl64): + return fromFloat(frame, channel, val); + break; + default: + throw std::runtime_error("BUG"); + } // End of switch +} +void DaqData::fromFloat(const us channel, const vd &vals) { + if (vals.size() != nframes) { + throw rte("Invalid number of frames in channel data"); + } + using DataType = DataTypeDescriptor::DataType; + switch (dtype) { + case (DataType::dtype_int8): + for (us frame = 0; frame < nframes; frame++) { + fromFloat(frame, channel, vals(frame)); + } + break; + case (DataType::dtype_int16): + for (us frame = 0; frame < nframes; frame++) { + fromFloat(frame, channel, vals(frame)); + } + break; + case (DataType::dtype_int32): + for (us frame = 0; frame < nframes; frame++) { + fromFloat(frame, channel, vals(frame)); + } + break; + case (DataType::dtype_fl32): + for (us frame = 0; frame < nframes; frame++) { + fromFloat(frame, channel, vals(frame)); + } + break; + case (DataType::dtype_fl64): + for (us frame = 0; frame < nframes; frame++) { + fromFloat(frame, channel, vals(frame)); + } + break; + default: + throw std::runtime_error("BUG"); + } // End of switch +} + void DaqData::print() const { cout << "Number of frames: " << nframes << endl; cout << "Number of channels: " << nchannels << endl; cout << "DataType: " << dtype_map.at(dtype).name << endl; - cout << "First sample of first channel (as float)" << toFloat(0,0) << endl; - cout << "Last sample of first channel (as float)" << toFloat(nframes-1,0) << endl; - cout << "Last sample of last channel (as float)" << toFloat(nframes-1,nchannels-1) << endl; + cout << "First sample of first channel (as float)" << toFloat(0, 0) << endl; + cout << "Last sample of first channel (as float)" << toFloat(nframes - 1, 0) + << endl; + cout << "Last sample of last channel (as float)" + << toFloat(nframes - 1, nchannels - 1) << endl; dmat data = toFloat(); vrd max = arma::max(data, 0); vrd min = arma::min(data, 0); cout << "Maximum value in buf: " << max << endl; cout << "Minumum value in buf: " << min << endl; - } diff --git a/src/lasp/device/lasp_daqdata.h b/src/lasp/device/lasp_daqdata.h index 034ebfa..06795f1 100644 --- a/src/lasp/device/lasp_daqdata.h +++ b/src/lasp/device/lasp_daqdata.h @@ -65,9 +65,9 @@ public: * @brief Initialize using no allocation */ DaqData(const DaqData &); - /* DaqData(DaqData &&); */ + DaqData(DaqData &&); DaqData &operator=(const DaqData &) = delete; - virtual ~DaqData(); + ~DaqData(); /** * @brief Return pointer to the raw data corresponding to a certain sample @@ -104,6 +104,15 @@ public: */ void copyInFromRaw(const std::vector &ptrs); + /** + * @brief Copy data from a set of raw pointers of *uninterleaved* data. + * Overwrites any existing available data. + * + * @param channel The channel to copy to + * @param ptr Pointers to data from channels + */ + void copyInFromRaw(const us channel, const byte_t *ptr); + /** * @brief Copy contents of DaqData for a certain channel to a raw pointer. * @@ -126,7 +135,7 @@ public: * column vector of floats. For data that is not already floating-point, * the data is scaled back from MAX_INT to +1.0. * - * @param channel The channel to convert + * @param channel_no The channel number to convert * * @return Array of floats */ @@ -137,13 +146,33 @@ public: * column vector of floats. For data that is not already floating-point, * the data is scaled back from MAX_INT to +1.0. * - * @param channel The channel to convert * @param frame_no The frame number to convert + * @param channel_no The channel number to convert * * @return Float value */ d toFloat(const us frame_no, const us channel_no) const; + /** + * @brief Convert to channel data of native type from floating point values. + * Useful for 'changing' raw data in any way. + * + * @param frame_no The frame + * @param channel_no The channel + * @param data The value + */ + void fromFloat(const us frame_no, const us channel_no, + const d data); + + /** + * @brief Convert to channel data of native type from floating point values. + * Useful for 'changing' raw data in any way. + * + * @param channel The channel to convert + * @param data Data to convert from float values + */ + void fromFloat(const us channel, const arma::Col &data); + // Return value based on type template T &value(const us frame, const us channel) { #if LASP_DEBUG == 1 @@ -188,6 +217,11 @@ protected: template arma::Col toFloat(const us channel_no) const; template d toFloat(const us frame_no, const us channel_no) const; + /* template void fromFloat(const dmat&); */ + template + void fromFloat(const us channel_no, const arma::Col &vals); + template + void fromFloat(const us frame_no, const us channel_no, const d val); /** * @brief Return a value as floating point. Does a conversion from integer to * float, if the stored type is integer. Also scales by its maximum value. diff --git a/src/lasp/device/lasp_device_common.py b/src/lasp/device/lasp_device_common.py deleted file mode 100644 index 2606215..0000000 --- a/src/lasp/device/lasp_device_common.py +++ /dev/null @@ -1,80 +0,0 @@ -__all__ = ['DaqChannel'] -import json, logging -from dataclasses import dataclass, field -from typing import List -import numpy as np -from dataclasses_json import dataclass_json -from ..lasp_common import Qty, SIQtys - - -@dataclass_json -@dataclass(eq=False) -class DaqChannel: - channel_enabled: bool - channel_name: str = 'Unnamed channel' - sensitivity: float = 1.0 - range_index: int = 0 - ACCoupling_enabled: bool = False - IEPE_enabled: bool = False - channel_metadata: str = '' - - def __post_init__(self): - # logging.debug(f'__post_init__({self.channel_name})') - self._qty = SIQtys.default() - - # Whether a digital high-pass filter should be used on input data. - # Negative means disabled. A positive number corresponds to the cut-on - # frequency of the installed highpass filter. - self._highpass = -1.0 - try: - meta = json.loads(self.channel_metadata) - if 'qty' in meta: - # The quantity itself is stored as a JSON string, in the JSON - # object called channel_metadata. - self._qty = Qty.from_json(meta['qty']) - if 'highpass' in meta: - self._highpass = meta['highpass'] - except json.JSONDecodeError: - logging.debug(f'No JSON data found in DaqChannel {self.channel_name}') - - @property - def qty(self): - return self._qty - - @qty.setter - def qty(self, newqty): - self._qty = newqty - self._store('qty', newqty.to_json()) - - @property - def highpass(self): - return self._highpass - - @highpass.setter - def highpass(self, newvalue: float): - newvalue = float(newvalue) - self._highpass = newvalue - self._store('highpass', newvalue) - - def _store(self, name, val): - try: - meta = json.loads(self.channel_metadata) - except json.JSONDecodeError: - meta = {} - - meta[name] = val - self.channel_metadata = json.dumps(meta) - - def __eq__(self, other): - """ - We overwrite the default equal-check as the floating point values - cannot be exactly checked due to conversion from/to JSON - """ - return (self.channel_enabled == other.channel_enabled and - self.channel_name == other.channel_name and - self.range_index == other.range_index and - self.ACCoupling_enabled == other.ACCoupling_enabled and - self.IEPE_enabled == other.IEPE_enabled and - self.channel_metadata == other.channel_metadata and - np.isclose(self.highpass,other.highpass)) - diff --git a/src/lasp/device/lasp_deviceinfo.cpp b/src/lasp/device/lasp_deviceinfo.cpp index ea382bc..f90e34e 100644 --- a/src/lasp/device/lasp_deviceinfo.cpp +++ b/src/lasp/device/lasp_deviceinfo.cpp @@ -12,8 +12,8 @@ #endif -std::vector DeviceInfo::getDeviceInfo() { - std::vector devs; +DeviceInfoList DeviceInfo::getDeviceInfo() { + DeviceInfoList devs; #if LASP_HAS_ULDAQ ==1 fillUlDaqDeviceInfo(devs); #endif diff --git a/src/lasp/device/lasp_deviceinfo.h b/src/lasp/device/lasp_deviceinfo.h index de3a168..81772ed 100644 --- a/src/lasp/device/lasp_deviceinfo.h +++ b/src/lasp/device/lasp_deviceinfo.h @@ -1,17 +1,34 @@ #pragma once #include "lasp_config.h" -#include "lasp_types.h" #include "lasp_daqconfig.h" +#include "lasp_types.h" +#include /** \addtogroup device * @{ */ +using DeviceInfoList = std::vector>; + /** * @brief Structure containing device info parameters */ class DeviceInfo { public: + /** + * @brief Virtual desctructor. Can be derived class. + */ + virtual ~DeviceInfo() {} + + /** + * @brief Clone a device info. + * + * @return Cloned device info + */ + virtual std::unique_ptr clone() const { + return std::make_unique(*this); + } + /** * @brief Backend API corresponding to device */ @@ -21,12 +38,6 @@ public: */ string device_name = ""; - /** - * @brief Specific for the device (Sub-API). Important for the RtAudio - * backend, as RtAudio is able to handle different API's. - */ - int api_specific_devindex = -1; - /** * @brief The available data types the device can output */ @@ -97,10 +108,29 @@ public: } /** - * @brief Whether the device has an internal monitor of the output signal. + * @brief Whether the device has an internal monitor of the output signal. If + * true, the device is able to monitor output signals internally and able to + * present output signals as virtual input signals. This only works together + * Daq's that are able to run in full duplex mode. */ bool hasInternalOutputMonitor = false; + /** + * @brief This flag is used to be able to indicate that the device cannot run + * input and output streams independently, without opening the device in + * duplex mode. This is for example true for the UlDaq: only one handle to + * the device can be given at the same time. + */ + bool duplexModeForced = false; + + /** + * @brief The physical quantity of the output signal. For 'normal' audio + * devices, this is typically a 'number' between +/- full scale. For some + * devices however, the output quantity corresponds to a physical signal, + * such a Volts. + */ + DaqChannel::Qty physicalOutputQty = DaqChannel::Qty::Number; + /** * @brief String representation of DeviceInfo * @@ -109,11 +139,9 @@ public: operator string() const { using std::endl; std::stringstream str; - str << api.apiname + " " << api_specific_devindex << endl - << " number of input channels: " << ninchannels << endl - /** - * @brief - */ + str << api.apiname + " " << endl + << " number of input channels: " << ninchannels + << endl << " number of output channels: " << noutchannels << endl; return str.str(); } @@ -123,7 +151,7 @@ public: * * @return The device info's for each found device. */ - static std::vector getDeviceInfo(); + static DeviceInfoList getDeviceInfo(); }; /** @} */ diff --git a/src/lasp/device/lasp_rtaudiodaq.cpp b/src/lasp/device/lasp_rtaudiodaq.cpp index 411ecda..7c0742c 100644 --- a/src/lasp/device/lasp_rtaudiodaq.cpp +++ b/src/lasp/device/lasp_rtaudiodaq.cpp @@ -17,9 +17,19 @@ using rte = std::runtime_error; using std::vector; using lck = std::scoped_lock; -DEBUGTRACE_VARIABLES; +class RtAudioDeviceInfo : public DeviceInfo { +public: + /** + * @brief Specific for the device (Sub-API). Important for the RtAudio + * backend, as RtAudio is able to handle different API's. + */ + int _api_devindex; + virtual std::unique_ptr clone() const override { + return std::make_unique(*this); + } +}; -void fillRtAudioDeviceInfo(vector &devinfolist) { +void fillRtAudioDeviceInfo(DeviceInfoList &devinfolist) { DEBUGTRACE_ENTER; vector apis; @@ -36,7 +46,7 @@ void fillRtAudioDeviceInfo(vector &devinfolist) { continue; } // "Our device info struct" - DeviceInfo d; + RtAudioDeviceInfo d; switch (api) { case RtAudio::LINUX_ALSA: d.api = rtaudioAlsaApi; @@ -60,13 +70,28 @@ void fillRtAudioDeviceInfo(vector &devinfolist) { } d.device_name = devinfo.name; - d.api_specific_devindex = devno; + d._api_devindex = devno; + + /// When 48k is available we overwrite the default sample rate with the 48 + /// kHz value, which is our preffered rate, + bool rate_48k_found = false; for (us j = 0; j < devinfo.sampleRates.size(); j++) { - us rate = devinfo.sampleRates[j]; - d.availableSampleRates.push_back((double)rate); - if (devinfo.preferredSampleRate == rate) { - d.prefSampleRateIndex = j; + + us rate_int = devinfo.sampleRates[j]; + + d.availableSampleRates.push_back((double)rate_int); + + if (!rate_48k_found) { + + if (devinfo.preferredSampleRate == rate_int) { + d.prefSampleRateIndex = j; + } + + if (rate_int == 48000) { + d.prefSampleRateIndex = j; + rate_48k_found = true; + } } } @@ -103,9 +128,9 @@ void fillRtAudioDeviceInfo(vector &devinfolist) { d.prefDataTypeIndex = d.availableDataTypes.size() - 1; d.availableFramesPerBlock = {512, 1024, 2048, 4096, 8192}; - d.prefFramesPerBlockIndex = 1; + d.prefFramesPerBlockIndex = 2; - devinfolist.push_back(d); + devinfolist.push_back(std::make_unique(d)); } } } @@ -130,9 +155,9 @@ class RtAudioDaq : public Daq { std::atomic _streamStatus{}; public: - RtAudioDaq(const DeviceInfo &devinfo, const DaqConfiguration &config) - : Daq(devinfo, config), - rtaudio(static_cast(devinfo.api.api_specific_subcode)), + RtAudioDaq(const DeviceInfo &devinfo_gen, const DaqConfiguration &config) + : Daq(devinfo_gen, config), rtaudio(static_cast( + devinfo_gen.api.api_specific_subcode)), nFramesPerBlock(Daq::framesPerBlock()) { DEBUGTRACE_ENTER; @@ -143,6 +168,8 @@ public: throw rte("RtAudio backend cannot run in duplex mode."); } assert(!monitorOutput); + const RtAudioDeviceInfo &devinfo = + static_cast(devinfo_gen); std::unique_ptr inParams, outParams; @@ -156,7 +183,7 @@ public: throw rte("Invalid input number of channels"); } inParams->firstChannel = 0; - inParams->deviceId = devinfo.api_specific_devindex; + inParams->deviceId = devinfo._api_devindex; } else { @@ -167,7 +194,7 @@ public: throw rte("Invalid output number of channels"); } outParams->firstChannel = 0; - outParams->deviceId = devinfo.api_specific_devindex; + outParams->deviceId = devinfo._api_devindex; } RtAudio::StreamOptions streamoptions; @@ -218,7 +245,7 @@ public: if (nFramesPerBlock_copy != nFramesPerBlock) { throw rte("Got different number of frames per block back from RtAudio " - "backend. Do not know what to do"); + "backend. I do not know what to do."); } } @@ -397,8 +424,10 @@ std::unique_ptr createRtAudioDevice(const DeviceInfo &devinfo, void myerrorcallback(RtAudioError::Type, const string &errorText) { cerr << "RtAudio backend stream error: " << errorText << endl; } -int mycallback(void *outputBuffer, void *inputBuffer, unsigned int nFrames, - double streamTime, RtAudioStreamStatus status, void *userData) { +int mycallback( + void *outputBuffer, void *inputBuffer, unsigned int nFrames, + __attribute__((unused)) double streamTime, // Not used parameter streamTime + RtAudioStreamStatus status, void *userData) { return static_cast(userData)->streamCallback( outputBuffer, inputBuffer, nFrames, status); diff --git a/src/lasp/device/lasp_rtaudiodaq.h b/src/lasp/device/lasp_rtaudiodaq.h index 47a8cbb..74b6d3e 100644 --- a/src/lasp/device/lasp_rtaudiodaq.h +++ b/src/lasp/device/lasp_rtaudiodaq.h @@ -10,5 +10,5 @@ std::unique_ptr createRtAudioDevice(const DeviceInfo& devinfo, * * @param devinfolist List to append to */ -void fillRtAudioDeviceInfo(std::vector &devinfolist); +void fillRtAudioDeviceInfo(DeviceInfoList &devinfolist); diff --git a/src/lasp/device/lasp_streammgr.cpp b/src/lasp/device/lasp_streammgr.cpp index dea3662..4ae3fa7 100644 --- a/src/lasp/device/lasp_streammgr.cpp +++ b/src/lasp/device/lasp_streammgr.cpp @@ -1,6 +1,7 @@ /* #define DEBUGTRACE_ENABLED */ #include "lasp_streammgr.h" #include "debugtrace.hpp" +#include "lasp_biquadbank.h" #include "lasp_thread.h" #include #include @@ -92,11 +93,46 @@ bool StreamMgr::inCallback(const DaqData &data) { std::scoped_lock lck(_inDataHandler_mtx); - for (auto &handler : _inDataHandlers) { + assert(_inputFilters.size() == data.nchannels); - bool res = handler->inCallback(data); - if (!res) - return false; + if (std::count_if(_inputFilters.cbegin(), _inputFilters.cend(), + [](const auto &a) { return bool(a); }) > 0) { + + /// Found a filter in vector of input filters. So we have to apply the + /// filters to each channel. + + DaqData input_filtered(data.nframes, data.nchannels, data.dtype); + + for (us ch = 0; ch < data.nchannels; ch++) { + if (_inputFilters[ch]) { + DEBUGTRACE_PRINT("Filter ch:"); + DEBUGTRACE_PRINT(ch); + vd inout = data.toFloat(ch); + _inputFilters[ch]->filter(inout); + input_filtered.fromFloat(ch, inout); + } else { + DEBUGTRACE_PRINT("No filter ch:"); + DEBUGTRACE_PRINT(ch); + input_filtered.copyInFromRaw(ch, data.raw_ptr(0, ch)); + } + } + + for (auto &handler : _inDataHandlers) { + bool res = handler->inCallback(input_filtered); + if (!res) { + return false; + } + } + + } else { + /// No input filters + for (auto &handler : _inDataHandlers) { + + bool res = handler->inCallback(data); + if (!res) { + return false; + } + } } return true; } @@ -225,22 +261,23 @@ void StreamMgr::startStream(const DaqConfiguration &config) { [](auto &i) { return i.enabled; }); // Find the first device that matches with the configuration + std::scoped_lock lck(_devices_mtx); - DeviceInfo devinfo; - { - std::scoped_lock lck(_devices_mtx); + DeviceInfo *devinfo = nullptr; + bool found = false; - auto devinfo2 = std::find_if( - _devices.cbegin(), _devices.cend(), - [&config](const DeviceInfo &d) { return config.match(d); }); - if (devinfo2 == std::cend(_devices)) { - throw rte("Could not find a device with name " + config.device_name + - " in list of devices."); + for (auto &devinfoi : _devices) { + if (config.match(*devinfoi)) { + devinfo = devinfoi.get(); + break; } - devinfo = *devinfo2; + } + if (!devinfo) { + throw rte("Could not find a device with name " + config.device_name + + " in list of devices."); } - isInput |= (config.monitorOutput && devinfo.hasInternalOutputMonitor); + isInput |= (config.monitorOutput && devinfo->hasInternalOutputMonitor); DEBUGTRACE_PRINT(isInput); bool isDuplex = isInput && isOutput; @@ -258,48 +295,92 @@ void StreamMgr::startStream(const DaqConfiguration &config) { "first stop existing stream"); } else if (_inputStream) { if (_inputStream->duplexMode() && isOutput) { - throw rte( - "Error: output stream is already running (in duplex mode). Please " - "first stop existing stream"); + throw rte("Error: output stream is already running (in duplex mode). " + "Please " + "first stop existing stream"); } } + if (_outputStream && isInput && _outputStream->duplexModeForced && + config.match(*_outputStream)) { + throw rte("This device is already opened for output. If input is also " + "required, please enable duplex mode for this device"); + } + + if (_inputStream && isOutput && _inputStream->duplexModeForced && + config.match(*_inputStream)) { + throw rte("This device is already opened for input. If output is also " + "required, please enable duplex mode for this device"); + } + InDaqCallback inCallback; OutDaqCallback outCallback; using namespace std::placeholders; - std::unique_ptr daq = Daq::createDaq(devinfo, config); + std::unique_ptr daq = Daq::createDaq(*devinfo, config); + + assert(daq); - std::unique_ptr *stream_placeholder; if (isInput) { + /// Give incallback as parameter to stream inCallback = std::bind(&StreamMgr::inCallback, this, _1); - stream_placeholder = std::addressof(_inputStream); + + /// Reset handlers in case of an input stream for (auto &handler : _inDataHandlers) { handler->reset(daq.get()); } + + d fs = daq->samplerate(); + /// Create input filters + _inputFilters.clear(); + /// No input filter for monitor channel. + if (config.monitorOutput && devinfo->hasInternalOutputMonitor) { + _inputFilters.push_back(nullptr); + } + for (auto &ch : daq->inchannel_config) { + if (ch.enabled) { + if (ch.digitalHighPassCutOn < 0) { + _inputFilters.push_back(nullptr); + } else if (ch.digitalHighPassCutOn == 0) { + throw rte("Digital highpass cuton should be > 0 if activated"); + } else { + // Put in a digital high-pass filter. + _inputFilters.emplace_back(std::make_unique( + SeriesBiquad::firstOrderHighPass(fs, ch.digitalHighPassCutOn))); + } + } + } // End of input filter creation } + if (isOutput) { + /// Give outcallback as parameter to stream + outCallback = std::bind(&StreamMgr::outCallback, this, _1); + + /// Reset signal generator in case of an output stream if (_siggen) { DEBUGTRACE_PRINT("Resetting _siggen with new samplerate of "); DEBUGTRACE_PRINT(daq->samplerate()); _siggen->reset(daq->samplerate()); } - outCallback = std::bind(&StreamMgr::outCallback, this, _1); - if(!isDuplex) { - stream_placeholder = std::addressof(_outputStream); - } } + /// Start the DAQ. If it fails, everything is still nicely cleaned up and + /// the daq unique_ptr cleans up resources nicely. daq->start(inCallback, outCallback); // Move daq ptr to right place - *stream_placeholder = std::move(daq); + if (isInput) { + _inputStream = std::move(daq); + } else { + _outputStream = std::move(daq); + } } void StreamMgr::stopStream(const StreamType t) { + DEBUGTRACE_ENTER; checkRightThread(); - switch (t) { - case (StreamType::input): { + + if (t == StreamType::input) { if (!_inputStream) { throw rte("Input stream is not running"); } @@ -309,21 +390,18 @@ void StreamMgr::stopStream(const StreamType t) { for (auto &handler : _inDataHandlers) { handler->reset(nullptr); } - } break; - case (StreamType::output): { + } else { + /// t == output + + /// Kill input stream in case that one is a duplex stream if (_inputStream && _inputStream->duplexMode()) { _inputStream.reset(); } else { - if (!_outputStream) { throw rte("Output stream is not running"); } _outputStream.reset(); } // end else - } break; - default: - throw rte("BUG"); - break; } } diff --git a/src/lasp/device/lasp_streammgr.h b/src/lasp/device/lasp_streammgr.h index ce155a1..fc69222 100644 --- a/src/lasp/device/lasp_streammgr.h +++ b/src/lasp/device/lasp_streammgr.h @@ -65,6 +65,8 @@ public: void stop(); }; +class SeriesBiquad; + /** * @brief Stream manager. Used to manage the input and output streams. * Implemented as a singleton: only one stream manager can be in existance for @@ -96,10 +98,16 @@ class StreamMgr { * implemented as to generate the same data for all output channels. */ std::shared_ptr _siggen; + + /** + * @brief Filters on input stream. For example, a digital high pass filter. + */ + std::vector> _inputFilters; + std::mutex _siggen_mtx; std::mutex _devices_mtx; - std::vector _devices; + DeviceInfoList _devices; StreamMgr(); @@ -137,9 +145,13 @@ class StreamMgr { * * @return A copy of the internal stored list of devices */ - std::vector getDeviceInfo() const { + DeviceInfoList getDeviceInfo() const { std::scoped_lock lck(const_cast(_devices_mtx)); - return _devices; + DeviceInfoList d2; + for(const auto& dev: _devices) { + d2.push_back(dev->clone()); + } + return d2; } /** diff --git a/src/lasp/device/lasp_uldaq.cpp b/src/lasp/device/lasp_uldaq.cpp index 8f0d4c4..41de7dc 100644 --- a/src/lasp/device/lasp_uldaq.cpp +++ b/src/lasp/device/lasp_uldaq.cpp @@ -1,593 +1,22 @@ /* #define DEBUGTRACE_ENABLED */ #include "debugtrace.hpp" #include "lasp_config.h" - #if LASP_HAS_ULDAQ == 1 -#include "lasp_daqconfig.h" #include "lasp_uldaq.h" -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include +#include "lasp_uldaq_impl.h" #include -#include -using namespace std::literals::chrono_literals; -using std::atomic; -using std::cerr; -using std::endl; -using rte = std::runtime_error; -#include "debugtrace.hpp" -DEBUGTRACE_VARIABLES; -const us MAX_DEV_COUNT_PER_API = 100; -/** - * @brief Reserve some space for an error message from UlDaq - */ -const us UL_ERR_MSG_LEN = 512; - -/** - * @brief Show the error to default error stream and return a string - * corresponding to the error - * - * @param err Error string - */ -string showErr(UlError err) { - string errstr; - errstr.reserve(UL_ERR_MSG_LEN); - if (err != ERR_NO_ERROR) { - char errmsg[UL_ERR_MSG_LEN]; - errstr = "UlDaq API Error: "; - ulGetErrMsg(err, errmsg); - errstr += errmsg; - - std::cerr << "\b\n**************** UlDAQ backend error **********\n"; - std::cerr << errstr << std::endl; - std::cerr << "***********************************************\n\n"; - return errstr; - } - return errstr; -} - -class DT9837A : public Daq { - - DaqDeviceHandle _handle = 0; - std::mutex _daqmutex; - - std::thread _thread; - atomic _stopThread{false}; - atomic _streamStatus; - - const us _nFramesPerBlock; - - void threadFcn(InDaqCallback inCallback, OutDaqCallback outcallback); - -public: - DaqDeviceHandle getHandle() const { return _handle; } - DT9837A(const DeviceInfo &devinfo, const DaqConfiguration &config); - - ~DT9837A() { - UlError err; - if (isRunning()) { - stop(); - } - - if (_handle) { - err = ulDisconnectDaqDevice(_handle); - showErr(err); - err = ulReleaseDaqDevice(_handle); - showErr(err); - } - } - - bool isRunning() const { - DEBUGTRACE_ENTER; - return _thread.joinable(); - } - virtual void start(InDaqCallback inCallback, - OutDaqCallback outCallback) override final; - - virtual StreamStatus getStreamStatus() const override { - return _streamStatus; - } - - void stop() override final { - DEBUGTRACE_ENTER; - StreamStatus status = _streamStatus; - if (!isRunning()) { - throw rte("No data acquisition running"); - } - - _stopThread = true; - if (_thread.joinable()) { - _thread.join(); - } - _stopThread = false; - status.isRunning = false; - _streamStatus = status; - } - - friend class InBufHandler; - friend class OutBufHandler; -}; - -void DT9837A::start(InDaqCallback inCallback, OutDaqCallback outCallback) { - DEBUGTRACE_ENTER; - if (isRunning()) { - throw rte("DAQ is already running"); - } - if (neninchannels() > 0) { - if (!inCallback) - throw rte("DAQ requires a callback for input data"); - } - if (nenoutchannels() > 0) { - if (!outCallback) - throw rte("DAQ requires a callback for output data"); - } - assert(neninchannels() + nenoutchannels() > 0); - _thread = std::thread(&DT9837A::threadFcn, this, inCallback, outCallback); -} - -class BufHandler { -protected: - DT9837A &daq; - const DataTypeDescriptor dtype_descr; - us nchannels, nFramesPerBlock; - double samplerate; - std::vector buf; - bool topenqueued, botenqueued; - - us increment = 0; - - us totalFramesCount = 0; - long long buffer_mid_idx; - -public: - BufHandler(DT9837A &daq, const us nchannels) - : daq(daq), dtype_descr(daq.dtypeDescr()), nchannels(nchannels), - nFramesPerBlock(daq.framesPerBlock()), samplerate(daq.samplerate()), - buf(2 * nchannels * - nFramesPerBlock, // Watch the two here, the top and the bottom! - 0), - buffer_mid_idx(nchannels * nFramesPerBlock) { - assert(nchannels > 0); - } -}; -class InBufHandler : public BufHandler { - bool monitorOutput; - InDaqCallback cb; - -public: - InBufHandler(DT9837A &daq, InDaqCallback cb) - : BufHandler(daq, daq.neninchannels()), cb(cb) - - { - DEBUGTRACE_ENTER; - assert(daq.getHandle() != 0); - - monitorOutput = daq.monitorOutput; - - DaqInScanFlag inscanflags = DAQINSCAN_FF_DEFAULT; - ScanOption scanoptions = SO_CONTINUOUS; - UlError err = ERR_NO_ERROR; - - std::vector indescs; - boolvec eninchannels_without_mon = daq.eninchannels(false); - - // Initialize input, if any - dvec ranges = daq.inputRangeForEnabledChannels(false); - for (us chin = 0; chin < 4; chin++) { - if (eninchannels_without_mon[chin] == true) { - DaqInChanDescriptor indesc; - indesc.type = DAQI_ANALOG_SE; - indesc.channel = chin; - - double rangeval = ranges.at(chin); - Range rangenum; - if (fabs(rangeval - 1.0) < 1e-8) { - rangenum = BIP1VOLTS; - } else if (fabs(rangeval - 10.0) < 1e-8) { - rangenum = BIP10VOLTS; - } else { - std::cerr << "Fatal: input range value is invalid" << endl; - return; - } - indesc.range = rangenum; - indescs.push_back(indesc); - } - } - - // Add possibly last channel as monitor - if (monitorOutput) { - DaqInChanDescriptor indesc; - indesc.type = DAQI_DAC; - indesc.channel = 0; - /// The output only has a range of 10V, therefore the monitor of the - /// output also has to be set to this value. - indesc.range = BIP10VOLTS; - indescs.push_back(indesc); - } - assert(indescs.size() == nchannels); - - DEBUGTRACE_MESSAGE("Starting input scan"); - - err = ulDaqInScan(daq.getHandle(), indescs.data(), nchannels, - 2 * nFramesPerBlock, // Watch the 2 here! - &samplerate, scanoptions, inscanflags, buf.data()); - if (err != ERR_NO_ERROR) { - showErr(err); - throw rte("Could not start input DAQ"); - } - } - void start() { - - ScanStatus status; - TransferStatus transferStatus; - UlError err = ulDaqInScanStatus(daq.getHandle(), &status, &transferStatus); - if (err != ERR_NO_ERROR) { - showErr(err); - throw rte("Unable to start input on DAQ"); - } - - totalFramesCount = transferStatus.currentTotalCount; - topenqueued = true; - botenqueued = true; - } - - /** - * @brief InBufHandler::operator()() - * - * @return true on success - */ - bool operator()() { - - /* DEBUGTRACE_ENTER; */ - - bool ret = true; - - auto runCallback = ([&](us totalOffset) { - /* DEBUGTRACE_ENTER; */ - - DaqData data(nFramesPerBlock, nchannels, - DataTypeDescriptor::DataType::dtype_fl64); - - us monitorOffset = monitorOutput ? 1 : 0; - /* /// Put the output monitor in front */ - if (monitorOutput) { - for (us frame = 0; frame < nFramesPerBlock; frame++) { - data.value(frame, 0) = - buf[totalOffset + (frame * nchannels) + (nchannels - 1)]; - } - } - - for (us channel = 0; channel < nchannels - monitorOffset; channel++) { - /* DEBUGTRACE_PRINT(channel); */ - for (us frame = 0; frame < nFramesPerBlock; frame++) { - data.value(frame, channel + monitorOffset) = - buf[totalOffset + (frame * nchannels) + channel]; - } - } - return cb(data); - }); - - ScanStatus status; - TransferStatus transferStatus; - - UlError err = ulDaqInScanStatus(daq.getHandle(), &status, &transferStatus); - if (err != ERR_NO_ERROR) { - showErr(err); - return false; - } - - increment = transferStatus.currentTotalCount - totalFramesCount; - totalFramesCount += increment; - - if (increment > nFramesPerBlock) { - cerr << "Error: overrun for input of DAQ!" << endl; - return false; - } - assert(status == SS_RUNNING); - - if (transferStatus.currentIndex < (long long)buffer_mid_idx) { - topenqueued = false; - if (!botenqueued) { - ret = runCallback(nchannels * nFramesPerBlock); - botenqueued = true; - } - } else { - botenqueued = false; - if (!topenqueued) { - ret = runCallback(0); - topenqueued = true; - } - } - return ret; - } - ~InBufHandler() { - // At exit of the function, stop scanning. - DEBUGTRACE_ENTER; - UlError err = ulDaqInScanStop(daq.getHandle()); - if (err != ERR_NO_ERROR) { - showErr(err); - } - } -}; - -class OutBufHandler : public BufHandler { - OutDaqCallback cb; - -public: - OutBufHandler(DT9837A &daq, OutDaqCallback cb) - : BufHandler(daq, daq.nenoutchannels()), cb(cb) { - - DEBUGTRACE_MESSAGE("Starting output scan"); - DEBUGTRACE_PRINT(nchannels); - AOutScanFlag outscanflags = AOUTSCAN_FF_DEFAULT; - ScanOption scanoptions = SO_CONTINUOUS; - UlError err = - ulAOutScan(daq.getHandle(), 0, 0, BIP10VOLTS, - 2 * nFramesPerBlock, // Watch the 2 here! - &samplerate, scanoptions, outscanflags, buf.data()); - - if (err != ERR_NO_ERROR) { - showErr(err); - throw rte("Unable to start output on DAQ"); - } - } - void start() { - - ScanStatus status; - TransferStatus transferStatus; - - UlError err = ulAOutScanStatus(daq.getHandle(), &status, &transferStatus); - if (err != ERR_NO_ERROR) { - showErr(err); - throw rte("Unable to start output on DAQ"); - } - if (status != SS_RUNNING) { - throw rte("Unable to start output on DAQ"); - } - totalFramesCount = transferStatus.currentTotalCount; - topenqueued = true; - botenqueued = true; - } - /** - * @brief OutBufHandler::operator() - * - * @return true on success - */ - bool operator()() { - - /* DEBUGTRACE_ENTER; */ - bool res = true; - assert(daq.getHandle() != 0); - - UlError err = ERR_NO_ERROR; - - ScanStatus status; - TransferStatus transferStatus; - - err = ulAOutScanStatus(daq.getHandle(), &status, &transferStatus); - if (err != ERR_NO_ERROR) { - showErr(err); - return false; - } - if (status != SS_RUNNING) { - return false; - } - increment = transferStatus.currentTotalCount - totalFramesCount; - totalFramesCount += increment; - - if (increment > nFramesPerBlock) { - cerr << "Error: underrun for output of DAQ!" << endl; - return false; - } - - if (transferStatus.currentIndex < buffer_mid_idx) { - topenqueued = false; - if (!botenqueued) { - DaqData d(nFramesPerBlock,1, - DataTypeDescriptor::DataType::dtype_fl64); - res = cb(d); - d.copyToRaw(0, reinterpret_cast(&(buf[buffer_mid_idx]))); - - botenqueued = true; - } - } else { - botenqueued = false; - if (!topenqueued) { - DaqData d(nFramesPerBlock,1, - DataTypeDescriptor::DataType::dtype_fl64); - res = cb(d); - d.copyToRaw(0, reinterpret_cast(&(buf[0]))); - - topenqueued = true; - } - } - return res; - } - - ~OutBufHandler() { - DEBUGTRACE_ENTER; - UlError err = ulAOutScanStop(daq.getHandle()); - if (err != ERR_NO_ERROR) { - showErr(err); - } - } -}; - -void DT9837A::threadFcn(InDaqCallback inCallback, OutDaqCallback outCallback) { +void fillUlDaqDeviceInfo(DeviceInfoList &devinfolist) { DEBUGTRACE_ENTER; - /* cerr << "******************\n" */ - /* "Todo: the current way of handling timing in this DAQ thread is not " */ - /* "really robust, due " */ - /* "to input / output callbacks that can be too time-consuming. We have " */ - /* "to fix the " */ - /* "sleep_for to properly deal with longer callbacks." */ - /* "\n*****************" */ - /* << endl; */ - try { - - std::unique_ptr obh; - std::unique_ptr ibh; - - StreamStatus status = _streamStatus; - status.isRunning = true; - _streamStatus = status; - - if (nenoutchannels() > 0) { - assert(outCallback); - obh = std::make_unique(*this, outCallback); - } - if (neninchannels() > 0) { - assert(inCallback); - ibh = std::make_unique(*this, inCallback); - } - if (obh) - obh->start(); - if (ibh) - ibh->start(); - - const double sleeptime_s = - static_cast(_nFramesPerBlock) / (16 * samplerate()); - const us sleeptime_us = static_cast(sleeptime_s * 1e6); - - while (!_stopThread) { - if (ibh) { - if (!(*ibh)()) { - _stopThread = true; - } - } - if (obh) { - if (!(*obh)()) { - _stopThread = true; - } - } else { - std::this_thread::sleep_for(std::chrono::microseconds(sleeptime_us)); - } - } - } catch (rte &e) { - - StreamStatus status = _streamStatus; - status.isRunning = false; - status.errorType = StreamStatus::StreamError::systemError; - _streamStatus = status; - - cerr << "\n******************\n"; - cerr << "Catched error in UlDAQ thread: " << e.what() << endl; - cerr << "\n******************\n"; - } - StreamStatus status = _streamStatus; - - status.isRunning = false; - _streamStatus = status; - _stopThread = false; -} - -std::unique_ptr createUlDaqDevice(const DeviceInfo &devinfo, - const DaqConfiguration &config) { - return std::make_unique(devinfo, config); -} - -DT9837A::DT9837A(const DeviceInfo &devinfo, const DaqConfiguration &config) - : Daq(devinfo, config), - _nFramesPerBlock(availableFramesPerBlock.at(framesPerBlockIndex)) { - - // Some sanity checks - if (inchannel_config.size() != 4) { - throw rte("Invalid length of enabled inChannels vector"); - } - - if (outchannel_config.size() != 1) { - throw rte("Invalid length of enabled outChannels vector"); - } - - if (_nFramesPerBlock < 24 || _nFramesPerBlock > 8192) { - throw rte("Unsensible number of samples per block chosen"); - } - - if (samplerate() < 10000 || samplerate() > 51000) { - throw rte("Invalid sample rate"); - } - - DaqDeviceDescriptor devdescriptors[MAX_DEV_COUNT_PER_API]; - DaqDeviceDescriptor descriptor; - DaqDeviceInterface interfaceType = ANY_IFC; UlError err; + unsigned int numdevs = MAX_ULDAQ_DEV_COUNT_PER_API; - us numdevs = MAX_DEV_COUNT_PER_API; - err = ulGetDaqDeviceInventory(interfaceType, devdescriptors, - (unsigned *)&numdevs); - if (err != ERR_NO_ERROR) { - throw rte("Device inventarization failed"); - } - - if ((us)api_specific_devindex >= numdevs) { - throw rte("Device number {deviceno} too high {err}. This could " - "happen when the device is currently not connected"); - } - - descriptor = devdescriptors[api_specific_devindex]; - - // get a handle to the DAQ device associated with the first descriptor - _handle = ulCreateDaqDevice(descriptor); - - if (_handle == 0) { - throw rte("Unable to create a handle to the specified DAQ " - "device. Is the device currently in use? Please make sure to set " - "the DAQ configuration in duplex mode if simultaneous input and " - "output is required."); - } - - err = ulConnectDaqDevice(_handle); - if (err != ERR_NO_ERROR) { - ulReleaseDaqDevice(_handle); - _handle = 0; - throw rte(string("Unable to connect to device: " + showErr(err))); - } - - for (us ch = 0; ch < 4; ch++) { - - err = ulAISetConfigDbl(_handle, AI_CFG_CHAN_SENSOR_SENSITIVITY, ch, 1.0); - showErr(err); - if (err != ERR_NO_ERROR) { - throw rte("Fatal: could normalize channel sensitivity"); - } - - CouplingMode cm = inchannel_config.at(ch).ACCouplingMode ? CM_AC : CM_DC; - err = ulAISetConfig(_handle, AI_CFG_CHAN_COUPLING_MODE, ch, cm); - if (err != ERR_NO_ERROR) { - showErr(err); - throw rte("Fatal: could not set AC/DC coupling mode"); - } - - IepeMode iepe = - inchannel_config.at(ch).IEPEEnabled ? IEPE_ENABLED : IEPE_DISABLED; - err = ulAISetConfig(_handle, AI_CFG_CHAN_IEPE_MODE, ch, iepe); - if (err != ERR_NO_ERROR) { - showErr(err); - throw rte("Fatal: could not set IEPE mode"); - } - } -} -void fillUlDaqDeviceInfo(std::vector &devinfolist) { - - DEBUGTRACE_ENTER; - - UlError err; - unsigned int numdevs = MAX_DEV_COUNT_PER_API; - - DaqDeviceDescriptor devdescriptors[MAX_DEV_COUNT_PER_API]; + DaqDeviceDescriptor devdescriptors[MAX_ULDAQ_DEV_COUNT_PER_API]; DaqDeviceDescriptor descriptor; DaqDeviceInterface interfaceType = ANY_IFC; @@ -602,11 +31,14 @@ void fillUlDaqDeviceInfo(std::vector &devinfolist) { descriptor = devdescriptors[i]; - DeviceInfo devinfo; + UlDaqDeviceInfo devinfo; + devinfo._uldaqDescriptor = descriptor; + devinfo.api = uldaqapi; string name, interface; - if (string(descriptor.productName) != "DT9837A") { - throw rte("Unknown UlDAQ type"); + string productname = descriptor.productName; + if (productname != "DT9837A") { + throw rte("Unknown UlDAQ type: " + productname); } switch (descriptor.devInterface) { @@ -629,7 +61,8 @@ void fillUlDaqDeviceInfo(std::vector &devinfolist) { name += string(descriptor.productName) + " " + string(descriptor.uniqueId); devinfo.device_name = std::move(name); - devinfo.api_specific_devindex = i; + devinfo.physicalOutputQty = DaqChannel::Qty::Voltage; + devinfo.availableDataTypes.push_back( DataTypeDescriptor::DataType::dtype_fl64); devinfo.prefDataTypeIndex = 0; @@ -654,8 +87,16 @@ void fillUlDaqDeviceInfo(std::vector &devinfolist) { devinfo.hasInternalOutputMonitor = true; + devinfo.duplexModeForced = true; + // Finally, this devinfo is pushed back in list - devinfolist.push_back(devinfo); + devinfolist.push_back(std::make_unique(devinfo)); } } + +std::unique_ptr createUlDaqDevice(const DeviceInfo &devinfo, + const DaqConfiguration &config) { + return std::make_unique(devinfo, config); +} + #endif // LASP_HAS_ULDAQ diff --git a/src/lasp/device/lasp_uldaq.h b/src/lasp/device/lasp_uldaq.h index 4dd7489..26e2305 100644 --- a/src/lasp/device/lasp_uldaq.h +++ b/src/lasp/device/lasp_uldaq.h @@ -1,14 +1,18 @@ #pragma once #include "lasp_daq.h" -std::unique_ptr createUlDaqDevice(const DeviceInfo& devinfo, - const DaqConfiguration& config); +/** + * @brief The maximum number of devices that can be enumerated when calling + * ulGetDaqDeviceInventory() + */ +const us MAX_ULDAQ_DEV_COUNT_PER_API = 100; + +std::unique_ptr createUlDaqDevice(const DeviceInfo &devinfo, + const DaqConfiguration &config); /** * @brief Fill device info list with UlDaq specific devices, if any. * - * @param devinfolist Info list to append to + * @param devinfolist Info list to append to. */ -void fillUlDaqDeviceInfo(std::vector &devinfolist); - - +void fillUlDaqDeviceInfo(DeviceInfoList& devinfolist); diff --git a/src/lasp/device/lasp_uldaq_impl.cpp b/src/lasp/device/lasp_uldaq_impl.cpp new file mode 100644 index 0000000..1afe301 --- /dev/null +++ b/src/lasp/device/lasp_uldaq_impl.cpp @@ -0,0 +1,480 @@ +/* #define DEBUGTRACE_ENABLED */ +#include "debugtrace.hpp" +#include "lasp_config.h" + +#if LASP_HAS_ULDAQ == 1 +#include "lasp_daqconfig.h" +#include "lasp_uldaq.h" +#include "lasp_uldaq_impl.h" + +using namespace std::literals::chrono_literals; + +/** + * @brief Reserve some space for an error message from UlDaq + */ +const us UL_ERR_MSG_LEN = 512; + +/** + * @brief Return a string corresponding to the UlDaq API error + * + * @param err error code + * + * @return Error string + */ +string getErrMsg(UlError err) { + string errstr; + errstr.reserve(UL_ERR_MSG_LEN); + char errmsg[UL_ERR_MSG_LEN]; + errstr = "UlDaq API Error: "; + ulGetErrMsg(err, errmsg); + errstr += errmsg; + return errstr; +} +inline void showErr(string errstr) { + std::cerr << "\b\n**************** UlDAQ backend error **********\n"; + std::cerr << errstr << std::endl; + std::cerr << "***********************************************\n\n"; +} +inline void showErr(UlError err) { + if (err != ERR_NO_ERROR) + showErr(getErrMsg(err)); +} + +DT9837A::~DT9837A() { + UlError err; + if (isRunning()) { + stop(); + } + + if (_handle) { + err = ulDisconnectDaqDevice(_handle); + showErr(err); + err = ulReleaseDaqDevice(_handle); + showErr(err); + } +} +DT9837A::DT9837A(const DeviceInfo &devinfo, const DaqConfiguration &config) + : Daq(devinfo, config), + _nFramesPerBlock(availableFramesPerBlock.at(framesPerBlockIndex)) { + + // Some sanity checks + if (inchannel_config.size() != 4) { + throw rte("Invalid length of enabled inChannels vector"); + } + + if (outchannel_config.size() != 1) { + throw rte("Invalid length of enabled outChannels vector"); + } + + if (_nFramesPerBlock < 24 || _nFramesPerBlock > 8192) { + throw rte("Unsensible number of samples per block chosen"); + } + + if (samplerate() < 10000 || samplerate() > 51000) { + throw rte("Invalid sample rate"); + } + + const UlDaqDeviceInfo *_info = + dynamic_cast(&devinfo); + if (_info == nullptr) { + throw rte("BUG: Could not cast DeviceInfo to UlDaqDeviceInfo"); + } + + // get a handle to the DAQ device associated with the first descriptor + _handle = ulCreateDaqDevice(_info->_uldaqDescriptor); + + if (_handle == 0) { + throw rte("Unable to create a handle to the specified DAQ " + "device. Is the device currently in use? Please make sure to set " + "the DAQ configuration in duplex mode if simultaneous input and " + "output is required."); + } + + UlError err = ulConnectDaqDevice(_handle); + + if (err != ERR_NO_ERROR) { + ulReleaseDaqDevice(_handle); + _handle = 0; + throw rte("Unable to connect to device: " + getErrMsg(err)); + } + + /// Loop over input channels, set parameters + for (us ch = 0; ch < 4; ch++) { + + err = ulAISetConfigDbl(_handle, AI_CFG_CHAN_SENSOR_SENSITIVITY, ch, 1.0); + showErr(err); + if (err != ERR_NO_ERROR) { + throw rte("Fatal: could normalize channel sensitivity"); + } + + CouplingMode cm = inchannel_config.at(ch).ACCouplingMode ? CM_AC : CM_DC; + err = ulAISetConfig(_handle, AI_CFG_CHAN_COUPLING_MODE, ch, cm); + if (err != ERR_NO_ERROR) { + showErr(err); + throw rte("Fatal: could not set AC/DC coupling mode"); + } + + IepeMode iepe = + inchannel_config.at(ch).IEPEEnabled ? IEPE_ENABLED : IEPE_DISABLED; + err = ulAISetConfig(_handle, AI_CFG_CHAN_IEPE_MODE, ch, iepe); + if (err != ERR_NO_ERROR) { + showErr(err); + throw rte("Fatal: could not set IEPE mode"); + } + } +} + +bool DT9837A::isRunning() const { + DEBUGTRACE_ENTER; + return _thread.joinable(); +} +void DT9837A::stop() { + DEBUGTRACE_ENTER; + StreamStatus status = _streamStatus; + if (!isRunning()) { + throw rte("No data acquisition running"); + } + + _stopThread = true; + if (_thread.joinable()) { + _thread.join(); + } + _stopThread = false; + status.isRunning = false; + _streamStatus = status; +} + +/** + * @brief Throws an appropriate stream exception based on the UlError number. + * The mapping is based on the error numbers as given in uldaq.h. There are a + * log of errors definded here (109 in total). Except for some, we will map + * most of them to a driver error. + * + * @param e + */ +inline void throwUlException(UlError err) { + if (err == ERR_NO_ERROR) { + return; + } + string errstr = getErrMsg(err); + showErr(errstr); + Daq::StreamStatus::StreamError serr; + if ((int)err == 18) { + serr = Daq::StreamStatus::StreamError::inputXRun; + } else if ((int)err == 19) { + serr = Daq::StreamStatus::StreamError::outputXRun; + } else { + serr = Daq::StreamStatus::StreamError::driverError; + } + + throw Daq::StreamException(serr, errstr); +} + +void DT9837A::start(InDaqCallback inCallback, OutDaqCallback outCallback) { + DEBUGTRACE_ENTER; + if (isRunning()) { + throw rte("DAQ is already running"); + } + if (neninchannels() > 0) { + if (!inCallback) + throw rte("DAQ requires a callback for input data"); + } + if (nenoutchannels() > 0) { + if (!outCallback) + throw rte("DAQ requires a callback for output data"); + } + assert(neninchannels() + nenoutchannels() > 0); + _thread = std::thread(&DT9837A::threadFcn, this, inCallback, outCallback); +} + +InBufHandler::InBufHandler(DT9837A &daq, InDaqCallback cb) + : BufHandler(daq, daq.neninchannels()), cb(cb) + +{ + DEBUGTRACE_ENTER; + assert(daq.getHandle() != 0); + + monitorOutput = daq.monitorOutput; + + DaqInScanFlag inscanflags = DAQINSCAN_FF_DEFAULT; + ScanOption scanoptions = SO_CONTINUOUS; + UlError err = ERR_NO_ERROR; + + std::vector indescs; + boolvec eninchannels_without_mon = daq.eninchannels(false); + + // Set ranges for each input. Below asks only channels that are not a + // monitor channel (hence the false flag). + dvec ranges = daq.inputRangeForEnabledChannels(false); + for (us chin = 0; chin < 4; chin++) { + if (eninchannels_without_mon[chin] == true) { + DaqInChanDescriptor indesc; + indesc.type = DAQI_ANALOG_SE; + indesc.channel = chin; + + double rangeval = ranges.at(chin); + Range rangenum; + if (fabs(rangeval - 1.0) < 1e-8) { + rangenum = BIP1VOLTS; + } else if (fabs(rangeval - 10.0) < 1e-8) { + rangenum = BIP10VOLTS; + } else { + throw Daq::StreamException(Daq::StreamStatus::StreamError::logicError); + std::cerr << "Fatal: input range value is invalid" << endl; + return; + } + indesc.range = rangenum; + indescs.push_back(indesc); + } + } + + // Add possibly last channel as monitor + if (monitorOutput) { + DaqInChanDescriptor indesc; + indesc.type = DAQI_DAC; + indesc.channel = 0; + /// The output only has a range of 10V, therefore the monitor of the + /// output also has to be set to this value. + indesc.range = BIP10VOLTS; + indescs.push_back(indesc); + } + assert(indescs.size() == nchannels); + + DEBUGTRACE_MESSAGE("Starting input scan"); + + err = ulDaqInScan(daq.getHandle(), indescs.data(), nchannels, + 2 * nFramesPerBlock, // Watch the 2 here! + &samplerate, scanoptions, inscanflags, buf.data()); + throwUlException(err); +} +void InBufHandler::start() { + + ScanStatus status; + TransferStatus transferStatus; + UlError err = ulDaqInScanStatus(daq.getHandle(), &status, &transferStatus); + throwUlException(err); + + totalFramesCount = transferStatus.currentTotalCount; + topenqueued = true; + botenqueued = true; +} + +bool InBufHandler::operator()() { + + /* DEBUGTRACE_ENTER; */ + + bool ret = true; + + auto runCallback = ([&](us totalOffset) { + /* DEBUGTRACE_ENTER; */ + + DaqData data(nFramesPerBlock, nchannels, dtype_descr.dtype); + + us monitorOffset = monitorOutput ? 1 : 0; + /* /// Put the output monitor in front */ + if (monitorOutput) { + for (us frame = 0; frame < nFramesPerBlock; frame++) { + data.value(frame, 0) = + buf[totalOffset + (frame * nchannels) + (nchannels - 1)]; + } + } + + for (us channel = 0; channel < nchannels - monitorOffset; channel++) { + /* DEBUGTRACE_PRINT(channel); */ + for (us frame = 0; frame < nFramesPerBlock; frame++) { + data.value(frame, channel + monitorOffset) = + buf[totalOffset + (frame * nchannels) + channel]; + } + } + return cb(data); + }); + + ScanStatus status; + TransferStatus transferStatus; + + UlError err = ulDaqInScanStatus(daq.getHandle(), &status, &transferStatus); + throwUlException(err); + + us increment = transferStatus.currentTotalCount - totalFramesCount; + totalFramesCount += increment; + + if (increment > nFramesPerBlock) { + throw Daq::StreamException(Daq::StreamStatus::StreamError::inputXRun); + } + assert(status == SS_RUNNING); + + if (transferStatus.currentIndex < (long long)buffer_mid_idx) { + topenqueued = false; + if (!botenqueued) { + ret = runCallback(nchannels * nFramesPerBlock); + botenqueued = true; + } + } else { + botenqueued = false; + if (!topenqueued) { + ret = runCallback(0); + topenqueued = true; + } + } + return ret; +} +InBufHandler::~InBufHandler() { + // At exit of the function, stop scanning. + DEBUGTRACE_ENTER; + UlError err = ulDaqInScanStop(daq.getHandle()); + if (err != ERR_NO_ERROR) { + showErr(err); + } +} + +OutBufHandler::OutBufHandler(DT9837A &daq, OutDaqCallback cb) + : BufHandler(daq, daq.nenoutchannels()), cb(cb) { + + DEBUGTRACE_MESSAGE("Starting output scan"); + DEBUGTRACE_PRINT(nchannels); + AOutScanFlag outscanflags = AOUTSCAN_FF_DEFAULT; + ScanOption scanoptions = SO_CONTINUOUS; + UlError err = ulAOutScan(daq.getHandle(), 0, 0, BIP10VOLTS, + 2 * nFramesPerBlock, // Watch the 2 here! + &samplerate, scanoptions, outscanflags, buf.data()); + + throwUlException(err); +} +void OutBufHandler::start() { + + ScanStatus status; + TransferStatus transferStatus; + + UlError err = ulAOutScanStatus(daq.getHandle(), &status, &transferStatus); + if (err != ERR_NO_ERROR) { + showErr(err); + throw rte("Unable to start output on DAQ"); + } + if (status != SS_RUNNING) { + throw rte("Unable to start output on DAQ"); + } + totalFramesCount = transferStatus.currentTotalCount; + topenqueued = true; + botenqueued = true; +} +bool OutBufHandler::operator()() { + + /* DEBUGTRACE_ENTER; */ + bool res = true; + assert(daq.getHandle() != 0); + + UlError err = ERR_NO_ERROR; + + ScanStatus status; + TransferStatus transferStatus; + + err = ulAOutScanStatus(daq.getHandle(), &status, &transferStatus); + throwUlException(err); + if (status != SS_RUNNING) { + return false; + } + us increment = transferStatus.currentTotalCount - totalFramesCount; + totalFramesCount += increment; + + if (increment > nFramesPerBlock) { + throw Daq::StreamException(Daq::StreamStatus::StreamError::outputXRun); + } + + if (transferStatus.currentIndex < buffer_mid_idx) { + topenqueued = false; + if (!botenqueued) { + DaqData d(nFramesPerBlock, 1, dtype_descr.dtype); + res = cb(d); + d.copyToRaw(0, reinterpret_cast(&(buf[buffer_mid_idx]))); + + botenqueued = true; + } + } else { + botenqueued = false; + if (!topenqueued) { + DaqData d(nFramesPerBlock, 1, dtype_descr.dtype); + res = cb(d); + d.copyToRaw(0, reinterpret_cast(&(buf[0]))); + + topenqueued = true; + } + } + return res; +} + +OutBufHandler::~OutBufHandler() { + DEBUGTRACE_ENTER; + UlError err = ulAOutScanStop(daq.getHandle()); + if (err != ERR_NO_ERROR) { + showErr(err); + } +} + +void DT9837A::threadFcn(InDaqCallback inCallback, OutDaqCallback outCallback) { + + DEBUGTRACE_ENTER; + + try { + + std::unique_ptr obh; + std::unique_ptr ibh; + + StreamStatus status = _streamStatus; + status.isRunning = true; + _streamStatus = status; + + if (nenoutchannels() > 0) { + assert(outCallback); + obh = std::make_unique(*this, outCallback); + } + if (neninchannels() > 0) { + assert(inCallback); + ibh = std::make_unique(*this, inCallback); + } + if (obh) + obh->start(); + if (ibh) + ibh->start(); + + const double sleeptime_s = + static_cast(_nFramesPerBlock) / (16 * samplerate()); + const us sleeptime_us = static_cast(sleeptime_s * 1e6); + + while (!_stopThread) { + if (ibh) { + if (!(*ibh)()) { + _stopThread = true; + break; + } + } + if (obh) { + if (!(*obh)()) { + _stopThread = true; + break; + } + } else { + std::this_thread::sleep_for(std::chrono::microseconds(sleeptime_us)); + } + } + + /// Update stream status that we are not running anymore + status.isRunning = false; + _streamStatus = status; + _stopThread = false; + + } catch (StreamException &e) { + + StreamStatus status = _streamStatus; + // Copy over error type + status.errorType = e.e; + _streamStatus = status; + + /* + cerr << "\n******************\n"; + cerr << "Catched error in UlDAQ thread: " << e.what() << endl; + cerr << "\n******************\n"; + */ + } +} + +#endif // LASP_HAS_ULDAQ diff --git a/src/lasp/device/lasp_uldaq_impl.h b/src/lasp/device/lasp_uldaq_impl.h new file mode 100644 index 0000000..9fce2e5 --- /dev/null +++ b/src/lasp/device/lasp_uldaq_impl.h @@ -0,0 +1,154 @@ +#pragma once +#include "lasp_daq.h" +#include +#include +#include +#include +#include +#include +#include +#include + +using std::atomic; +using std::cerr; +using std::endl; +using rte = std::runtime_error; + +/** + * @brief UlDaq-specific device information. Adds a copy of the underlying + * DaqDeDaqDeviceDescriptor. + */ +class UlDaqDeviceInfo : public DeviceInfo { + +public: + DaqDeviceDescriptor _uldaqDescriptor; + virtual std::unique_ptr clone() const { + return std::make_unique(*this); + } +}; + +class DT9837A : public Daq { + + DaqDeviceHandle _handle = 0; + std::mutex _daqmutex; + + std::thread _thread; + atomic _stopThread{false}; + atomic _streamStatus; + + const us _nFramesPerBlock; + + void threadFcn(InDaqCallback inCallback, OutDaqCallback outcallback); + +public: + DaqDeviceHandle getHandle() const { return _handle; } + /** + * @brief Create a DT9837A instance. + * + * @param devinfo DeviceInfo to connect to + * @param config DaqConfiguration settings + */ + DT9837A(const DeviceInfo &devinfo, const DaqConfiguration &config); + + virtual ~DT9837A(); + + bool isRunning() const; + + void stop() override final; + + friend class InBufHandler; + friend class OutBufHandler; + + virtual void start(InDaqCallback inCallback, + OutDaqCallback outCallback) override final; + + virtual StreamStatus getStreamStatus() const override { + return _streamStatus; + } +}; + +/** + * @brief Helper class for managing input and output samples of the DAQ device. + */ +class BufHandler { +protected: + /** + * @brief Reference to underlying Daq + */ + DT9837A &daq; + /** + * @brief The type of data, in this case always double precision floats + */ + const DataTypeDescriptor dtype_descr = dtype_desc_fl64; + /** + * @brief The number of channels, number of frames per callback (block). + */ + us nchannels, nFramesPerBlock; + /** + * @brief Sampling frequency in Hz + */ + double samplerate; + std::vector buf; + /** + * @brief Whether the top / bottom part of the buffer are ready to be + * enqueued + */ + bool topenqueued, botenqueued; + + /** + * @brief Counter for the total number of frames acquired / sent since the + * start of the stream. + */ + us totalFramesCount = 0; + long long buffer_mid_idx; + +public: + /** + * @brief Initialize bufhandler + * + * @param daq + * @param nchannels + */ + BufHandler(DT9837A &daq, const us nchannels) + : daq(daq), dtype_descr(daq.dtypeDescr()), nchannels(nchannels), + nFramesPerBlock(daq.framesPerBlock()), samplerate(daq.samplerate()), + buf(2 * nchannels * + nFramesPerBlock, // Watch the two here, the top and the bottom! + 0), + buffer_mid_idx(nchannels * nFramesPerBlock) { + assert(nchannels > 0); + } +}; +/** + * @brief Specific helper for the input buffer + */ +class InBufHandler : public BufHandler { + bool monitorOutput; + InDaqCallback cb; + +public: + InBufHandler(DT9837A &daq, InDaqCallback cb); + void start(); + /** + * @brief InBufHandler::operator()() + * + * @return true on success + */ + bool operator()(); + ~InBufHandler(); +}; +class OutBufHandler : public BufHandler { + OutDaqCallback cb; + +public: + OutBufHandler(DT9837A &daq, OutDaqCallback cb); + void start(); + /** + * @brief OutBufHandler::operator() + * + * @return true on success + */ + bool operator()(); + + ~OutBufHandler(); +}; diff --git a/src/lasp/dsp/CMakeLists.txt b/src/lasp/dsp/CMakeLists.txt index 61a8f59..c3332d6 100644 --- a/src/lasp/dsp/CMakeLists.txt +++ b/src/lasp/dsp/CMakeLists.txt @@ -7,6 +7,7 @@ set(lasp_dsp_files lasp_window.cpp lasp_fft.cpp lasp_rtaps.cpp + lasp_rtsignalviewer.cpp lasp_avpowerspectra.cpp lasp_biquadbank.cpp lasp_thread.cpp @@ -14,6 +15,7 @@ set(lasp_dsp_files lasp_slm.cpp lasp_threadedindatahandler.cpp lasp_ppm.cpp + lasp_clip.cpp ) diff --git a/src/lasp/dsp/lasp_biquadbank.cpp b/src/lasp/dsp/lasp_biquadbank.cpp index 115f6e3..09ec834 100644 --- a/src/lasp/dsp/lasp_biquadbank.cpp +++ b/src/lasp/dsp/lasp_biquadbank.cpp @@ -41,6 +41,45 @@ SeriesBiquad::SeriesBiquad(const vd &filter_coefs) { "Filter coefficients should have fourth element (a0) equal to 1.0"); } } + +SeriesBiquad SeriesBiquad::firstOrderHighPass(const d fs, const d cuton_Hz) { + + if(fs <= 0) { + throw rte("Invalid sampling frequency: " + std::to_string(fs) + " [Hz]"); + } + if(cuton_Hz <= 0) { + throw rte("Invalid cuton frequency: " + std::to_string(cuton_Hz) + " [Hz]"); + } + if(cuton_Hz >= 0.98*fs/2) { + throw rte("Invalid cuton frequency. We limit this to 0.98* fs / 2. Given value" + std::to_string(cuton_Hz) + " [Hz]"); + } + + const d tau = 1/(2*arma::datum::pi*cuton_Hz); + const d facnum = 2*fs*tau/(1+2*fs*tau); + const d facden = (1-2*fs*tau)/(1+2*fs*tau); + + vd coefs(6); + // b0 + coefs(0) = facnum; + // b1 + coefs(1) = -facnum; + // b2 + coefs(2) = 0; + + // a0 + coefs(3) = 1; + + // a1 + coefs(4) = facden; + + // a2 + coefs(5) = 0; + + return SeriesBiquad(coefs); + +} + + std::unique_ptr SeriesBiquad::clone() const { // sos.as_col() concatenates all columns, exactly what we want. return std::make_unique(sos.as_col()); diff --git a/src/lasp/dsp/lasp_biquadbank.h b/src/lasp/dsp/lasp_biquadbank.h index c37df12..1426cb5 100644 --- a/src/lasp/dsp/lasp_biquadbank.h +++ b/src/lasp/dsp/lasp_biquadbank.h @@ -39,6 +39,17 @@ public: virtual ~SeriesBiquad() override {} void reset() override final; std::unique_ptr clone() const override final; + + /** + * @brief Create a SeriesBiquad object for a first order high-pass filter + * + * @param fs Sampling frequency [Hz] + * @param cuton_Hz Cuton-frequency [Hz] + * + * @return SeriesBiquad object with a single biquad, corresponding to a first + * order high pass filter. + */ + static SeriesBiquad firstOrderHighPass(const d fs, const d cuton_Hz); }; /** diff --git a/src/lasp/dsp/lasp_clip.cpp b/src/lasp/dsp/lasp_clip.cpp new file mode 100644 index 0000000..a3b8c8f --- /dev/null +++ b/src/lasp/dsp/lasp_clip.cpp @@ -0,0 +1,94 @@ +/* #define DEBUGTRACE_ENABLED */ +#include "debugtrace.hpp" +#include "lasp_clip.h" +#include + +using std::cerr; +using std::endl; + +using Lck = std::scoped_lock; +using rte = std::runtime_error; + +ClipHandler::ClipHandler(StreamMgr &mgr) + : ThreadedInDataHandler(mgr){ + + DEBUGTRACE_ENTER; + } + +bool ClipHandler::inCallback_threaded(const DaqData &d) { + + DEBUGTRACE_ENTER; + Lck lck(_mtx); + + dmat data = d.toFloat(); + + const us nchannels = d.nchannels; + assert(data.n_cols == nchannels); + + if (nchannels != _clip_time.size()) { + DEBUGTRACE_PRINT("Resizing clip indication"); + _clip_time = vd(nchannels, arma::fill::value(-1)); + } + + /// Obtain max abs values for each column, convert this row-vec to a column + /// vector, and scale with the max-range for each channel + vd maxabs = arma::max(arma::abs(data), 0).as_col() / _max_range; + + /// Find indices for channels that have a clip + arma::uvec clips = maxabs > clip_point; + + for (us i = 0; i < nchannels; i++) { + if (clips(i)) { + /// Reset clip counter + _clip_time(i) = 0; + } else if (_clip_time(i) > clip_indication_time) { + /// Reset to 'unclipped' + _clip_time(i) = -1; + } else if (_clip_time(i) >= 0) { + /// Add a bit of clip time + _clip_time(i) += _dt; + } + } + return true; +} + +arma::uvec ClipHandler::getCurrentValue() const { + + DEBUGTRACE_ENTER; + Lck lck(_mtx); + + arma::uvec clips(_clip_time.size(), arma::fill::zeros); + clips.elem(arma::find(_clip_time >= 0)).fill(1); + + return {clips}; +} + +void ClipHandler::reset(const Daq *daq) { + + DEBUGTRACE_ENTER; + Lck lck(_mtx); + + if (daq) { + + const us nchannels = daq->neninchannels(); + _max_range.resize(nchannels); + + dvec ranges = daq->inputRangeForEnabledChannels(); + assert(ranges.size() == nchannels); + for(us i=0;ineninchannels();i++) { + _max_range[i] = ranges[i]; + } + + _clip_time.fill(-1); + + const d fs = daq->samplerate(); + /* DEBUGTRACE_PRINT(fs); */ + _dt = daq->framesPerBlock() / fs; + } +} + +ClipHandler::~ClipHandler() { + DEBUGTRACE_ENTER; + Lck lck(_mtx); + stop(); +} diff --git a/src/lasp/dsp/lasp_clip.h b/src/lasp/dsp/lasp_clip.h new file mode 100644 index 0000000..e5f76e5 --- /dev/null +++ b/src/lasp/dsp/lasp_clip.h @@ -0,0 +1,77 @@ +// lasp_clip.h +// +// Author: J.A. de Jong, C.R.D. Jansen - ASCEE +// +// Description: Clip handler +#pragma once +#include +#include "lasp_filter.h" +#include "lasp_mathtypes.h" +#include "lasp_threadedindatahandler.h" + +/** + * \addtogroup dsp + * @{ + * + * \addtogroup rt + * @{ + */ + + +/** + * @brief Clipping detector (Clip). Detects when a signal overdrives the input + * */ +class ClipHandler: public ThreadedInDataHandler { + + /** + * @brief Assuming full scale of a signal is +/- 1.0. If a value is found + */ + static inline const d clip_point = 0.98; + + /** + * @brief How long it takes in [s] after a clip event has happened, that we + * are actually still in 'clip' mode. + */ + static inline const d clip_indication_time = 2.0; + + /** + * @brief Inverse of block sampling frequency [s]: (framesPerBlock/fs) + */ + d _dt; + + + mutable std::mutex _mtx; + /** + * @brief How long ago the last clip has happened. Negative in case no clip + * has happened. + */ + vd _clip_time; + + /** + * @brief Storage for maximum values + */ + vd _max_range; + + public: + /** + * @brief Constructs Clipping indicator + * + * @param mgr Stream Mgr to operate on + */ + ClipHandler(StreamMgr& mgr); + ~ClipHandler(); + + /** + * @brief Get the current values of the clip handler. Returns a vector of of True's. + * + * @return clipping indication + */ + arma::uvec getCurrentValue() const; + + bool inCallback_threaded(const DaqData& ) override final; + void reset(const Daq*) override final; + +}; + +/** @} */ +/** @} */ diff --git a/src/lasp/dsp/lasp_fft.cpp b/src/lasp/dsp/lasp_fft.cpp index 829a26e..f5c556a 100644 --- a/src/lasp/dsp/lasp_fft.cpp +++ b/src/lasp/dsp/lasp_fft.cpp @@ -135,7 +135,7 @@ cmat Fft::fft(const dmat &freqdata) { /// * WARNING *. This was source of a serious bug. It is not possible to run /// FFT's and IFFT's on the same _impl, as it overwrites the same memory. /// Uncommenting the line below results in faulty results. - /// #pragma omp parallel for + // #pragma omp parallel for for (us colno = 0; colno < freqdata.n_cols; colno++) { res.col(colno) = _impl->fft(freqdata.col(colno)); } @@ -152,7 +152,7 @@ dmat Fft::ifft(const cmat &freqdata) { /// * WARNING *. This was source of a serious bug. It is not possible to run /// FFT's and IFFT's on the same _impl, as it overwrites the same memory. /// Uncommenting the line below results in faulty results. - /// #pragma omp parallel for + // #pragma omp parallel for for (us colno = 0; colno < freqdata.n_cols; colno++) { res.col(colno) = _impl->ifft(freqdata.col(colno)); diff --git a/src/lasp/dsp/lasp_ppm.h b/src/lasp/dsp/lasp_ppm.h index 0499f83..51b0872 100644 --- a/src/lasp/dsp/lasp_ppm.h +++ b/src/lasp/dsp/lasp_ppm.h @@ -84,7 +84,14 @@ class PPMHandler: public ThreadedInDataHandler { */ std::tuple getCurrentValue() const; - bool inCallback_threaded(const DaqData& ) override final; + /** + * @brief Implements callback on thread + * + * @param d DaqData to work with + * + * @return true when stream should continue. + */ + bool inCallback_threaded(const DaqData& d) override final; void reset(const Daq*) override final; }; diff --git a/src/lasp/dsp/lasp_rtaps.h b/src/lasp/dsp/lasp_rtaps.h index 5e95a2a..5c86bb0 100644 --- a/src/lasp/dsp/lasp_rtaps.h +++ b/src/lasp/dsp/lasp_rtaps.h @@ -19,6 +19,10 @@ * @{ */ +/** + * @brief Real time spectral estimator using Welch method of spectral + * estimation. + */ class RtAps : public ThreadedInDataHandler { std::unique_ptr _filterPrototype; @@ -42,7 +46,8 @@ public: * @param mgr StreamMgr singleton reference * @param freqWeightingFilter Optionally: the frequency weighting filter. * Nullptr should be given for Z-weighting. - * @param For all other arguments, see constructor of AvPowerSpectra + * + * For all other arguments, see constructor of AvPowerSpectra */ RtAps(StreamMgr &mgr, const Filter *freqWeightingFilter, const us nfft = 2048, const Window::WindowType w = Window::WindowType::Hann, @@ -57,7 +62,14 @@ public: */ ccube getCurrentValue() const; - bool inCallback_threaded(const DaqData &) override final; + /** + * @brief Implements the work to to when new DaqData arrives + * + * @param d DaqData to use for computing/updating spectra + * + * @return true if stream should continue. + */ + bool inCallback_threaded(const DaqData & d) override final; void reset(const Daq *) override final; }; diff --git a/src/lasp/dsp/lasp_rtsignalviewer.cpp b/src/lasp/dsp/lasp_rtsignalviewer.cpp new file mode 100644 index 0000000..fa16ca6 --- /dev/null +++ b/src/lasp/dsp/lasp_rtsignalviewer.cpp @@ -0,0 +1,104 @@ +/* #define DEBUGTRACE_ENABLED */ +#include "debugtrace.hpp" +#include "lasp_rtsignalviewer.h" +#include +#include + +using std::cerr; +using std::endl; +using Lck = std::scoped_lock; +using rte = std::runtime_error; + +RtSignalViewer::RtSignalViewer(StreamMgr &mgr, const d approx_time_hist, + const us resolution, const us channel) + : ThreadedInDataHandler(mgr), _approx_time_hist(approx_time_hist), + _resolution(resolution), _channel(channel) { + + DEBUGTRACE_ENTER; + + if (approx_time_hist <= 0) { + throw rte("Invalid time history. Should be > 0"); + } + if (resolution <= 1) { + throw rte("Invalid resolution. Should be > 1"); + } + } + +bool RtSignalViewer::inCallback_threaded(const DaqData &data) { + + DEBUGTRACE_ENTER; + + Lck lck(_sv_mtx); + + /// Push new data in time buffer, scaled by sensitivity + _tb.push(data.toFloat(_channel) / _sens(_channel)); + + _dat.col(0) += _nsamples_per_point / _fs; + /// Required number of samples for a new 'point' + while (_tb.size() > _nsamples_per_point) { + // Roll forward time column + _dat.col(0) += _nsamples_per_point / _fs; + vd newvals = _tb.pop(_nsamples_per_point); + d newmax = newvals.max(); + d newmin = newvals.min(); + + vd mincol = _dat.col(1); + vd maxcol = _dat.col(2); + _dat(arma::span(0, _resolution-2),1) = mincol(arma::span(1, _resolution-1)); + _dat(arma::span(0, _resolution-2),2) = maxcol(arma::span(1, _resolution-1)); + _dat(_resolution-1, 1) = newmin; + _dat(_resolution-1, 2) = newmax; + } + + return true; +} + +RtSignalViewer::~RtSignalViewer() { + Lck lck(_sv_mtx); + stop(); +} +void RtSignalViewer::reset(const Daq *daq) { + + DEBUGTRACE_ENTER; + Lck lck(_sv_mtx); + + // Reset time buffer + _tb.reset(); + + if (daq) { + _fs = daq->samplerate(); + + /// Update sensitivity values + _sens.resize(daq->neninchannels()); + us i = 0; + for (const auto &ch : daq->enabledInChannels()) { + _sens[i] = ch.sensitivity; + i++; + } + + /// The number of samples per point, based on which minmax are computed + _nsamples_per_point = + std::max(static_cast(_approx_time_hist / _resolution * _fs), 1); + + const d dt_points = _nsamples_per_point / _fs; + const d tend = dt_points * _resolution; + + // Initialize data array + _dat.resize(_resolution, 3); + _dat.col(0) = arma::linspace(0, tend, _resolution); + _dat.col(1).zeros(); + _dat.col(2).zeros(); + + } else { + _sens.zeros(); + _fs = 0; + _dat.reset(); + } +} + +dmat RtSignalViewer::getCurrentValue() const { + + DEBUGTRACE_ENTER; + Lck lck(_sv_mtx); + return _dat; +} diff --git a/src/lasp/dsp/lasp_rtsignalviewer.h b/src/lasp/dsp/lasp_rtsignalviewer.h new file mode 100644 index 0000000..85fefd1 --- /dev/null +++ b/src/lasp/dsp/lasp_rtsignalviewer.h @@ -0,0 +1,93 @@ +// lasp_threadedaps.h +// +// Author: J.A. de Jong - ASCEE +// +// Description: Real Time Signal Viewer. +#pragma once +#include "lasp_avpowerspectra.h" +#include "lasp_filter.h" +#include "lasp_mathtypes.h" +#include "lasp_threadedindatahandler.h" +#include "lasp_timebuffer.h" +#include +#include + +/** + * \addtogroup dsp + * @{ + * + * \addtogroup rt + * @{ + */ + +/** + * @brief Real time signal viewer. Shows envelope of the signal based on amount + * of history shown. + */ +class RtSignalViewer : public ThreadedInDataHandler { + + /** + * @brief Storage for sensitivity values + */ + vd _sens; + + /** + * @brief Storage for time samples + */ + TimeBuffer _tb; + + /** + * @brief The amount of time history to show + */ + const d _approx_time_hist; + /** + * @brief The number of points to show + */ + const us _resolution; + + /** + * @brief The channel to show + */ + const us _channel; + + us _nsamples_per_point; + d _fs; + + dmat _dat; + + /** + * @brief Mutex only for _signalviewer. + */ + mutable std::mutex _sv_mtx; + +public: + /** + * @brief + * + * @param mgr Stream manager + * @param approx_time_hist The *approximate* time history to show. The exact + * time history will be calculated such that there fits an integer number of + * samples in a single 'point'. + * @param resolution Number of time points + * @param channel The channel number + */ + RtSignalViewer(StreamMgr &mgr, const d approx_time_hist, const us resolution, + const us channel); + + ~RtSignalViewer(); + + /** + * @brief Returns a 2D array, with: + * - first column: time instances. + * - second column: minimum value of signal + * - third column: maximum value of signal + * + */ + dmat getCurrentValue() const; + + bool inCallback_threaded(const DaqData &) override final; + void reset(const Daq *) override final; +}; + +/** @} */ +/** @} */ diff --git a/src/lasp/dsp/lasp_siggen_impl.cpp b/src/lasp/dsp/lasp_siggen_impl.cpp index 0900418..8aa72a6 100644 --- a/src/lasp/dsp/lasp_siggen_impl.cpp +++ b/src/lasp/dsp/lasp_siggen_impl.cpp @@ -6,10 +6,10 @@ // Signal generators implementation ////////////////////////////////////////////////////////////////////// /* #define DEBUGTRACE_ENABLED */ -#include "debugtrace.hpp" #include "lasp_siggen_impl.h" #include "debugtrace.hpp" #include "lasp_mathtypes.h" +#include using rte = std::runtime_error; @@ -26,13 +26,13 @@ vd Noise::genSignalUnscaled(us nframes) { } void Noise::resetImpl() {} -Sine::Sine(const d freq) : omg(2 * arma::datum::pi * freq) { DEBUGTRACE_ENTER;} +Sine::Sine(const d freq) : omg(2 * arma::datum::pi * freq) { DEBUGTRACE_ENTER; } vd Sine::genSignalUnscaled(const us nframes) { /* DEBUGTRACE_ENTER; */ const d pi = arma::datum::pi; vd phase_vec = - arma::linspace(phase, phase + omg * (nframes - 1) / _fs, nframes); + arma::linspace(phase, phase + omg * (nframes - 1) / _fs, nframes); phase += omg * nframes / _fs; while (phase > 2 * arma::datum::pi) { phase -= 2 * pi; @@ -43,10 +43,10 @@ vd Sine::genSignalUnscaled(const us nframes) { vd Periodic::genSignalUnscaled(const us nframes) { vd res(nframes); - if(_signal.size() == 0) { + if (_signal.size() == 0) { throw rte("No signal defined while calling"); } - for(us i=0;i pre_filter, - std::vector> bandpass) - : _pre_filter(std::move(pre_filter)), _bandpass(std::move(bandpass)), - _alpha(exp(-1 / (fs * tau))), - _sp_storage(_bandpass.size(), arma::fill::zeros), // Storage for - // components of - // single pole low pass - // filter - Lrefsq(Lref*Lref), // Reference level - downsampling_fac(downsampling_fac), + std::unique_ptr pre_filter, + std::vector> bandpass) + : _pre_filter(std::move(pre_filter)), _bandpass(std::move(bandpass)), + _alpha(exp(-1 / (fs * tau))), + _sp_storage(_bandpass.size(), arma::fill::zeros), // Storage for + // components of + // single pole low pass + // filter + Lrefsq(Lref * Lref), // Reference level + downsampling_fac(downsampling_fac), - // Initalize mean square - Pm(_bandpass.size(), arma::fill::zeros), + // Initalize mean square + Pm(_bandpass.size(), arma::fill::zeros), - // Initalize max - Pmax(_bandpass.size(), arma::fill::zeros), + // Initalize max + Pmax(_bandpass.size(), arma::fill::zeros), - // Initalize peak - Ppeak(_bandpass.size(), arma::fill::zeros) + // Initalize peak + Ppeak(_bandpass.size(), arma::fill::zeros) { DEBUGTRACE_ENTER; + DEBUGTRACE_PRINT(_alpha); // Make sure thread pool is running getPool(); @@ -42,7 +43,7 @@ SLM::SLM(const d fs, const d Lref, const us downsampling_fac, const d tau, if (Lref <= 0) { throw rte("Invalid reference level"); } - if (tau <= 0) { + if (tau < 0) { throw rte("Invalid time constant for Single pole lowpass filter"); } if (fs <= 0) { @@ -66,43 +67,43 @@ std::vector> createBandPass(const dmat &coefs) { } return bf; } -us SLM::suggestedDownSamplingFac(const d fs,const d tau) { - if(tau<0) throw rte("Invalid time weighting time constant"); - if(fs<=0) throw rte("Invalid sampling frequency"); +us SLM::suggestedDownSamplingFac(const d fs, const d tau) { + if (fs <= 0) + throw rte("Invalid sampling frequency"); // A reasonable 'framerate' for the sound level meter, based on the // filtering time constant. if (tau > 0) { d fs_slm = 10 / tau; - if(fs_slm < 30) { + if (fs_slm < 30) { fs_slm = 30; } - return std::max((us) 1, static_cast(fs / fs_slm)); + return std::max((us)1, static_cast(fs / fs_slm)); } else { return 1; } } SLM SLM::fromBiquads(const d fs, const d Lref, const us downsampling_fac, - const d tau, const vd &pre_filter_coefs, - const dmat &bandpass_coefs) { + const d tau, const vd &pre_filter_coefs, + const dmat &bandpass_coefs) { DEBUGTRACE_ENTER; return SLM(fs, Lref, downsampling_fac, tau, - std::make_unique(pre_filter_coefs), - createBandPass(bandpass_coefs)); + std::make_unique(pre_filter_coefs), + createBandPass(bandpass_coefs)); } SLM SLM::fromBiquads(const d fs, const d Lref, const us downsampling_fac, - const d tau, const dmat &bandpass_coefs) { + const d tau, const dmat &bandpass_coefs) { DEBUGTRACE_ENTER; return SLM(fs, Lref, downsampling_fac, tau, - nullptr, // Pre-filter - createBandPass(bandpass_coefs) // Bandpass coefficients - ); + nullptr, // Pre-filter + createBandPass(bandpass_coefs) // Bandpass coefficients + ); } -vd SLM::run_single(vd work,const us i) { +vd SLM::run_single(vd work, const us i) { // Filter input in-place _bandpass[i]->filter(work); @@ -126,7 +127,7 @@ vd SLM::run_single(vd work,const us i) { for (us j = 0; j < work.n_rows; j++) { // Update mean square of signal, work is here still signal power Pm(i) = (Pm(i) * static_cast(N_local) + work(j)) / - (static_cast(N_local) + 1); + (static_cast(N_local) + 1); N_local++; @@ -142,7 +143,7 @@ vd SLM::run_single(vd work,const us i) { Pmax(i) = std::max(Pmax(i), arma::max(work)); // Convert to levels in dB - work = 10*arma::log10((work+arma::datum::eps)/Lrefsq); + work = 10 * arma::log10((work + arma::datum::eps) / Lrefsq); return work; } @@ -170,8 +171,8 @@ dmat SLM::run(const vd &input_orig) { /* DEBUGTRACE_PRINT(res.n_rows); */ /* DEBUGTRACE_PRINT(futs.size()); */ - // Update the total number of samples harvested so far. NOTE: *This should be - // done AFTER the threads are done!!!* + // Update the total number of samples harvested so far. NOTE: *This should + // be done AFTER the threads are done!!!* } N += input.n_rows; diff --git a/src/lasp/dsp/lasp_threadedindatahandler.h b/src/lasp/dsp/lasp_threadedindatahandler.h index f628f44..20f5df2 100644 --- a/src/lasp/dsp/lasp_threadedindatahandler.h +++ b/src/lasp/dsp/lasp_threadedindatahandler.h @@ -1,5 +1,4 @@ #pragma once -#include #include "lasp_streammgr.h" const us RINGBUFFER_SIZE = 1024; @@ -52,11 +51,11 @@ class ThreadedInDataHandler: public InDataHandler { * @brief This function should be overridden with an actual implementation, * of what should happen on a different thread. * - * @param DaqData Input daq data + * @param d Input daq data * * @return true on succes. False when an error occured. */ - virtual bool inCallback_threaded(const DaqData&) = 0; + virtual bool inCallback_threaded(const DaqData& d) = 0; }; diff --git a/src/lasp/filter/filterbank_design.py b/src/lasp/filter/filterbank_design.py index 777b720..7ca0903 100644 --- a/src/lasp/filter/filterbank_design.py +++ b/src/lasp/filter/filterbank_design.py @@ -85,8 +85,8 @@ class FilterBankDesigner: if self.nominal_txt(x) == nom_txt: return x raise ValueError( - f'Could not find a nominal frequency corresponding to {nom_txt}. Hint: use \'5k\' instead of \'5000\'.' - ) + f'Could not find a nominal frequency corresponding to {nom_txt}. ' + 'Hint: use \'5k\' instead of \'5000\'.') def sanitize_input(self, input_): if isinstance(input_, int): @@ -172,7 +172,6 @@ class FilterBankDesigner: Args: x: Band designator """ - SOS_ORDER = 5 fs = self.fs @@ -186,6 +185,18 @@ class FilterBankDesigner: x = self.sanitize_input(x) fu_n = fu / fnyq + # Handle fl & fu >= fNyquist: return 'no pass' + if fl_n >= 1: + return np.tile(np.asarray([0, 0, 0, 1, 0, 0]), (SOS_ORDER, 1)) + + # Handle fu >= fNyquist: return high pass + if fu_n >= 1: + highpass = butter(SOS_ORDER, fl_n, output='sos', btype='highpass') + allpass = np.tile(np.asarray([1, 0, 0, 1, 0, 0]), + (SOS_ORDER-highpass.shape[0], 1)) + return np.vstack((highpass, allpass)) # same shape=(SOS_ORDER, 6) + + # Regular bands return butter(SOS_ORDER, [fl_n, fu_n], output='sos', btype='band') def firFreqResponse(self, x, freq): @@ -254,8 +265,8 @@ class FilterBankDesigner: for i, x in enumerate(range(xl, xu + 1)): fl = self.fl(x) fu = self.fu(x) - # Find the indices in the frequency array which correspond to the - # frequency band x + # Find the indices in the frequency array which correspond to + # the frequency band x if x != xu: indices_cur = np.where((freq >= fl) & (freq < fu)) else: @@ -330,20 +341,20 @@ class OctaveBankDesigner(FilterBankDesigner): @property def xs(self): """All possible band designators for an octave band filter.""" - return list(range(-6, 5)) + return list(range(-7, 5)) def band_limits(self, x, filter_class=0): """Returns the octave band filter limits for filter designator x. Args: x: Filter offset power from the reference frequency of 1000 Hz. - filter_class: Either 0 or 1, defines the tolerances on the frequency - response + filter_class: Either 0 or 1, defines the tolerances on the + frequency response Returns: - freq, llim, ulim: Tuple of Numpy arrays containing the frequencies of - the corner points of the filter frequency response limits, lower limits - in *deciBell*, upper limits in *deciBell*, respectively. + freq, llim, ulim: Tuple of Numpy arrays containing the frequencies + of the corner points of the filter frequency response limits, lower + limits in *deciBell*, upper limits in *deciBell*, respectively. """ b = 1 @@ -393,7 +404,8 @@ class OctaveBankDesigner(FilterBankDesigner): -3: '125', -4: '63', -5: '31.5', - -6: '16' + -6: '16', + -7: '8' } assert len(nominals) == len(self.xs) return nominals[x] @@ -450,8 +462,7 @@ class OctaveBankDesigner(FilterBankDesigner): if int(self.fs) in [44100, 48000, 96000]: return 1.0 else: - raise ValueError('Unimplemented sampling frequency for SOS' - 'filter design') + return 1.0 def sosFac_u(self, x): """Right side percentage of change in cut-on frequency for designing @@ -463,8 +474,7 @@ class OctaveBankDesigner(FilterBankDesigner): if int(self.fs) in [44100, 48000, 96000]: return 1.0 else: - raise ValueError('Unimplemented sampling frequency for SOS' - 'filter design') + return 1.0 class ThirdOctaveBankDesigner(FilterBankDesigner): @@ -475,18 +485,16 @@ class ThirdOctaveBankDesigner(FilterBankDesigner): one-third octave bands Args: - fs: Sampling frequency in [Hz] - + fs: Sampling frequency in [Hz] - can be None """ super().__init__(fs) - self.xs = list(range(-16, 14)) + self.xs = list(range(-19, 14)) # Text corresponding to the nominal frequency - self._nominal_txt = [ + self._nominal_txt = ['12.5', '16', '20', '25', '31.5', '40', '50', '63', '80', '100', '125', '160', '200', '250', '315', '400', '500', '630', '800', '1k', '1.25k', '1.6k', '2k', '2.5k', '3.15k', '4k', '5k', '6.3k', '8k', '10k', '12.5k', - '16k', '20k' - ] + '16k', '20k'] assert len(self.xs) == len(self._nominal_txt) @@ -511,13 +519,13 @@ class ThirdOctaveBankDesigner(FilterBankDesigner): Args: x: Filter offset power from the reference frequency of 1000 Hz. - filter_class: Either 0 or 1, defines the tolerances on the frequency - response + filter_class: Either 0 or 1, defines the tolerances on the + frequency response Returns: - freq, llim, ulim: Tuple of Numpy arrays containing the frequencies of - the corner points of the filter frequency response limits, lower limits - in *deciBell*, upper limits in *deciBell*, respectively. + freq, llim, ulim: Tuple of Numpy arrays containing the frequencies + of the corner points of the filter frequency response limits, lower + limits in *deciBell*, upper limits in *deciBell*, respectively. """ fm = self.G**(x / self.b) * self.fr @@ -603,25 +611,15 @@ class ThirdOctaveBankDesigner(FilterBankDesigner): """Left side percentage of change in cut-on frequency for designing the filter.""" # Idea: correct for frequency warping: - if np.isclose(self.fs, 48000): - return 1.00 - elif np.isclose(self.fs, 41000): - warnings.warn( - f'Frequency {self.fs} might not result in correct filters') - - return 1.00 - elif np.isclose(self.fs, 32768): - return 1.00 + if int(self.fs) in [48000]: + return 1.0 else: - raise ValueError('Unimplemented sampling frequency for SOS' - 'filter design') + return 1.0 + def sosFac_u(self, x): """Right side percentage of change in cut-on frequency for designing the filter.""" - if np.isclose(self.fs, 48000): - return 1 - elif np.isclose(self.fs, 32768): - return 1.00 + if int(self.fs) in [48000]: + return 1.0 else: - raise ValueError('Unimplemented sampling frequency for SOS' - 'filter design') + return 1.0 diff --git a/src/lasp/lasp_common.py b/src/lasp/lasp_common.py index 773ea41..696ac76 100644 --- a/src/lasp/lasp_common.py +++ b/src/lasp/lasp_common.py @@ -73,15 +73,21 @@ class Qty: return f'{self.name} [{self.unit_symb}]' def __eq__(self, other): - # logging.debug('eq()') """ Comparison breaks for the other objects, level unit, level ref name, etc as these are tuples / a single string. """ + # logging.debug(f'eq() {self.name} {other.name}') return (self.name == other.name and self.unit_name == other.unit_name) + def toInt(self): + """ + Convert quantity to its specific enum integer value + """ + return self.cpp_enum.value + @unique @@ -147,10 +153,13 @@ class CalSetting: qty: Qty class CalibrationSettings: - one = CalSetting('94 dB SPL', 94.0 , 1.0, SIQtys.AP) - two = CalSetting('114 dB SPL', 114.0 , 10.0, SIQtys.AP) + one = CalSetting('94 dB SPL', 94.0 , 10**(94/20)*2e-5, SIQtys.AP.value) + two = CalSetting('1 Pa rms', 93.98, 1.0, SIQtys.AP.value) + three = CalSetting('114 dB SPL', 114.0 , 10**(114/20)*2e-5, SIQtys.AP.value) + four = CalSetting('10 Pa rms', 113.98, 10.0, SIQtys.AP.value) + five = CalSetting('93.7 dB SPL', 93.7 , 1.0, SIQtys.AP.value) - types = (one, two) + types = (one, two, three, four, five) default = one default_index = 0 @@ -277,46 +286,57 @@ class this_lasp_shelve(Shelve): return os.path.join(lasp_appdir, f'{node}_config.shelve') class TimeWeighting: - none = (-1, 'Raw (no time weighting)') + # This is not actually a time weighting + none = (0, 'Raw (no time weighting)') + uufast = (1e-4, '0.1 ms') ufast = (35e-3, 'Impulse (35 ms)') fast = (0.125, 'Fast (0.125 s)') slow = (1.0, 'Slow (1.0 s)') tens = (10., '10 s') - infinite = (0, 'Infinite') - types_realtime = (ufast, fast, slow, tens, infinite) - types_all = (none, uufast, ufast, fast, slow, tens, infinite) + + # All-averaging: compute the current average of all history. Kind of the + # equivalent level. Only useful for power spectra. For Sound Level Meter, + # equivalent levels does that thing. + averaging = (-1, 'All-averaging') + + types_all = (none, uufast, ufast, fast, slow, tens, averaging) + + # Used for combobox of realt time power spectra + types_realtimeaps = (ufast, fast, slow, tens, averaging) + # Used for combobox of realt time sound level meter + types_realtimeslm = (ufast, fast, slow, tens) + + # Used for combobox of sound level meter - time dependent + types_slmt = (none, uufast, ufast, fast, slow, tens) + + # Used for combobox of sound level meter - Lmax statistics + types_slmstats = (uufast, ufast, fast, slow, tens) default = fast - default_index = 3 - default_index_realtime = 1 @staticmethod - def fillComboBox(cb, realtime=False): + def fillComboBox(cb, types): """ Fill TimeWeightings to a combobox Args: cb: QComboBox to fill + types: The types to fill it with """ cb.clear() - if realtime: - types = TimeWeighting.types_realtime - defindex = TimeWeighting.default_index_realtime - else: - types = TimeWeighting.types_all - defindex = TimeWeighting.default_index + + logging.debug(f'{types}') for tw in types: cb.addItem(tw[1], tw) - - cb.setCurrentIndex(defindex) + if TimeWeighting.default == tw: + cb.setCurrentIndex(cb.count()-1) @staticmethod def getCurrent(cb): - if cb.count() == len(TimeWeighting.types_realtime): - return TimeWeighting.types_realtime[cb.currentIndex()] - else: - return TimeWeighting.types_all[cb.currentIndex()] + for type in TimeWeighting.types_all: + if cb.currentText() == type[1]: + return type class FreqWeighting: diff --git a/src/lasp/lasp_cpp.cpp b/src/lasp/lasp_cpp.cpp index 91ce669..f541f1b 100644 --- a/src/lasp/lasp_cpp.cpp +++ b/src/lasp/lasp_cpp.cpp @@ -8,14 +8,23 @@ * Welcome to the LASP (Library for Acoustic Signal Processing) code * documentation. The code comprises a part which is written in C++, a part * that is written in Python, and a part that functions as glue, which is - * Pybind11 C++ glue code. An example of such a file is the current one. + * Pybind11 C++ glue code. An example of such a file is lasp_cpp.cpp. + * This is the internal documentation of LASP. It serves as background + * information for programmers. * * \section Installation * - * For the installation manual, please refer to the README.md of the Git - * repository. It is recommended to install the software from source. + * For the installation manual, please refer to the README + * of the Git repository. * * + * \section Usage + * + * Some usage examples are given in the examples + * directory of the repository. + * * */ /** @@ -48,7 +57,6 @@ PYBIND11_MODULE(lasp_cpp, m) { m.attr("__version__") = std::to_string(LASP_VERSION_MAJOR) + "." + std::to_string(LASP_VERSION_MINOR); - m.attr("LASP_VERSION_MAJOR") = LASP_VERSION_MAJOR; m.attr("LASP_VERSION_MINOR") = LASP_VERSION_MINOR; } diff --git a/src/lasp/lasp_daqconfigs.py b/src/lasp/lasp_daqconfigs.py index fe80f4f..c7bf381 100644 --- a/src/lasp/lasp_daqconfigs.py +++ b/src/lasp/lasp_daqconfigs.py @@ -29,7 +29,7 @@ class DaqConfigurations: output_config: DaqConfiguration, ): """ - Initialize set of DaqConfigurations. + Initialize set of DaqConfigurations. Args: duplex_mode: If true, the input configuration is used for output as @@ -92,6 +92,18 @@ class DaqConfigurations: output_config = DaqConfiguration.fromTOML(config_str[2]) return DaqConfigurations(duplex_mode, input_config, output_config) + @staticmethod + def loadRaw(): + """ + Returns configurations presets in the raw form they are stored. + + Returns: + all configurations, raw + """ + with lasp_shelve() as sh: + configs_raw = sh.load(f"daqconfigs_v{LASP_VERSION_MAJOR}", {}) + return configs_raw + def save(self, name: str): """ Save the current set of configurations to the shelve store. @@ -109,6 +121,17 @@ class DaqConfigurations: configs_str[name] = [self.duplex_mode, input_str, output_str] sh.store(f"daqconfigs_v{LASP_VERSION_MAJOR}", configs_str) + @staticmethod + def saveRaw(configs_raw): + """ + Save configurations presets, using already formatted data + + Arg: + all configurations, raw data format in form they are stored + """ + with lasp_shelve() as sh: + sh.store(f"daqconfigs_v{LASP_VERSION_MAJOR}", configs_raw) + @staticmethod def delete(name: str): """ diff --git a/src/lasp/lasp_measurement.py b/src/lasp/lasp_measurement.py index 1b1607c..ad0e8e9 100644 --- a/src/lasp/lasp_measurement.py +++ b/src/lasp/lasp_measurement.py @@ -50,7 +50,8 @@ from .lasp_config import LASP_NUMPY_FLOAT_TYPE from scipy.io import wavfile import os, time, wave, logging from .lasp_common import SIQtys, Qty, getFreq -from .lasp_cpp import Window, DaqChannel, LASP_VERSION_MAJOR +from .lasp_cpp import Window, DaqChannel, LASP_VERSION_MAJOR, AvPowerSpectra +from typing import List def getSampWidth(dtype): @@ -195,7 +196,7 @@ class Measurement: # Open the h5 file in read-plus mode, to allow for changing the # measurement comment. - with h5.File(fn, 'r+') as f: + with h5.File(fn, 'r') as f: # Check for video data try: f['video'] @@ -233,10 +234,9 @@ class Measurement: self._channelNames = [f'Unnamed {i}' for i in range(self.nchannels)] # comment = read-write thing - try: + if 'comment' in f.keys(): self._comment = f.attrs['comment'] - except KeyError: - f.attrs['comment'] = '' + else: self._comment = '' # Sensitivity @@ -250,28 +250,27 @@ class Measurement: self._time = f.attrs['time'] + # Quantity stored as channel. self._qtys = None - try: - qtys_json = f.attrs['qtys'] - # Load quantity data - self._qtys = [Qty.from_json(qty_json) for qty_json in qtys_json] - except KeyError: - # If quantity data is not available, this is an 'old' - # measurement file. - pass try: qtys_enum_idx = f.attrs['qtys_enum_idx'] self._qtys = [SIQtys.fromInt(idx) for idx in qtys_enum_idx] except KeyError: - pass + try: + qtys_json = f.attrs['qtys'] + # Load quantity data + self._qtys = [Qty.from_json(qty_json) for qty_json in qtys_json] + except KeyError: + # If quantity data is not available, this is an 'old' + # measurement file. + pass if self._qtys is None: self._qtys = [SIQtys.default() for i in range(self.nchannels)] logging.debug(f'Physical quantity data not available in measurement file. Assuming {SIQtys.default}') - - logging.debug(f'Physical quantity data not available in measurement file. Assuming {SIQtys.default}') + def setAttribute(self, atrname, value): """ Set an attribute in the measurement file, and keep a local copy in @@ -311,14 +310,14 @@ class Measurement: return chcfg @channelConfig.setter - def channelConfig(self, chcfg): + def channelConfig(self, chcfg: List[DaqChannel]): chname = [] sens = [] qtys = [] for ch in chcfg: chname.append(ch.name) sens.append(ch.sensitivity) - qtys.append(SIQtys.fromInt(ch.qty)) + qtys.append(SIQtys.fromCppEnum(ch.qty)) self.channelNames = chname self.sensitivity = sens @@ -332,11 +331,14 @@ class Measurement: def qtys(self, newqtys): if not len(newqtys) == len(self._qtys): raise ValueError('Invalid number of quantities') - qtys_json = [qty.to_json() for qty in newqtys] + qtys_int = [qty.toInt() for qty in newqtys] # Use setAttribute here, but thos store the jsonified version as well, # which we have to overwrite again with the deserialized ones. This is # actually not a very nice way of coding. - self.setAttribute('qtys', qtys_json) + with self.file('r+') as f: + # Update comment attribute in the file + f.attrs['qtys_enum_idx'] = qtys_int + self._qtys = newqtys @contextmanager @@ -477,20 +479,20 @@ class Measurement: """ nfft = kwargs.pop('nfft', 2048) - window = kwargs.pop('window', Window.hann) + window = kwargs.pop('windowType', Window.WindowType.Hann) overlap = kwargs.pop('overlap', 50.0) if channels is None: channels = list(range(self.nchannels)) nchannels = len(channels) - aps = AvPowerSpectra(nfft, nchannels, overlap, window) + aps = AvPowerSpectra(nfft, window, overlap) freq = getFreq(self.samplerate, nfft) for data in self.iterData(channels, **kwargs): - CS = aps.addTimeData(data) + CS = aps.compute(data) - return freq, CS + return freq, aps.get_est() def periodicAverage(self, N, channels=None, noiseCorrection=True, **kwargs): """ diff --git a/src/lasp/lasp_measurementset.py b/src/lasp/lasp_measurementset.py new file mode 100644 index 0000000..6010b7a --- /dev/null +++ b/src/lasp/lasp_measurementset.py @@ -0,0 +1,66 @@ +""" +Provides class MeasurementSet, a class used to perform checks and adjustments +on a group of measurements at the same time. + +""" +__all__ = ['MeasurementSet'] +from .lasp_measurement import Measurement +from typing import List + + +class MeasurementSet(list): + """ + Group of measurements that have some correspondence to one another. Class + is used to operate on multiple measurements at once. + + """ + def __init__(self, mlist: List[Measurement] =[]): + """ + Initialize a measurement set + + Args: + mlist: Measurement list + + """ + if any([not isinstance(i, Measurement) for i in mlist]): + raise TypeError('Object in list should be of Measurement type') + + super().__init__(mlist) + + def measTimeSame(self): + """ + Returns True if all measurements have the same measurement + time (recorded time) + + """ + if len(self) > 0: + first = self[0].N + return all([first == meas.N for meas in self]) + else: + return False + + def measSimilar(self): + """ + Similar means: channel metadata is the same, and the measurement time + is the same. It means that the recorded data is, of course, different. + + Returns: + True if measChannelsSame() and measTimeSame() else False + + """ + return self.measTimeSame() and self.measChannelsSame() + + + def measChannelsSame(self): + """ + This method is used to check whether a set of measurements can be + accessed in a loop, i.e. for computing power spectra or sound levels on + a set of measurements, simultaneously. If the channel data is the same + (name, sensitivity, ...) it returns True. + """ + if len(self) > 0: + first = self[0].channelConfig + return all([first == meas.channelConfig for meas in self]) + else: + return False + diff --git a/src/lasp/lasp_octavefilter.py b/src/lasp/lasp_octavefilter.py index 632cd78..2a183ac 100644 --- a/src/lasp/lasp_octavefilter.py +++ b/src/lasp/lasp_octavefilter.py @@ -266,6 +266,8 @@ class SosFilterBank: return {} filtered_data = self._fb.filter(data) + if filtered_data.ndim == 1: + filtered_data = filtered_data[:, None] # Output given as a dictionary with nom_txt as the key output = {} diff --git a/src/lasp/lasp_record.py b/src/lasp/lasp_record.py index 8f83ae5..627c1ef 100644 --- a/src/lasp/lasp_record.py +++ b/src/lasp/lasp_record.py @@ -135,7 +135,8 @@ class Recording: f.attrs["channelNames"] = [ch.name for ch in in_ch] # Add the start delay here, as firstFrames() is called right after the - # constructor is called. + # constructor is called. time.time() returns a floating point + # number of seconds after epoch. f.attrs["time"] = time.time() + self.startDelay # In V2, we do not store JSON metadata anymore, but just an enumeration @@ -143,6 +144,7 @@ class Recording: f.attrs["qtys_enum_idx"] = [ch.qty.value for ch in in_ch] # Measured physical quantity metadata + # This was how it was in LASP version < 1.0 # f.attrs['qtys'] = [ch.qty.to_json() for ch in in_ch] def firstFrames(self, adata): @@ -174,18 +176,15 @@ class Recording: def inCallback(self, adata): """ - This method should be called to grab data from the input queue, which - is filled by the stream, and put it into a file. It should be called at - a regular interval to prevent overflowing of the queue. It is called - within the start() method of the recording, if block is set to True. - Otherwise, it should be called from its parent at regular intervals. - For example, in Qt this can be done using a QTimer. + This method is called when a block of audio data from the stream is + available. It should return either True or False. + When returning False, it will stop the stream. """ if self.stop(): # Stop flag is raised. We do not add any data anymore. - return + return True with self.file_mtx: diff --git a/src/lasp/lasp_slm.py b/src/lasp/lasp_slm.py index f00af31..9b605dc 100644 --- a/src/lasp/lasp_slm.py +++ b/src/lasp/lasp_slm.py @@ -27,35 +27,35 @@ class SLM: Multi-channel Sound Level Meter. Input data: time data with a certain sampling frequency. Output: time-weighted (fast/slow) sound pressure levels in dB(A/C/Z). Possibly in octave bands. - """ def __init__(self, fs, - fbdesigner=None, - tw: TimeWeighting =TimeWeighting.fast, - fw: FreqWeighting =FreqWeighting.A, - xmin = None, - xmax = None, + fbdesigner=None, + tw: TimeWeighting = TimeWeighting.fast, + fw: FreqWeighting = FreqWeighting.A, + xmin=None, + xmax=None, include_overall=True, - level_ref_value=P_REF): + level_ref_value=P_REF, + offset_t=0): """ Initialize a sound level meter object. Args: fs: Sampling frequency of input data [Hz] - fbdesigner: FilterBankDesigner to use for creating the + fbdesigner: FilterBankDesigner to use for creating the (fractional) octave bank filters. Set this one to None to only do overalls fs: Sampling frequency [Hz] tw: Time Weighting to apply fw: Frequency weighting to apply - xmin: Filter designator of lowest band. + xmin: Filter designator of lowest band. xmax: Filter designator of highest band include_overall: If true, a non-functioning filter is added which is used to compute the overall level. level_ref_value: Reference value for computing the levels in dB - + offset_t: Offset to be added to output time data [s] """ self.fbdesigner = fbdesigner @@ -63,8 +63,8 @@ class SLM: xmin = fbdesigner.xs[0] elif fbdesigner is None: xmin = 0 - - if xmax is None and self.fbdesigner is not None: + + if xmax is None and self.fbdesigner is not None: xmax = fbdesigner.xs[-1] elif fbdesigner is None: xmax = 0 @@ -72,9 +72,11 @@ class SLM: self.xs = list(range(xmin, xmax + 1)) nfilters = len(self.xs) - if include_overall: nfilters +=1 + if include_overall: + nfilters += 1 self.include_overall = include_overall - + self.offset_t = offset_t + spld = SPLFilterDesigner(fs) if fw == FreqWeighting.A: prefilter = spld.A_Sos_design().flatten() @@ -90,7 +92,7 @@ class SLM: # This is a bit of a hack, as the 5 is hard-encoded here, but should in # fact be coming from somewhere else.. - sos_overall = np.array([1,0,0,1,0,0]*5, dtype=float) + sos_overall = np.array([1, 0, 0, 1, 0, 0]*5, dtype=float) if fbdesigner is not None: assert fbdesigner.fs == fs @@ -104,22 +106,22 @@ class SLM: self.nom_txt.append(fbdesigner.nominal_txt(x)) if include_overall: - # Create a unit impulse response filter, every third index equals + # Create a unit impulse response filter, every third index equals # 1, so b0 = 1 and a0 is 1 (by definition) # a0 = 1, b0 = 1, rest is zero - sos[:,-1] = sos_overall + sos[:, -1] = sos_overall self.nom_txt.append('overall') else: # No filterbank, means we do only compute the overall values. This # means that in case of include_overall, it creates two overall # channels. That would be confusing, so we do not allow it. - sos = sos_overall[:,np.newaxis] + sos = sos_overall[:, np.newaxis] self.nom_txt.append('overall') # Downsampling factor, determine from single pole low pass filter time # constant, such that aliasing is ~ allowed at 20 dB lower value - # and + # and dsfac = cppSLM.suggestedDownSamplingFac(fs, tw[0]) if prefilter is not None: @@ -136,7 +138,6 @@ class SLM: # Initialize counter to 0 self.N = 0 - def run(self, data): """ Add new fresh timedata to the Sound Level Meter @@ -153,14 +154,14 @@ class SLM: Ncur = levels.shape[0] tend = tstart + Ncur / self.fs_slm - t = np.linspace(tstart, tend, Ncur, endpoint=False) + t = np.linspace(tstart, tend, Ncur, endpoint=False) + self.offset_t self.N += Ncur output = {} for i, x in enumerate(self.xs): # '31.5' to '16k' - output[self.nom_txt[i]] = {'t': t, + output[self.nom_txt[i]] = {'t': t, 'data': levels[:, [i]], 'x': x} if self.include_overall and self.fbdesigner is not None: @@ -168,7 +169,6 @@ class SLM: return output - def return_as_dict(self, dat): """ Helper function used to return resulting data in a proper way. diff --git a/src/lasp/pybind11/lasp_daq.cpp b/src/lasp/pybind11/lasp_daq.cpp index fa91280..4192f6a 100644 --- a/src/lasp/pybind11/lasp_daq.cpp +++ b/src/lasp/pybind11/lasp_daq.cpp @@ -25,18 +25,17 @@ void init_daq(py::module &m) { .value("driverError", Daq::StreamStatus::StreamError::driverError) .value("systemError", Daq::StreamStatus::StreamError::systemError) .value("threadError", Daq::StreamStatus::StreamError::threadError) - .value("logicError", Daq::StreamStatus::StreamError::logicError) - .value("apiSpecificError", - Daq::StreamStatus::StreamError::apiSpecificError); + .value("logicError", Daq::StreamStatus::StreamError::logicError); ss.def("errorMsg", &Daq::StreamStatus::errorMsg); /// Daq - daq.def("neninchannels", &Daq::neninchannels, py::arg("include_monitor") = true); + daq.def("neninchannels", &Daq::neninchannels, + py::arg("include_monitor") = true); + daq.def("nenoutchannels", &Daq::nenoutchannels); daq.def("samplerate", &Daq::samplerate); daq.def("dataType", &Daq::dataType); daq.def("framesPerBlock", &Daq::framesPerBlock); daq.def("getStreamStatus", &Daq::getStreamStatus); - } diff --git a/src/lasp/pybind11/lasp_deviceinfo.cpp b/src/lasp/pybind11/lasp_deviceinfo.cpp index 3d0ad4a..ad3b9de 100644 --- a/src/lasp/pybind11/lasp_deviceinfo.cpp +++ b/src/lasp/pybind11/lasp_deviceinfo.cpp @@ -1,4 +1,3 @@ -#include #include #include #include "lasp_deviceinfo.h" @@ -36,5 +35,7 @@ void init_deviceinfo(py::module& m) { devinfo.def_readonly("hasInputACCouplingSwitch", &DeviceInfo::hasInputACCouplingSwitch); + devinfo.def_readonly("physicalOutputQty", &DeviceInfo::physicalOutputQty); + } diff --git a/src/lasp/pybind11/lasp_pyindatahandler.cpp b/src/lasp/pybind11/lasp_pyindatahandler.cpp index 69af06f..d8639b3 100644 --- a/src/lasp/pybind11/lasp_pyindatahandler.cpp +++ b/src/lasp/pybind11/lasp_pyindatahandler.cpp @@ -1,11 +1,13 @@ -#include -#define DEBUGTRACE_ENABLED +/* #define DEBUGTRACE_ENABLED */ #include "arma_npy.h" #include "debugtrace.hpp" #include "lasp_ppm.h" +#include "lasp_clip.h" #include "lasp_rtaps.h" +#include "lasp_rtsignalviewer.h" #include "lasp_streammgr.h" #include "lasp_threadedindatahandler.h" +#include #include #include #include @@ -185,11 +187,21 @@ void init_datahandler(py::module &m) { ppm.def(py::init()); ppm.def("getCurrentValue", [](const PPMHandler &ppm) { - std::tuple tp = ppm.getCurrentValue(); return py::make_tuple(ColToNpy(std::get<0>(tp)), - ColToNpy(std::get<1>(tp))); + ColToNpy(std::get<1>(tp))); + }); + + /// Clip Detector + py::class_ clip(m, "ClipHandler"); + clip.def(py::init()); + + clip.def("getCurrentValue", [](const ClipHandler &clip) { + + arma::uvec cval = clip.getCurrentValue(); + + return ColToNpy(cval); // something goes wrong here }); /// Real time Aps @@ -216,6 +228,20 @@ void init_datahandler(py::module &m) { rtaps.def("getCurrentValue", [](RtAps &rt) { ccube val = rt.getCurrentValue(); - return CubeToNpy(val); + return CubeToNpy(val); + }); + + /// Real time Signal Viewer + /// + py::class_ rtsv(m, "RtSignalViewer"); + rtsv.def(py::init()); + + rtsv.def("getCurrentValue", [](RtSignalViewer &rt) { + dmat val = rt.getCurrentValue(); + return MatToNpy(val); }); } diff --git a/src/lasp/pybind11/lasp_siggen.cpp b/src/lasp/pybind11/lasp_siggen.cpp index 9a6e424..9704889 100644 --- a/src/lasp/pybind11/lasp_siggen.cpp +++ b/src/lasp/pybind11/lasp_siggen.cpp @@ -1,10 +1,6 @@ -#include "lasp_avpowerspectra.h" -#include "lasp_biquadbank.h" -#include "lasp_fft.h" #include "lasp_siggen.h" +#include "arma_npy.h" #include "lasp_siggen_impl.h" -#include "lasp_slm.h" -#include "lasp_window.h" #include #include @@ -24,25 +20,30 @@ namespace py = pybind11; */ void init_siggen(py::module &m) { - /// Siggen + /// Siggen py::class_> siggen(m, "Siggen"); siggen.def("setLevel", &Siggen::setLevel, "Set the level of the signal generator"); siggen.def("setDCOffset", &Siggen::setDCOffset); siggen.def("setMute", &Siggen::setMute); siggen.def("setLevel", &Siggen::setLevel); - siggen.def("setLevel", &Siggen::setLevel, py::arg("newLevel"), py::arg("dB") = true); + siggen.def("setLevel", &Siggen::setLevel, py::arg("newLevel"), + py::arg("dB") = true); + siggen.def("reset", &Siggen::reset); siggen.def("setFilter", &Siggen::setFilter); - py::class_> noise(m, "Noise", siggen); + py::class_> noise(m, "Noise", siggen); noise.def(py::init<>()); py::class_> sine(m, "Sine", siggen); sine.def(py::init()); sine.def("setFreq", &Sine::setFreq); - py::class_> periodic(m, "Periodic", siggen); + py::class_> periodic(m, "Periodic", + siggen); + periodic.def("getSequence", + [](const Sweep &s) { return ColToNpy(s.getSequence()); }); py::class_> sweep(m, "Sweep", periodic); sweep.def(py::init()); @@ -50,6 +51,5 @@ void init_siggen(py::module &m) { sweep.def_readonly_static("BackwardSweep", &Sweep::BackwardSweep); sweep.def_readonly_static("LinearSweep", &Sweep::LinearSweep); sweep.def_readonly_static("LogSweep", &Sweep::LogSweep); - } /** @} */ diff --git a/src/lasp/tools/tools.py b/src/lasp/tools/tools.py index 02c3201..8329f46 100644 --- a/src/lasp/tools/tools.py +++ b/src/lasp/tools/tools.py @@ -79,6 +79,70 @@ class SmoothingType: # TO DO: add possibility to insert data that is not lin spaced in frequency +def smoothCalcMatrix(freq, sw: SmoothingWidth): + """ + Args: + freq: array of frequencies of data points [Hz] - equally spaced + sw: SmoothingWidth + + Returns: + freq: array frequencies of data points [Hz] + Q: matrix to smooth power: {fsm} = [Q] * {fraw} + + Warning: this method does not work on levels (dB) + """ + # Settings + tr = 2 # truncate window after 2x std; shorter is faster and less accurate + Noct = sw.value[0] + assert Noct > 0, "'Noct' must be absolute positive" + if Noct < 1: + raise Warning('Check if \'Noct\' is entered correctly') + + # Initialize + L = len(freq) + Q = np.zeros(shape=(L, L), dtype=np.float16) # float16: keep size small + + x0 = 1 if freq[0] == 0 else 0 # Skip first data point if zero frequency + # Loop over indices of raw frequency vector + for x in range(x0, L): + # Find indices of data points to calculate current (smoothed) magnitude + # + # Indices beyond [0, L] point to non-existing data. Beyond 0 does not + # occur in this implementation. Beyond L occurs when the smoothing + # window nears the end of the series. + # If one end of the window is truncated, the other end + # could be truncated as well, to prevent an error on magnitude data + # with a slope. It however results in unsmoothed looking data at the + # end. + fc = freq[x] # center freq. of smoothing window + fl = fc / np.sqrt(2**(tr/Noct)) # lower cutoff + fu = fc * np.sqrt(2**(tr/Noct)) # upper cutoff + + # If the upper (frequency) side of the window is truncated because + # there is no data beyond the Nyquist frequency, also truncate the + # other side to keep it symmetric in a log(frequency) scale. + # So: fu / fc = fc / fl + fNq = freq[-1] + if fu > fNq: + fu = fNq # no data beyond fNq + fl = fc**2 / fu # keep window symmetric + + # Find indices corresponding to frequencies + xl = bisect.bisect_left(freq, fl) # index corresponding to fl + xu = bisect.bisect_left(freq, fu) + + xl = xu-1 if xu-xl <= 0 else xl # Guarantee window length of at least 1 + + # Calculate window + xg = np.arange(xl, xu) # indices + fg = freq[xg] # [Hz] corresponding freq + gs = np.sqrt( 1/ ((1+((fg/fc - fc/fg)*(1.507*Noct))**6)) ) # raw windw + gs /= np.sum(gs) # normalize: integral=1 + Q[x, xl:xu] = gs # add to matrix + + return Q + + def smoothSpectralData(freq, M, sw: SmoothingWidth, st: SmoothingType = SmoothingType.levels): """ @@ -92,6 +156,10 @@ def smoothSpectralData(freq, M, sw: SmoothingWidth, side. The deviation is largest when Noct is small (e.g. coarse smoothing). Casper Jansen, 07-05-2021 + Update 16-01-2023: speed up algorithm + - Smoothing is performed using matrix multiplication + - The smoothing matrix is not calculated if it already exists + Args: freq: array of frequencies of data points [Hz] - equally spaced M: array of either power, transfer functin or dB points. Depending on @@ -104,9 +172,6 @@ def smoothSpectralData(freq, M, sw: SmoothingWidth, """ # TODO: Make this function multi-dimensional array aware. - # Settings - tr = 2 # truncate window after 2x std; shorter is faster and less accurate - # Safety MM = copy.deepcopy(M) Noct = sw.value[0] @@ -135,169 +200,33 @@ def smoothSpectralData(freq, M, sw: SmoothingWidth, # P is power while smoothing. x are indices of P. Psm = np.zeros_like(P) # Smoothed power - to be calculated - x0 = 1 if freq[0] == 0 else 0 # Skip first data point if zero frequency - Psm[0] = P[0] # Reuse old value in case first data.. + if freq[0] == 0: + Psm[0] = P[0] # Reuse old value in case first data.. # ..point is skipped. Not plotted any way. - # Loop through data points - for x in range(x0, L): - # Find indices of data points to calculate current (smoothed) magnitude - # - # Indices beyond [0, L] point to non-existing data. Beyond 0 does not - # occur in this implementation. Beyond L occurs when the smoothing - # window nears the end of the series. - # If one end of the window is truncated, the other end - # could be truncated as well, to prevent an error on magnitude data - # with a slope. It however results in unsmoothed looking data at the - # end. - fc = freq[x] # center freq. of smoothing window - fl = fc / np.sqrt(2**(tr/Noct)) # lower cutoff - fu = fc * np.sqrt(2**(tr/Noct)) # upper cutoff + # Re-use smoothing matrix Q if available. Otherwise, calculate. + # Store in dict 'Qdict' + nfft = int(2*(len(freq)-1)) + key = f"nfft{nfft}_Noct{Noct}" # matrix name - # If the upper (frequency) side of the window is truncated because - # there is no data beyond the Nyquist frequency, also truncate the - # other side to keep it symmetric in a log(frequency) scale. - # So: fu / fc = fc / fl - fNq = freq[-1] - if fu > fNq: - fu = fNq # no data beyond fNq - fl = fc**2 / fu # keep window symmetric + if 'Qdict' not in globals(): # Guarantee Qdict exists + global Qdict + Qdict = {} - # Find indices corresponding to frequencies - xl = bisect.bisect_left(freq, fl) # index corresponding to fl - xu = bisect.bisect_left(freq, fu) + if key in Qdict: + Q = Qdict[key] + else: + Q = smoothCalcMatrix(freq, sw) + Qdict[key] = Q - # Guarantee window length of at least 1 - if xu - xl <= 0: - xl = xu - 1 - - # Calculate window - g = np.zeros(xu-xl) - for n, xi in enumerate(range(xl, xu)): - fi = freq[xi] # current frequency - g[n] = np.sqrt( 1/ ((1+((fi/fc - fc/fi)*(1.507*Noct))**6)) ) - g /= np.sum(g) # normalize: integral=1 - - # Apply smoothing - Psm[x] = np.dot(g, P[xl:xu]) + # Apply smoothing + Psm = np.matmul(Q, P) if st == SmoothingType.levels: Psm = 10*np.log10(Psm) return Psm -## OLD VERSION -#from scipy.signal.windows import gaussian -#def smoothSpectralData(freq, M, sw: SmoothingWidth, -# st: SmoothingType = SmoothingType.levels): -# """ -# Apply fractional octave smoothing to magnitude data in frequency domain. -# Smoothing is performed to power, using a sliding Gaussian window with -# variable length. The window is truncated after 2x std at either side. -# -# The implementation is not exact, because f is linearly spaced and -# fractional octave smoothing is related to log spaced data. In this -# implementation, the window extends with a fixed frequency step to either -# side. The deviation is largest when Noct is small (e.g. coarse smoothing). -# Casper Jansen, 07-05-2021 -# -# Args: -# freq: array of frequencies of data points [Hz] - equally spaced -# M: array of either power, transfer functin or dB points. Depending on -# the smoothing type `st`, the smoothing is applied. -# -# Returns: -# freq : array frequencies of data points [Hz] -# Msm : float smoothed magnitude of data points -# -# """ -# # TODO: Make this function multi-dimensional array aware. -# -# # Settings -# tr = 2 # truncate window after 2x std; shorter is faster and less accurate -# -# # Safety -# Noct = sw.value[0] -# assert Noct > 0, "'Noct' must be absolute positive" -# if Noct < 1: raise Warning('Check if \'Noct\' is entered correctly') -# assert len(freq)==len(M), 'f and M should have equal length' -# -# if st == SmoothingType.ps: -# assert np.min(M) >= 0, 'absolute magnitude M cannot be negative' -# if st == SmoothingType.levels and isinstance(M.dtype, complex): -# raise RuntimeError('Decibel input should be real-valued') -# -# # Initialize -# L = M.shape[0] # number of data points -# -# P = M -# if st == SmoothingType.levels: -# P = 10**(P/10) -# # TODO: This does not work due to complex numbers. Should be split up in -# # magnitude and phase. -# # elif st == SmoothingType.tf: -# # P = P**2 -# -# # P is power while smoothing. x are indices of P. -# Psm = np.zeros_like(P) # Smoothed power - to be calculated -# x0 = 1 if freq[0]==0 else 0 # Skip first data point if zero frequency -# Psm[0] = P[0] # Reuse old value in case first data.. -# # ..point is skipped. Not plotted any way. -# df = freq[1] - freq[0] # Frequency step -# -# # Loop through data points -# for x in range(x0, L): -# # Find indices of data points to calculate current (smoothed) magnitude -# fc = freq[x] # center freq. of smoothing window -# Df = tr * fc / Noct # freq. range of smoothing window -# -# # desired lower index of frequency array to be used during smoothing -# xl = int(np.ceil(x - 0.5*Df/df)) -# -# # upper index + 1 (because half-open interval) -# xu = int(np.floor(x + 0.5*Df/df)) + 1 -# -# # Create window, suitable for frequency lin-spaced data points -# Np = xu - xl # number of points -# std = Np / (2 * tr) -# wind = gaussian(Np, std) # Gaussian window -# -# # Clip indices to valid range -# # -# # Indices beyond [0, L] point to non-existing data. This occurs when -# # the smoothing windows nears the beginning or end of the series. -# # Optional: if one end of the window is clipped, the other end -# # could be clipped as well, to prevent an error on magnitude data with -# # a slope. It however results in unsmoothed looking data at the ends. -# if xl < 0: -# rl = 0 - xl # remove this number of points at the lower end -# xl = xl + rl # .. from f -# wind = wind[rl:] # .. and from window -# -# # rl = 0 - xl # remove this number of points at the lower end -# # xl = xl + rl # .. from f -# # xu = xu - rl -# # wind = wind[rl:-rl] # .. and from window -# -# if xu > L: -# ru = xu - L # remove this number of points at the upper end -# xu = xu - ru -# wind = wind[:-ru] -# -# # ru = xu - L # remove this number of points at the upper end -# # xl = xl + ru -# # xu = xu - ru -# # wind = wind[ru:-ru] -# -# # Apply smoothing -# wind_int = np.sum(wind) # integral -# Psm[x] = np.dot(wind, P[xl:xu]) / wind_int # apply window -# -# if st == SmoothingType.levels: -# Psm = 10*np.log10(Psm) -# -# return Psm - def smoothSpectralDataAlt(freq, Mdb, sw: SmoothingWidth): Noct = 1/sw.value[0] @@ -382,7 +311,7 @@ if __name__ == "__main__": Noct = 2 # Noct = 6 for 1/6 oct. smoothing # Create dummy data set 1: noise - fmin = 3e3 # [Hz] min freq + fmin = 1e3 # [Hz] min freq fmax = 24e3 # [Hz] max freq Ndata = 200 # number of data points freq = np.linspace(fmin, fmax, Ndata) # frequency points @@ -402,10 +331,13 @@ if __name__ == "__main__": class sw: value = [Noct] st = SmoothingType.levels # so data is given in dB + # st = SmoothingType.ps # so data is given in power + # Smooth Msm = smoothSpectralData(freq, M, sw, st) fsm = freq + # Plot - lin frequency plt.figure() plt.plot(freq, M, '.b') @@ -414,12 +346,14 @@ if __name__ == "__main__": plt.ylabel('magnitude') plt.xlim([100, fmax]) plt.title('lin frequency') + plt.legend(['Raw', 'Smooth']) # Plot - log frequency plt.figure() plt.semilogx(freq, M, '.b') - plt.semilogx(fsm, Msm, 'r') + plt.semilogx(fsm, Msm, 'r') plt.xlabel('f (Hz)') plt.ylabel('magnitude') plt.xlim([100, fmax]) plt.title('log frequency') + plt.legend(['Raw', 'Smooth 1']) diff --git a/test/scipy_sos_thirdoctavebank.py b/test/scipy_sos_thirdoctavebank.py new file mode 100644 index 0000000..70413eb --- /dev/null +++ b/test/scipy_sos_thirdoctavebank.py @@ -0,0 +1,77 @@ +""" +This script checks whether the 3rd octave equalizer is norm compliant. + +Created on Tue Dec 31 14:43:19 2019 +@author: anne +""" + +from lasp.filter import (ThirdOctaveBankDesigner) +import matplotlib.pyplot as plt +import numpy as np +from scipy.signal import freqz, iirdesign, sosfreqz, butter, cheby2 + + +plt.close('all') + + +# %% Run check +def check_one(fs): + """Check for fs = fs""" + # Create filter bank + nsections = 5 # filter property + designer = ThirdOctaveBankDesigner(fs) + + # Test + xmin = designer.nominal_txt_tox('25') # lowest highest band, nominal name + xmax = designer.nominal_txt_tox('20k') # highest band + + compliant = True # remains True if all filters comply + fnyq = fs/2 + fll = designer.fl(xmin) / 1.1 # plot limits + fuu = designer.fu(xmax) * 1.1 + fig = plt.figure() + plt.xlabel('f (Hz)') + plt.ylabel('Magnitude (dB)') + for x in range(xmin, xmax+1): # loop over bands + + fl = designer.fl(x) # lower cut off + fu = designer.fu(x) # upper + + fln = fl/fnyq # lower cutoff, normalized [0, 1] + fun = fu/fnyq # upper + + sos = designer.createSOSFilter(x) + fnorm, h = sosfreqz(sos, worN=2**18, whole=False) + freq = fnorm*fnyq/np.pi + + h_dB = 20*np.log10(np.abs(h) + np.finfo(float).eps) + if not designer.testStandardCompliance(x, freq, h_dB): + print('==============================') + print(designer.nominal_txt(x), ' does not comply\n') + compliant = False + + plt.semilogx(freq, h_dB) + + # Upper / lower limits + freqlim, ulim, llim = designer.band_limits(x) + plt.semilogx(freqlim, ulim) + plt.semilogx(freqlim, llim) + + plt.ylim(-5, 1) + plt.xlim(fll, fuu) + + print(f"Passed test : {compliant} for fs = {fs} Hz") + return compliant + + +if __name__ == "__main__": + # Make list of fs to be tested + # Taken from https://en.wikipedia.org/wiki/Sampling_(signal_processing) + fs_all = [8000, 11025, 16000, 22050, 32000, 37800, 44056, 44100, 47250, + 48000, 50000, 50400, 64000, 88200, 96000, 176400, 192000, 352800] + allcompliant = True + for fs in fs_all: + compliant = check_one(fs) + if not compliant: + allcompliant = False +print(f"\n\nAll passed test : {allcompliant}") diff --git a/test/test_aps.py b/test/test_aps.py index 4ec248f..1e78498 100644 --- a/test/test_aps.py +++ b/test/test_aps.py @@ -8,8 +8,6 @@ Created on Mon Jan 15 19:45:33 2018 import numpy as np from lasp import AvPowerSpectra, Window -import matplotlib.pyplot as plt -# plt.close('all') def test_aps1(): nfft = 16384 diff --git a/test/test_cppslm.py b/test/test_cppslm.py index 5890b8d..ab3dbb8 100644 --- a/test/test_cppslm.py +++ b/test/test_cppslm.py @@ -2,34 +2,35 @@ import numpy as np from lasp import cppSLM from lasp.filter import SPLFilterDesigner -import matplotlib.pyplot as plt + def test_cppslm1(): """ Generate a sine wave """ fs = 48000 - omg = 2*np.pi*1000 - - slm = cppSLM.fromBiquads(fs, 2e-5, 1, 0.125, [1.,0,0,1,0,0]) - - t = np.linspace(0, 10, 10*fs, endpoint=False) + omg = 2 * np.pi * 1000 + + slm = cppSLM.fromBiquads(fs, 2e-5, 1, 0.125, + np.array([[1., 0, 0, 1, 0, 0]]).T) + + t = np.linspace(0, 10, 10 * fs, endpoint=False) + + # Input signal with an rms of 1 Pa + in_ = np.sin(omg * t) * np.sqrt(2) - # Input signal with an rms of 1 Pa - in_ = np.sin(omg*t)*np.sqrt(2) - # Compute overall RMS - rms = np.sqrt(np.sum(in_**2)/in_.size) + rms = np.sqrt(np.sum(in_**2) / in_.size) # Compute overall level - level = 20*np.log10(rms/2e-5) - + level = 20 * np.log10(rms / 2e-5) + # Output of SLM out = slm.run(in_) - + # Output of SLM should be close to theoretical # level, at least for reasonable time constants # (Fast, Slow etc) - assert(np.isclose(out[-1,0], level)) + assert (np.isclose(out[-1, 0], level)) def test_cppslm2(): @@ -37,53 +38,57 @@ def test_cppslm2(): Generate a sine wave, now A-weighted """ fs = 48000 - omg = 2*np.pi*1000 - + omg = 2 * np.pi * 1000 + filt = SPLFilterDesigner(fs).A_Sos_design() - slm = cppSLM.fromBiquads(fs, 2e-5, 0, 0.125, filt.flatten(), [1.,0,0,1,0,0]) - - t = np.linspace(0, 10, 10*fs, endpoint=False) - - # Input signal with an rms of 1 Pa - in_ = np.sin(omg*t) *np.sqrt(2) - + slm = cppSLM.fromBiquads(fs, 2e-5, 0, 0.125, + filt.flatten(), # Pre-filter coefs + np.array([[1., 0, 0, 1, 0, 0]]).T # Bandpass coefs + ) + + t = np.linspace(0, 10, 10 * fs, endpoint=False) + + # Input signal with an rms of 1 Pa + in_ = np.sin(omg * t) * np.sqrt(2) + # Compute overall RMS - rms = np.sqrt(np.sum(in_**2)/in_.size) + rms = np.sqrt(np.sum(in_**2) / in_.size) # Compute overall level - level = 20*np.log10(rms/2e-5) - + level = 20 * np.log10(rms / 2e-5) + # Output of SLM out = slm.run(in_) # Output of SLM should be close to theoretical # level, at least for reasonable time constants # (Fast, Slow etc) - assert np.isclose(out[-1,0], level, atol=1e-2) + assert np.isclose(out[-1, 0], level, atol=1e-2) + def test_cppslm3(): fs = 48000 - omg = 2*np.pi*1000 - + omg = 2 * np.pi * 1000 + filt = SPLFilterDesigner(fs).A_Sos_design() - slm = cppSLM.fromBiquads(fs, 2e-5, 0, 0.125, filt.flatten(), [1.,0,0,1,0,0]) - t = np.linspace(0, 10, 10*fs, endpoint=False) - - in_ = 10*np.sin(omg*t) * np.sqrt(2)+np.random.randn() + slm = cppSLM.fromBiquads(fs, 2e-5, 0, 0.125, + filt.flatten(), + np.array([[1., 0, 0, 1, 0, 0]]).T) + t = np.linspace(0, 10, 10 * fs, endpoint=False) + + in_ = 10 * np.sin(omg * t) * np.sqrt(2) + np.random.randn() # Compute overall RMS - rms = np.sqrt(np.sum(in_**2)/in_.size) + rms = np.sqrt(np.sum(in_**2) / in_.size) # Compute overall level - level = 20*np.log10(rms/2e-5) - - + level = 20 * np.log10(rms / 2e-5) + # Output of SLM out = slm.run(in_) - - Lpeak = 20*np.log10(np.max(np.abs(in_)/2e-5)) + + Lpeak = 20 * np.log10(np.max(np.abs(in_) / 2e-5)) Lpeak slm.Lpeak() - assert np.isclose(out[-1,0], slm.Leq()[0][0], atol=1e-2) - assert np.isclose(Lpeak, slm.Lpeak()[0][0], atol=2e0) - + assert np.isclose(out[-1, 0], slm.Leq()[0], atol=1e-2) + assert np.isclose(Lpeak, slm.Lpeak()[0], atol=2e0) if __name__ == '__main__': diff --git a/test/test_ps.py b/test/test_ps.py index b320da9..cb238a5 100644 --- a/test/test_ps.py +++ b/test/test_ps.py @@ -5,8 +5,6 @@ Testing code for power spectra """ import numpy as np from lasp import PowerSpectra, Window -# import matplotlib.pyplot as plt -# plt.close('all') def test_ps(): """ diff --git a/test/test_slm.py b/test/test_slm.py deleted file mode 100644 index 1596262..0000000 --- a/test/test_slm.py +++ /dev/null @@ -1,37 +0,0 @@ -#!/usr/bin/python3 -import numpy as np -from lasp import SLM - - - -nframes = 0 -samplerate = 48000 -omg = 2*np.pi*1000 - - - -# def mycallback(input_, nframes, streamtime): -# t = np.linspace(streamtime, streamtime + nframes/samplerate, -# nframes)[np.newaxis,:] -# outp = 0.1*np.sin(omg*t) -# return outp, 0 - -# if __name__ == '__main__': -# pa = RtAudio() -# count = pa.getDeviceCount() -# # dev = pa.getDeviceInfo(0) -# for i in range(count): -# dev = pa.getDeviceInfo(i) -# print(dev) - -# outputparams = {'deviceid': 0, 'nchannels': 1, 'firstchannel': 0} -# pa.openStream(outputparams, None , Format_FLOAT64,samplerate, 512, mycallback) -# pa.startStream() - -# input() - -# pa.stopStream() -# pa.closeStream() - - -