From bd67c71c20dbd355d245b8a27df5607d009cb055 Mon Sep 17 00:00:00 2001 From: David Schote Date: Fri, 24 Jul 2026 08:44:08 +0200 Subject: [PATCH 1/3] Added .gitattributes --- .gitattributes | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) create mode 100644 .gitattributes diff --git a/.gitattributes b/.gitattributes new file mode 100644 index 00000000..ca876417 --- /dev/null +++ b/.gitattributes @@ -0,0 +1,18 @@ +# Store everything as LF in the repo, check out as LF everywhere +* text=auto eol=lf + +# Windows batch files genuinely need CRLF — cmd.exe mis-parses LF-only labels/goto +*.bat text eol=crlf +*.cmd text eol=crlf + +# Belt and braces for binaries (Git usually detects these, but be explicit) +*.png binary +*.jpg binary +*.gif binary +*.ico binary +*.pdf binary +*.zip binary +*.whl binary +*.pyd binary +*.so binary +*.dll binary From 874b4428dcd5dd0c5e3e5f5d11af4fd7c6770d5a Mon Sep 17 00:00:00 2001 From: David Schote Date: Fri, 24 Jul 2026 08:44:16 +0200 Subject: [PATCH 2/3] Added .editorconfig --- .editorconfig | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) create mode 100644 .editorconfig diff --git a/.editorconfig b/.editorconfig new file mode 100644 index 00000000..78b8832d --- /dev/null +++ b/.editorconfig @@ -0,0 +1,18 @@ +root = true + +[*] +end_of_line = lf +insert_final_newline = true +charset = utf-8 +trim_trailing_whitespace = true +indent_style = space +indent_size = 4 + +[*.{bat,cmd}] +end_of_line = crlf + +[*.{yml,yaml,json,toml}] +indent_size = 2 + +[*.md] +trim_trailing_whitespace = false From 604ed2b312ebd4d34cc4f6e242f688987672fc6e Mon Sep 17 00:00:00 2001 From: David Schote Date: Fri, 24 Jul 2026 08:47:34 +0200 Subject: [PATCH 3/3] Normalized line endings to LF --- .github/workflows/deployment.yml | 220 +-- .github/workflows/docs.yml | 98 +- .github/workflows/pytest.yml | 130 +- .github/workflows/static-tests.yml | 68 +- docs/make.bat | 70 +- examples/example_device_config.yaml | 134 +- src/console/interfaces/rx_data.py | 468 +++--- .../pulseq_interpreter/sequence_provider.py | 1324 ++++++++--------- src/console/spcm_control/rx_device.py | 1066 ++++++------- src/console/spcm_control/rx_processor.py | 224 +-- src/console/spcm_control/spcm/tools.py | 254 ++-- src/console/utilities/json_encoder.py | 52 +- tests/acquisition/test_ddc.py | 112 +- 13 files changed, 2110 insertions(+), 2110 deletions(-) diff --git a/.github/workflows/deployment.yml b/.github/workflows/deployment.yml index 09267991..59664c5e 100644 --- a/.github/workflows/deployment.yml +++ b/.github/workflows/deployment.yml @@ -1,111 +1,111 @@ -name: Build and Deploy - -on: - push: - branches: - - "**" # Match all branches, excludes tag pushes - -permissions: - contents: write - -jobs: - # Step 1: Build and upload package, runs on any git push - build-package: - name: Build package and upload artifact - runs-on: ubuntu-latest - outputs: - version: ${{ steps.get_version.outputs.version }} - - steps: - - name: Checkout code - uses: actions/checkout@v4 - with: - fetch-depth: 0 # Needed for full Git history and tags - - - name: Set up Python - uses: actions/setup-python@v5 - with: - python-version: '3.13' - cache: 'pip' - - - name: Install build tools - run: | - python -m pip install --upgrade build setuptools-git-versioning - - - name: Get version - id: get_version - run: | - # Assuming pyproject.toml is at the root of the repository. - VERSION=$(python -c "import setuptools_git_versioning; print(setuptools_git_versioning.get_version(root='.'))") - echo "Calculated Version: $VERSION" - echo "version=$VERSION" >> $GITHUB_OUTPUT - echo "VERSION=$VERSION" >> $GITHUB_ENV - - - name: Build package - run: python -m build - - - name: Upload distribution artifact - uses: actions/upload-artifact@v4 - with: - name: nexus-console-dist - path: dist/ - - # Step 2: Tag and Publish to PyPI, only runs on main - publish: - name: Tag and publish package to PyPI - runs-on: ubuntu-latest - needs: build-package - # Only run if on main branch and a version was determined for building - if: github.ref == 'refs/heads/main' && startsWith(needs.build-package.outputs.version, 'v') - - steps: - - name: Checkout code - uses: actions/checkout@v4 - with: - fetch-depth: 0 - - - name: Download package artifact - uses: actions/download-artifact@v4 - with: - name: nexus-console-dist - path: dist/ - - - name: Set Git user for tagging - run: | - git config user.name "github-actions" - git config user.email "github-actions@github.com" - - - name: Extract base version for tagging - id: tag_version_prep - run: | - RAW_VERSION="${{ needs.build-package.outputs.version }}" - TAG_VERSION="v$RAW_VERSION" - echo "Attempting to create tag: $TAG_VERSION" - echo "tag_version=$TAG_VERSION" >> $GITHUB_OUTPUT - - - name: Check if Git tag already exists - id: check_tag - run: | - TAG_VERSION="${{ steps.tag_version_prep.outputs.tag_version }}" - if git rev-parse "$TAG_VERSION" >/dev/null 2>&1; then - echo "Tag $TAG_VERSION already exists" - echo "tag_exists=true" >> $GITHUB_OUTPUT - else - echo "Tag $TAG_VERSION does not exist" - echo "tag_exists=false" >> $GITHUB_OUTPUT - fi - - - name: Create and push Git tag - if: steps.check_tag.outputs.tag_exists == 'false' - env: - GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} - run: | - TAG_VERSION="${{ steps.tag_version_prep.outputs.tag_version }}" - git tag $TAG_VERSION - git push origin $TAG_VERSION - - - name: Publish to PyPI - if: steps.check_tag.outputs.tag_exists == 'false' - uses: pypa/gh-action-pypi-publish@release/v1 - with: +name: Build and Deploy + +on: + push: + branches: + - "**" # Match all branches, excludes tag pushes + +permissions: + contents: write + +jobs: + # Step 1: Build and upload package, runs on any git push + build-package: + name: Build package and upload artifact + runs-on: ubuntu-latest + outputs: + version: ${{ steps.get_version.outputs.version }} + + steps: + - name: Checkout code + uses: actions/checkout@v4 + with: + fetch-depth: 0 # Needed for full Git history and tags + + - name: Set up Python + uses: actions/setup-python@v5 + with: + python-version: '3.13' + cache: 'pip' + + - name: Install build tools + run: | + python -m pip install --upgrade build setuptools-git-versioning + + - name: Get version + id: get_version + run: | + # Assuming pyproject.toml is at the root of the repository. + VERSION=$(python -c "import setuptools_git_versioning; print(setuptools_git_versioning.get_version(root='.'))") + echo "Calculated Version: $VERSION" + echo "version=$VERSION" >> $GITHUB_OUTPUT + echo "VERSION=$VERSION" >> $GITHUB_ENV + + - name: Build package + run: python -m build + + - name: Upload distribution artifact + uses: actions/upload-artifact@v4 + with: + name: nexus-console-dist + path: dist/ + + # Step 2: Tag and Publish to PyPI, only runs on main + publish: + name: Tag and publish package to PyPI + runs-on: ubuntu-latest + needs: build-package + # Only run if on main branch and a version was determined for building + if: github.ref == 'refs/heads/main' && startsWith(needs.build-package.outputs.version, 'v') + + steps: + - name: Checkout code + uses: actions/checkout@v4 + with: + fetch-depth: 0 + + - name: Download package artifact + uses: actions/download-artifact@v4 + with: + name: nexus-console-dist + path: dist/ + + - name: Set Git user for tagging + run: | + git config user.name "github-actions" + git config user.email "github-actions@github.com" + + - name: Extract base version for tagging + id: tag_version_prep + run: | + RAW_VERSION="${{ needs.build-package.outputs.version }}" + TAG_VERSION="v$RAW_VERSION" + echo "Attempting to create tag: $TAG_VERSION" + echo "tag_version=$TAG_VERSION" >> $GITHUB_OUTPUT + + - name: Check if Git tag already exists + id: check_tag + run: | + TAG_VERSION="${{ steps.tag_version_prep.outputs.tag_version }}" + if git rev-parse "$TAG_VERSION" >/dev/null 2>&1; then + echo "Tag $TAG_VERSION already exists" + echo "tag_exists=true" >> $GITHUB_OUTPUT + else + echo "Tag $TAG_VERSION does not exist" + echo "tag_exists=false" >> $GITHUB_OUTPUT + fi + + - name: Create and push Git tag + if: steps.check_tag.outputs.tag_exists == 'false' + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + run: | + TAG_VERSION="${{ steps.tag_version_prep.outputs.tag_version }}" + git tag $TAG_VERSION + git push origin $TAG_VERSION + + - name: Publish to PyPI + if: steps.check_tag.outputs.tag_exists == 'false' + uses: pypa/gh-action-pypi-publish@release/v1 + with: password: ${{ secrets.PYPI_API_TOKEN }} \ No newline at end of file diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml index d2fb1881..1be8ec10 100644 --- a/.github/workflows/docs.yml +++ b/.github/workflows/docs.yml @@ -1,49 +1,49 @@ -name: Sphinx Docs - -on: - push: - branches: - - main - - docu - pull_request: - -permissions: - contents: write - -jobs: - Documentation: - name: 'Build and deploy Documentation' - - runs-on: ubuntu-latest - - steps: - - name: Checkout - uses: actions/checkout@v4 - - - name: Set up Python - uses: actions/setup-python@v5 - with: - python-version: '3.13' - cache: 'pip' - - - name: Install dependencies - run: | - pip install --upgrade pip - pip install ".[docs]" - - - name: Build docs - run: sphinx-build -b html ./docs/source ./docs/build/html - - - name: Deploy to github pages - uses: peaceiris/actions-gh-pages@v3 - if: github.event_name != 'pull_request' - with: - publish_branch: github-pages - github_token: ${{ secrets.GITHUB_TOKEN }} - publish_dir: ./docs/build/html - - - name: Save documentation - uses: actions/upload-artifact@v4 - with: - name: Documentation - path: docs/build/html/ +name: Sphinx Docs + +on: + push: + branches: + - main + - docu + pull_request: + +permissions: + contents: write + +jobs: + Documentation: + name: 'Build and deploy Documentation' + + runs-on: ubuntu-latest + + steps: + - name: Checkout + uses: actions/checkout@v4 + + - name: Set up Python + uses: actions/setup-python@v5 + with: + python-version: '3.13' + cache: 'pip' + + - name: Install dependencies + run: | + pip install --upgrade pip + pip install ".[docs]" + + - name: Build docs + run: sphinx-build -b html ./docs/source ./docs/build/html + + - name: Deploy to github pages + uses: peaceiris/actions-gh-pages@v3 + if: github.event_name != 'pull_request' + with: + publish_branch: github-pages + github_token: ${{ secrets.GITHUB_TOKEN }} + publish_dir: ./docs/build/html + + - name: Save documentation + uses: actions/upload-artifact@v4 + with: + name: Documentation + path: docs/build/html/ diff --git a/.github/workflows/pytest.yml b/.github/workflows/pytest.yml index 3e1065bb..8a154a8d 100644 --- a/.github/workflows/pytest.yml +++ b/.github/workflows/pytest.yml @@ -1,65 +1,65 @@ -name: Pytest - -on: - pull_request: - push: - branches: - - main - -jobs: - test: - - runs-on: ubuntu-latest - - permissions: - pull-requests: write - contents: write - - steps: - - name: Checkout - uses: actions/checkout@v4 - - - name: Set up Python - uses: actions/setup-python@v5 - with: - python-version: '3.13' - cache: 'pip' - - - name: Install dependencies - run: | - pip install --upgrade pip - pip install .[test] - - - name: Install pytest-github-actions-annotate-failures plugin - run: pip install pytest-github-actions-annotate-failures - - - name: Run pytest - run: | - python -m pytest -n 4 --junitxml=pytest.xml --cov-report=term-missing:skip-covered --cov=console | tee pytest-coverage.txt - - - name: Pytest coverage comment - id: coverageComment - uses: MishaKav/pytest-coverage-comment@main - with: - pytest-coverage-path: ./pytest-coverage.txt - junitxml-path: ./pytest.xml - - - name: Create the Badge - uses: schneegans/dynamic-badges-action@v1.7.0 - with: - auth: ${{ secrets.GIST_SECRET }} - gistID: 4d47c22492a23337a79400f4859a4c25 - filename: coverage.json - label: Coverage Report - message: ${{ steps.coverageComment.outputs.coverage }} - color: ${{ steps.coverageComment.outputs.color }} - namedLogo: python - - - name: Set pipeline status - run: | - if [[ ${{ steps.coverageComment.outputs.errors }} -ne 0 || ${{ steps.coverageComment.outputs.failures }} -ne 0 ]]; then - echo "Errors or failures detected, marking pipeline as failure." - exit 1 - fi - - +name: Pytest + +on: + pull_request: + push: + branches: + - main + +jobs: + test: + + runs-on: ubuntu-latest + + permissions: + pull-requests: write + contents: write + + steps: + - name: Checkout + uses: actions/checkout@v4 + + - name: Set up Python + uses: actions/setup-python@v5 + with: + python-version: '3.13' + cache: 'pip' + + - name: Install dependencies + run: | + pip install --upgrade pip + pip install .[test] + + - name: Install pytest-github-actions-annotate-failures plugin + run: pip install pytest-github-actions-annotate-failures + + - name: Run pytest + run: | + python -m pytest -n 4 --junitxml=pytest.xml --cov-report=term-missing:skip-covered --cov=console | tee pytest-coverage.txt + + - name: Pytest coverage comment + id: coverageComment + uses: MishaKav/pytest-coverage-comment@main + with: + pytest-coverage-path: ./pytest-coverage.txt + junitxml-path: ./pytest.xml + + - name: Create the Badge + uses: schneegans/dynamic-badges-action@v1.7.0 + with: + auth: ${{ secrets.GIST_SECRET }} + gistID: 4d47c22492a23337a79400f4859a4c25 + filename: coverage.json + label: Coverage Report + message: ${{ steps.coverageComment.outputs.coverage }} + color: ${{ steps.coverageComment.outputs.color }} + namedLogo: python + + - name: Set pipeline status + run: | + if [[ ${{ steps.coverageComment.outputs.errors }} -ne 0 || ${{ steps.coverageComment.outputs.failures }} -ne 0 ]]; then + echo "Errors or failures detected, marking pipeline as failure." + exit 1 + fi + + diff --git a/.github/workflows/static-tests.yml b/.github/workflows/static-tests.yml index ccb737f8..91c8960c 100644 --- a/.github/workflows/static-tests.yml +++ b/.github/workflows/static-tests.yml @@ -1,34 +1,34 @@ -name: 'Static Tests' - -on: - pull_request: - -jobs: - linting: - - - runs-on: ubuntu-latest - - steps: - - name: Checkout - uses: actions/checkout@v4 - - - name: Set up Python - uses: actions/setup-python@v5 - with: - python-version: '3.13' - cache: 'pip' - - - name: Install dependencies - run: | - pip install --upgrade pip - pip install ".[lint]" - - - name: Run ruff check - run: ruff check - - - name: Run mypy - run: mypy src - - - +name: 'Static Tests' + +on: + pull_request: + +jobs: + linting: + + + runs-on: ubuntu-latest + + steps: + - name: Checkout + uses: actions/checkout@v4 + + - name: Set up Python + uses: actions/setup-python@v5 + with: + python-version: '3.13' + cache: 'pip' + + - name: Install dependencies + run: | + pip install --upgrade pip + pip install ".[lint]" + + - name: Run ruff check + run: ruff check + + - name: Run mypy + run: mypy src + + + diff --git a/docs/make.bat b/docs/make.bat index 747ffb7b..dc1312ab 100644 --- a/docs/make.bat +++ b/docs/make.bat @@ -1,35 +1,35 @@ -@ECHO OFF - -pushd %~dp0 - -REM Command file for Sphinx documentation - -if "%SPHINXBUILD%" == "" ( - set SPHINXBUILD=sphinx-build -) -set SOURCEDIR=source -set BUILDDIR=build - -%SPHINXBUILD% >NUL 2>NUL -if errorlevel 9009 ( - echo. - echo.The 'sphinx-build' command was not found. Make sure you have Sphinx - echo.installed, then set the SPHINXBUILD environment variable to point - echo.to the full path of the 'sphinx-build' executable. Alternatively you - echo.may add the Sphinx directory to PATH. - echo. - echo.If you don't have Sphinx installed, grab it from - echo.https://www.sphinx-doc.org/ - exit /b 1 -) - -if "%1" == "" goto help - -%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% -goto end - -:help -%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% - -:end -popd +@ECHO OFF + +pushd %~dp0 + +REM Command file for Sphinx documentation + +if "%SPHINXBUILD%" == "" ( + set SPHINXBUILD=sphinx-build +) +set SOURCEDIR=source +set BUILDDIR=build + +%SPHINXBUILD% >NUL 2>NUL +if errorlevel 9009 ( + echo. + echo.The 'sphinx-build' command was not found. Make sure you have Sphinx + echo.installed, then set the SPHINXBUILD environment variable to point + echo.to the full path of the 'sphinx-build' executable. Alternatively you + echo.may add the Sphinx directory to PATH. + echo. + echo.If you don't have Sphinx installed, grab it from + echo.https://www.sphinx-doc.org/ + exit /b 1 +) + +if "%1" == "" goto help + +%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% +goto end + +:help +%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% + +:end +popd diff --git a/examples/example_device_config.yaml b/examples/example_device_config.yaml index 00ff632a..c4fdd651 100644 --- a/examples/example_device_config.yaml +++ b/examples/example_device_config.yaml @@ -1,67 +1,67 @@ -# ======================================== -# Configuration file for the Nexus-console -# ======================================== - -# NOTE: This is just an example, parameters need to be adjusted for individual system - -# >> TX DEVICE: M2p.6546-x4 -> /dev/spcm1 (AWG) -# Total number of (analog) transmit channels: 4 -TxConfiguration: - # Could be extended to a list - device_path: "/dev/spcm1" - # Number of channels - max_available_channels: 4 - # Device sampling rate in MHz - sampling_rate: 20 - # Max. output amplitude per channel in mV - channel_max_amplitude: [200, 6000, 6000, 6000] - # Filter configuration of each transmit channel - channel_filter_type: [0, 2, 2, 2] - # Configure if RF/gradients are terminated into 50 ohm impedance (true == 50 ohm) - rf_terminated_50ohm: True - gradients_terminated_50ohm: False - # Calculate grad_to_volt per gradient channel: - # Gradient efficiency in T/m/A - gradient_efficiency: [0.37e-3, 0.451e-3, 0.4e-3] - # Gradient power amplifier in A/V - gpa_gain: [7.10, 7.10, 7.10] - # Scales RF in Hz (pulseq) to mV - rf_to_mvolt: 5.e-3 - # rf_gain_lut_path: "/path/to/lut.npy" - - -# >> RX DEVICE: M2p.5933-x4 -> /dev/spcm0 (Digitizer) -# Total number of receive channels: 8 -RxConfiguration: - # Could be extended to a list - device_path: "/dev/spcm0" - # Number of channels - max_available_channels: 8 - # Device sampling rate in MHz - sampling_rate: 20 - # Enable the receive channels - channel_enable: [1, 0, 0, 0, 0, 0, 0, 0] - # Set the max. amplitude per channel in mV, list length must at least match the number of enabled channels - channel_max_amplitude: [200, 200, 200, 200, 1000, 1000, 1000, 1000] - # Configure which channels are terminated into 50 ohm impedance (true == 50 ohm), list length must at least match the number of enabled channels - channel_terminated_50ohm: [True, True, True, True, False, False, False, False] - -SystemLimits: - # Maximum gradient amplitude in Hz/m (PyPulseq system stores this unit) - max_grad: 1703040.0 # equals 40 mT/m - # Maximum slew rate in Hz/m/s (PyPulseq system stores this unit) - max_slew: 2128800000.0 # equals 100 T/m/s - # Dead time at the beginning of RF event in s, covered by rf_delay - rf_dead_time: 20.e-6 - # Time delay at the end of an RF event in s - rf_ringdown_time: 30.e-6 - # Time delay at the beginning of an ADC event in s - adc_dead_time: 0. - # Raster time for block durations - block_duration_raster: 1.e-6 - # Raster time for RF pulses. - rf_raster_time: 1.e-6 - # Raster time for gradient waveforms. - grad_raster_time: 1.e-6 - # Raster time for ADC readout. - adc_raster_time: 1.e-6 +# ======================================== +# Configuration file for the Nexus-console +# ======================================== + +# NOTE: This is just an example, parameters need to be adjusted for individual system + +# >> TX DEVICE: M2p.6546-x4 -> /dev/spcm1 (AWG) +# Total number of (analog) transmit channels: 4 +TxConfiguration: + # Could be extended to a list + device_path: "/dev/spcm1" + # Number of channels + max_available_channels: 4 + # Device sampling rate in MHz + sampling_rate: 20 + # Max. output amplitude per channel in mV + channel_max_amplitude: [200, 6000, 6000, 6000] + # Filter configuration of each transmit channel + channel_filter_type: [0, 2, 2, 2] + # Configure if RF/gradients are terminated into 50 ohm impedance (true == 50 ohm) + rf_terminated_50ohm: True + gradients_terminated_50ohm: False + # Calculate grad_to_volt per gradient channel: + # Gradient efficiency in T/m/A + gradient_efficiency: [0.37e-3, 0.451e-3, 0.4e-3] + # Gradient power amplifier in A/V + gpa_gain: [7.10, 7.10, 7.10] + # Scales RF in Hz (pulseq) to mV + rf_to_mvolt: 5.e-3 + # rf_gain_lut_path: "/path/to/lut.npy" + + +# >> RX DEVICE: M2p.5933-x4 -> /dev/spcm0 (Digitizer) +# Total number of receive channels: 8 +RxConfiguration: + # Could be extended to a list + device_path: "/dev/spcm0" + # Number of channels + max_available_channels: 8 + # Device sampling rate in MHz + sampling_rate: 20 + # Enable the receive channels + channel_enable: [1, 0, 0, 0, 0, 0, 0, 0] + # Set the max. amplitude per channel in mV, list length must at least match the number of enabled channels + channel_max_amplitude: [200, 200, 200, 200, 1000, 1000, 1000, 1000] + # Configure which channels are terminated into 50 ohm impedance (true == 50 ohm), list length must at least match the number of enabled channels + channel_terminated_50ohm: [True, True, True, True, False, False, False, False] + +SystemLimits: + # Maximum gradient amplitude in Hz/m (PyPulseq system stores this unit) + max_grad: 1703040.0 # equals 40 mT/m + # Maximum slew rate in Hz/m/s (PyPulseq system stores this unit) + max_slew: 2128800000.0 # equals 100 T/m/s + # Dead time at the beginning of RF event in s, covered by rf_delay + rf_dead_time: 20.e-6 + # Time delay at the end of an RF event in s + rf_ringdown_time: 30.e-6 + # Time delay at the beginning of an ADC event in s + adc_dead_time: 0. + # Raster time for block durations + block_duration_raster: 1.e-6 + # Raster time for RF pulses. + rf_raster_time: 1.e-6 + # Raster time for gradient waveforms. + grad_raster_time: 1.e-6 + # Raster time for ADC readout. + adc_raster_time: 1.e-6 diff --git a/src/console/interfaces/rx_data.py b/src/console/interfaces/rx_data.py index 37121f0a..67b279f3 100644 --- a/src/console/interfaces/rx_data.py +++ b/src/console/interfaces/rx_data.py @@ -1,234 +1,234 @@ -""""Define the dataclass and processing of receiver data.""" -from dataclasses import asdict, dataclass, field -from multiprocessing.shared_memory import SharedMemory - -import numpy as np -from scipy import signal - -from console.interfaces.acquisition_parameter import DDCMethod -from console.utilities import ddc - - -@dataclass -class RxData: - """Receive data object containing both the data and metadata of each receive event.""" - - # Rx data event number - index: int - - # Number of samples after decimation - num_samples: int - # Number of raw samples before decimation - num_samples_raw: int - # Number of samples to be discarded before and after ADC, defined by dead time - num_samples_discard: int - # Dwell time of decimated data in s - dwell_time: float - # Dwell time of undecimated data in s - dwell_time_raw: float - # Phase offset from sequence definition in rad - phase_offset: float - # Frequency offset from sequence definition in Hz - freq_offset: float - - # Averages tracking - total_averages: int - average_index: int = 0 - - # ADC labels - labels: dict[str, int | None] | None = None - - # Used for demodulation, value set in post init - decimation_factor: int = field(init=False) - - # Set the larmor frequency in Hz for each object - larmor_frequency: None | float = None - - # Frequency in Hz with which the data are demodulated - demod_frequency: None | float = None - - # Frequency in Hz with which the reference signal is demodulated - phase_ref_frequency: float | None = None - - # Set the default demodulation method to FIR - ddc_method: DDCMethod = DDCMethod.FIR - - # Scaling factor for each receive channel - scaling_factor: None | np.ndarray | list[float] = None - - # Raw data is the raw data coming from the Rx cards, prior to demodulation and decimation - # Shape of Raw data is (num_channels_enabled, raw number of samples) - raw_data: None | np.ndarray = None - - # Phase reference signal - phase_reference: None | np.ndarray = None - - # Timestamp in s of start of data acquisition relative to sequence execution start - time_stamp: None | float = None - - # Processed data is the demodulated, phased and decimated data - processed_data: None | np.ndarray = None - - # Shared memory handle and (name, shape) needed to reattach after pickling. - # track=False on both sides: the worker unlinks explicitly, avoiding a - # spurious 'leaked shared_memory' warning from the resource tracker. - _shm: SharedMemory | None = field(default=None, init=False, repr=False, compare=False) - _shm_meta: tuple[str, tuple[int, ...]] | None = field(default=None, init=False, repr=False, compare=False) - - def __post_init__(self) -> None: - """Post init method to calculate the decimation factor.""" - self.decimation_factor = round(self.dwell_time / self.dwell_time_raw) - - def __getstate__(self) -> dict: - """Strip the live handle and raw_data view; keep the shm name/shape (custom pickle).""" - state = self.__dict__.copy() - if self._shm is not None: - state["raw_data"] = None - state["_shm"] = None - return state - - def __setstate__(self, state: dict) -> None: - """Reattach to shared memory if a name/shape pair is present (custom pickle).""" - self.__dict__.update(state) - if self._shm_meta is not None: - name, shape = self._shm_meta - try: - self._shm = SharedMemory(name=name, create=False, track=False) - self.raw_data = np.ndarray(shape, dtype=np.int16, buffer=self._shm.buf) - except Exception as e: - raise RuntimeError(f"Failed to reattach to shared memory '{name}': {e}") from e - - def write_raw_data(self, data: np.ndarray) -> None: - """Copy ADC data into shared memory, allocating it on the first call. - - Parameters - ---------- - data - Raw int16 ADC data with shape (num_channels, num_samples_raw). - """ - if self._shm is None: - self._shm = SharedMemory(create=True, size=data.nbytes, track=False) - self._shm_meta = (self._shm.name, data.shape) - self.raw_data = np.ndarray(data.shape, dtype=np.int16, buffer=self._shm.buf) - if self.raw_data is None: - raise RuntimeError("Shared memory buffer not initialized.") - self.raw_data[:] = data - - def materialize(self, keep: bool = False) -> None: - """Finalize raw_data after processing, releasing shared memory. - - Parameters - ---------- - keep - If True, copy raw_data to a regular numpy array before releasing - shared memory so the data survives. If False, raw_data is set to None. - """ - if self._shm is None: - return - regular_array = self.raw_data.copy() if keep and self.raw_data is not None else None - self._shm.close() - self._shm.unlink() - self._shm = None - self._shm_meta = None - self.raw_data = regular_array - - def __str__(self) -> str: - """Return string representation of information contained within RxData class.""" - lines = ["RxData:"] - lines.append("-" * 7) - for key, value in self.dict().items(): - lines.append(f"{key:<20}: {value}") - return "\n".join(lines) - - def dict(self) -> dict: - """Return RxData meta information as string.""" - return { - key: ( - value.shape if key in ("processed_data", "raw_data") and value is not None - else "None" if value is None - else value - ) - for key, value in asdict(self).items() - if not key.startswith("_") - } - - def decimate_data(self, data) -> np.ndarray: - """Decimate the data using the defined `DDCMethod` method.""" - if self.decimation_factor <= 1 or not isinstance(self.decimation_factor, int): - raise ValueError(f"Invalid decimation factor {self.decimation_factor}") - - # Recover 50% amplitude loss from filtering the 2*f_Larmor mixing term; at f=0, LO is 1 (identity) - scaling = 1. - if self.larmor_frequency is not None and self.larmor_frequency > 0.: - scaling = 2. - - match self.ddc_method: - case DDCMethod.CIC: - return scaling * ddc.filter_cic_fir_comp(data, decimation=self.decimation_factor, number_of_stages=5) - case DDCMethod.AVG: - return scaling * ddc.filter_moving_average(data, decimation=self.decimation_factor, overlap=8) - case _: - # Default case is FIR decimation - return scaling * signal.decimate(data, q=self.decimation_factor, ftype="fir", axis=-1) - - def demod_and_phase_data(self, data) -> np.ndarray: - """Demodulate and phase the data contained in raw_data. - - This step first demodulates the acquired data using the demodulation frequency, - which is usually the Larmor frequency. If a phase reference has been acquired, - the phase reference signal is demodulated at the phase reference frequency. - The phase correction term calculated from the phase reference is used to correct the - acquired MR data. In a last step the phase offset defined by the sequence is applied. - """ - if self.demod_frequency is None: - raise RuntimeError("Demodulation frequency not set") - # Demodulate the data - time = np.arange(np.size(data, -1)) * self.dwell_time_raw - data_demod = data * np.exp(-2j * np.pi * time * self.demod_frequency) - - # Demodulate the reference signal if available and correct acquired data - if self.phase_reference is not None and self.phase_ref_frequency is not None: - # Demodulation of the phase reference signal - time_reference = np.arange(self.phase_reference.size) * self.dwell_time_raw - ref_demod = self.phase_reference * np.exp(-2j * np.pi * self.phase_ref_frequency * time_reference) - # Calculation of the phase correction term for the acquired MR data - phase_correction = np.angle(np.sum(ref_demod)) * (self.demod_frequency / self.phase_ref_frequency) - # Apply phase correction in place - data_demod *= np.exp(-1j * phase_correction) - - # Apply receive phase offset to data and return data - return data_demod * np.exp(1j * self.phase_offset) - - def scale_data(self, data) -> np.ndarray: - """Scale the receive data to go from ADC units to mV.""" - if self.scaling_factor is not None: - return data * np.expand_dims(self.scaling_factor, axis=-1) - else: - # If no scaling data is provided then just return the array as an array of floats for consistency - return data.astype(float) - - def process_data(self, store_unprocessed: bool = True) -> None: - """Process (demodulate, phase and downsample) the raw data contained in the rx object.""" - if self.larmor_frequency is None: - raise RuntimeError("Larmor frequency not set, please set prior to processing data") - - if self.raw_data is None: - raise RuntimeError("Can't process data; No raw data present in RxData object") - - if np.size(self.raw_data, axis=-1) != self.num_samples_raw: - raise ValueError(f"Number of collected samples is different from expected: " - f"{np.size(self.raw_data, axis = -1)} collected vs {self.num_samples_raw} expected") - - self.demod_frequency = self.larmor_frequency + self.freq_offset - scaled_data = self.scale_data(self.raw_data) - demod_data = self.demod_and_phase_data(scaled_data) - - # Creating the processed data output array first and copying the values of the output of the decimation - # avoids an apparent memory leak when using the scipy.decimate with the 'iir' ftype - # Note that the processed data may contain samples from pre and post sampling - output_shape = (*np.shape(demod_data)[:-1], self.num_samples + int(2*self.num_samples_discard)) - self.processed_data = np.zeros(output_shape, dtype=complex) - self.processed_data[:] = self.decimate_data(demod_data) - - if not store_unprocessed: - self.raw_data = None +""""Define the dataclass and processing of receiver data.""" +from dataclasses import asdict, dataclass, field +from multiprocessing.shared_memory import SharedMemory + +import numpy as np +from scipy import signal + +from console.interfaces.acquisition_parameter import DDCMethod +from console.utilities import ddc + + +@dataclass +class RxData: + """Receive data object containing both the data and metadata of each receive event.""" + + # Rx data event number + index: int + + # Number of samples after decimation + num_samples: int + # Number of raw samples before decimation + num_samples_raw: int + # Number of samples to be discarded before and after ADC, defined by dead time + num_samples_discard: int + # Dwell time of decimated data in s + dwell_time: float + # Dwell time of undecimated data in s + dwell_time_raw: float + # Phase offset from sequence definition in rad + phase_offset: float + # Frequency offset from sequence definition in Hz + freq_offset: float + + # Averages tracking + total_averages: int + average_index: int = 0 + + # ADC labels + labels: dict[str, int | None] | None = None + + # Used for demodulation, value set in post init + decimation_factor: int = field(init=False) + + # Set the larmor frequency in Hz for each object + larmor_frequency: None | float = None + + # Frequency in Hz with which the data are demodulated + demod_frequency: None | float = None + + # Frequency in Hz with which the reference signal is demodulated + phase_ref_frequency: float | None = None + + # Set the default demodulation method to FIR + ddc_method: DDCMethod = DDCMethod.FIR + + # Scaling factor for each receive channel + scaling_factor: None | np.ndarray | list[float] = None + + # Raw data is the raw data coming from the Rx cards, prior to demodulation and decimation + # Shape of Raw data is (num_channels_enabled, raw number of samples) + raw_data: None | np.ndarray = None + + # Phase reference signal + phase_reference: None | np.ndarray = None + + # Timestamp in s of start of data acquisition relative to sequence execution start + time_stamp: None | float = None + + # Processed data is the demodulated, phased and decimated data + processed_data: None | np.ndarray = None + + # Shared memory handle and (name, shape) needed to reattach after pickling. + # track=False on both sides: the worker unlinks explicitly, avoiding a + # spurious 'leaked shared_memory' warning from the resource tracker. + _shm: SharedMemory | None = field(default=None, init=False, repr=False, compare=False) + _shm_meta: tuple[str, tuple[int, ...]] | None = field(default=None, init=False, repr=False, compare=False) + + def __post_init__(self) -> None: + """Post init method to calculate the decimation factor.""" + self.decimation_factor = round(self.dwell_time / self.dwell_time_raw) + + def __getstate__(self) -> dict: + """Strip the live handle and raw_data view; keep the shm name/shape (custom pickle).""" + state = self.__dict__.copy() + if self._shm is not None: + state["raw_data"] = None + state["_shm"] = None + return state + + def __setstate__(self, state: dict) -> None: + """Reattach to shared memory if a name/shape pair is present (custom pickle).""" + self.__dict__.update(state) + if self._shm_meta is not None: + name, shape = self._shm_meta + try: + self._shm = SharedMemory(name=name, create=False, track=False) + self.raw_data = np.ndarray(shape, dtype=np.int16, buffer=self._shm.buf) + except Exception as e: + raise RuntimeError(f"Failed to reattach to shared memory '{name}': {e}") from e + + def write_raw_data(self, data: np.ndarray) -> None: + """Copy ADC data into shared memory, allocating it on the first call. + + Parameters + ---------- + data + Raw int16 ADC data with shape (num_channels, num_samples_raw). + """ + if self._shm is None: + self._shm = SharedMemory(create=True, size=data.nbytes, track=False) + self._shm_meta = (self._shm.name, data.shape) + self.raw_data = np.ndarray(data.shape, dtype=np.int16, buffer=self._shm.buf) + if self.raw_data is None: + raise RuntimeError("Shared memory buffer not initialized.") + self.raw_data[:] = data + + def materialize(self, keep: bool = False) -> None: + """Finalize raw_data after processing, releasing shared memory. + + Parameters + ---------- + keep + If True, copy raw_data to a regular numpy array before releasing + shared memory so the data survives. If False, raw_data is set to None. + """ + if self._shm is None: + return + regular_array = self.raw_data.copy() if keep and self.raw_data is not None else None + self._shm.close() + self._shm.unlink() + self._shm = None + self._shm_meta = None + self.raw_data = regular_array + + def __str__(self) -> str: + """Return string representation of information contained within RxData class.""" + lines = ["RxData:"] + lines.append("-" * 7) + for key, value in self.dict().items(): + lines.append(f"{key:<20}: {value}") + return "\n".join(lines) + + def dict(self) -> dict: + """Return RxData meta information as string.""" + return { + key: ( + value.shape if key in ("processed_data", "raw_data") and value is not None + else "None" if value is None + else value + ) + for key, value in asdict(self).items() + if not key.startswith("_") + } + + def decimate_data(self, data) -> np.ndarray: + """Decimate the data using the defined `DDCMethod` method.""" + if self.decimation_factor <= 1 or not isinstance(self.decimation_factor, int): + raise ValueError(f"Invalid decimation factor {self.decimation_factor}") + + # Recover 50% amplitude loss from filtering the 2*f_Larmor mixing term; at f=0, LO is 1 (identity) + scaling = 1. + if self.larmor_frequency is not None and self.larmor_frequency > 0.: + scaling = 2. + + match self.ddc_method: + case DDCMethod.CIC: + return scaling * ddc.filter_cic_fir_comp(data, decimation=self.decimation_factor, number_of_stages=5) + case DDCMethod.AVG: + return scaling * ddc.filter_moving_average(data, decimation=self.decimation_factor, overlap=8) + case _: + # Default case is FIR decimation + return scaling * signal.decimate(data, q=self.decimation_factor, ftype="fir", axis=-1) + + def demod_and_phase_data(self, data) -> np.ndarray: + """Demodulate and phase the data contained in raw_data. + + This step first demodulates the acquired data using the demodulation frequency, + which is usually the Larmor frequency. If a phase reference has been acquired, + the phase reference signal is demodulated at the phase reference frequency. + The phase correction term calculated from the phase reference is used to correct the + acquired MR data. In a last step the phase offset defined by the sequence is applied. + """ + if self.demod_frequency is None: + raise RuntimeError("Demodulation frequency not set") + # Demodulate the data + time = np.arange(np.size(data, -1)) * self.dwell_time_raw + data_demod = data * np.exp(-2j * np.pi * time * self.demod_frequency) + + # Demodulate the reference signal if available and correct acquired data + if self.phase_reference is not None and self.phase_ref_frequency is not None: + # Demodulation of the phase reference signal + time_reference = np.arange(self.phase_reference.size) * self.dwell_time_raw + ref_demod = self.phase_reference * np.exp(-2j * np.pi * self.phase_ref_frequency * time_reference) + # Calculation of the phase correction term for the acquired MR data + phase_correction = np.angle(np.sum(ref_demod)) * (self.demod_frequency / self.phase_ref_frequency) + # Apply phase correction in place + data_demod *= np.exp(-1j * phase_correction) + + # Apply receive phase offset to data and return data + return data_demod * np.exp(1j * self.phase_offset) + + def scale_data(self, data) -> np.ndarray: + """Scale the receive data to go from ADC units to mV.""" + if self.scaling_factor is not None: + return data * np.expand_dims(self.scaling_factor, axis=-1) + else: + # If no scaling data is provided then just return the array as an array of floats for consistency + return data.astype(float) + + def process_data(self, store_unprocessed: bool = True) -> None: + """Process (demodulate, phase and downsample) the raw data contained in the rx object.""" + if self.larmor_frequency is None: + raise RuntimeError("Larmor frequency not set, please set prior to processing data") + + if self.raw_data is None: + raise RuntimeError("Can't process data; No raw data present in RxData object") + + if np.size(self.raw_data, axis=-1) != self.num_samples_raw: + raise ValueError(f"Number of collected samples is different from expected: " + f"{np.size(self.raw_data, axis = -1)} collected vs {self.num_samples_raw} expected") + + self.demod_frequency = self.larmor_frequency + self.freq_offset + scaled_data = self.scale_data(self.raw_data) + demod_data = self.demod_and_phase_data(scaled_data) + + # Creating the processed data output array first and copying the values of the output of the decimation + # avoids an apparent memory leak when using the scipy.decimate with the 'iir' ftype + # Note that the processed data may contain samples from pre and post sampling + output_shape = (*np.shape(demod_data)[:-1], self.num_samples + int(2*self.num_samples_discard)) + self.processed_data = np.zeros(output_shape, dtype=complex) + self.processed_data[:] = self.decimate_data(demod_data) + + if not store_unprocessed: + self.raw_data = None diff --git a/src/console/pulseq_interpreter/sequence_provider.py b/src/console/pulseq_interpreter/sequence_provider.py index e9c33968..f611c7b5 100644 --- a/src/console/pulseq_interpreter/sequence_provider.py +++ b/src/console/pulseq_interpreter/sequence_provider.py @@ -1,662 +1,662 @@ -"""Sequence provider class.""" -import logging -from collections.abc import Callable -from dataclasses import dataclass -from math import floor -from pathlib import Path -from types import SimpleNamespace -from typing import Any - -import numpy as np -from pypulseq.opts import Opts -from pypulseq.Sequence.sequence import Sequence -from scipy.signal import resample - -from console.interfaces.acquisition_parameter import AcquisitionParameter -from console.interfaces.dimensions import Dimensions -from console.interfaces.rx_data import RxData -from console.interfaces.unrolled_sequence import UnrolledSequence - -try: - from line_profiler import profile -except ImportError: - def profile(func: Callable[..., Any]) -> Callable[..., Any]: - """Define placeholder for profile decorator.""" - return func - - -INT16_MAX = np.iinfo(np.int16).max -INT16_MIN = np.iinfo(np.int16).min - -NUM_REFERENCE_SAMPLES = 1000 -REFERENCE_FREQUENCY = 1.095e6 - - -@dataclass -class ADCGate: - """Define precalculated attributes of an ADC gate.""" - - start: int - num_samples_discard: int - num_samples_raw: int - -class SequenceProvider(Sequence): - """Sequence provider class. - - This object is inherited from pulseq sequence object, so that all methods of the - pypulseq ``Sequence`` object can be accessed. - - The main functionality of the ``SequenceProvider`` is to unroll a given pulseq sequence. - Usually the first step is to read a sequence file. The unrolling step can be achieved using - the ``unroll_sequence()`` function. - - Example - ------- - >>> seq_provider = SequenceProvider() - >>> seq_provider.read("./seq_file.seq") - >>> unrolled = seq_provider.unroll_sequence(acquisition_parameter) - """ - - __name__: str = "SequenceProvider" - - def __init__( - self, - gradient_efficiency: tuple[float, float, float], - gpa_gain: tuple[float, float, float], - gradient_output_limits: tuple[int, int, int], - gradients_50ohms: bool, - rf_output_limit: int, - rf_50ohms: bool, - rf_to_mvolt: float, - spcm_dwell_time: float, - system: Opts, - rf_gain_lut_path: Path | None = None, - ): - """Initialize sequence provider class which is used to unroll a pulseq sequence. - - Parameters - ---------- - gradient_efficiency - Efficiency of the gradient coils in mT/m/A, e.g. [0.4e-3, 0.4e-3, 0.4e-3]. - gpa_gain - Gain factor of the GPA per gradient channel, e.g. [4.7, 4.7, 4.7]. - gradient_output_limits - Integer output limit per gradient channel in mV, e.g. [6000, 6000, 6000]. - gradients_50ohms - Boolean flag which indicates if the gradient output is terminated into 50 ohms or high impedance. - If terminated into high impedance, the card output doubles, - what needs to be considered when calculating the sequence. - rf_output_limit - Integer output limit of the RF channel in mV. - rf_50ohms - Boolean flag which indicates if the rf output is terminated into 50 ohms (see gradients_50ohms). - rf_to_mvolt - Translation of RF waveform from pulseq (Hz) to mV. - spcm_dwell_time - Sampling time raster of the output waveform (depends on spectrum card). - system_limits - Absolute maximum system limits defined in the device configuration. - Used to instantiate the pypulseq `Opts()` class. - """ - if not isinstance(system, Opts): - raise AttributeError("Invalid system: Pypulseq `Opts` definition required.") - super().__init__(system=system) - self.log = logging.getLogger("SeqProv") - - # Set class instance attributes - self.rf_to_mvolt = rf_to_mvolt - self.spcm_dwell_time = spcm_dwell_time - self.spcm_freq = 1 / spcm_dwell_time - self.gpa_gain = gpa_gain - self.grad_eff = gradient_efficiency - - # Scale output limit dependent on high impedance flags: - # If output is terminated into high impedance (flag is true), the channel output is doubled. - # Otherwise, if output is terminated into 50 ohms impedance, output limit remains unchanged. - self.rf_out_limit = rf_output_limit if rf_50ohms else int(2 * rf_output_limit) - # Ensure tuple[int, int, int] for mypy typing - self.gradient_out_limits = ( - gradient_output_limits[0] if gradients_50ohms else int(2 * gradient_output_limits[0]), - gradient_output_limits[1] if gradients_50ohms else int(2 * gradient_output_limits[1]), - gradient_output_limits[2] if gradients_50ohms else int(2 * gradient_output_limits[2]), - ) - - # Setup phase reference signal - time = np.arange(NUM_REFERENCE_SAMPLES) * self.spcm_dwell_time - signal = np.exp(2j * np.pi * REFERENCE_FREQUENCY * time) - self.phase_reference = np.zeros(NUM_REFERENCE_SAMPLES, dtype=np.uint16) - self.phase_reference[signal > 0] = np.uint16(2**15) - - # Load LUT for RF gain correction if path to LUT is provided - self._rf_gain_lut: np.ndarray | None = None - self._load_rfpa_lut(rf_gain_lut_path) - - - # -------- PyPulseq interface -------- # - - def from_pypulseq(self, seq: Sequence) -> None: - """Read a pypulseq sequence to sequence provider. - - Parameters - ---------- - seq - Pypulseq ``Sequence`` instance - - Raises - ------ - AttributeError - seq is not a valid pypulseq ``Sequence`` instance - """ - if not isinstance(seq, Sequence): - raise AttributeError("Invalid sequence.") - # Re-initialize the parent to start from a clean pypulseq sequence - super().__init__(system=self.system) - for block_index, _ in seq.block_events.items(): - block = seq.get_block(block_index) - self.add_block(block) - # Set definitions without overwriting existing ones - for key, value in seq.definitions.items(): - self.definitions.setdefault(key, value) - - def to_pypulseq(self) -> Sequence | None: - """Create a pypulseq sequence from sequence provider.""" - seq = Sequence(system=self.system) - for block_index, _ in self.block_events.items(): - block = self.get_block(block_index) - seq.add_block(block) - seq.definitions = self.definitions - return seq - - # -------- Public interface -------- # - - def dict(self) -> dict: - """Abstract method which returns variables for logging in dictionary.""" - return { - "rf_to_mvolt": self.rf_to_mvolt, - "spcm_freq": self.spcm_freq, - "spcm_dwell_time": self.spcm_dwell_time, - "gpa_gain": self.gpa_gain, - "gradient_efficiency": self.grad_eff, - "output_limits": self.gradient_out_limits, - } - - @profile - def unroll_sequence(self, parameter: AcquisitionParameter) -> UnrolledSequence: - """Unroll the pypulseq sequence description. - - TODO: Reduce complexity. - - Parameters - ---------- - parameter - Instance of AcquisitionParameter containing all necessary parameters to - calculate the sequence waveforms, i.e. larmor frequency, gradient offsets, etc. - - Returns - ------- - UnrolledSequence - Instance of an unrolled sequence object which contains a list of numpy arrays with - the block-wise calculated sample points in correct spectrum card order (Fortran). - - The list of unrolled sequence arrays is returned as uint16 values which contain a digital - signal encoded by 15th bit. Only the RF channel does not contain a digital signal. - In addition, all receive events are described and returned in a list within the unrolled - sequence object. - - Examples - -------- - For channels ch0, ch1, ch2, ch3, data values n = 0, 1, ..., N are ordered the following way. - - >>> data = [ch0_0, ch1_0, ch2_0, ch3_0, ch0_1, ch1_1, ..., ch0_n, ..., ch3_N] - - Per channel data can be extracted by the following code. - - >>> rf = seq[0::4] - >>> gx = (seq[1::4] << 1).astype(np.int16) - >>> gy = (seq[2::4] << 1).astype(np.int16) - >>> gz = (seq[3::4] << 1).astype(np.int16) - - All the gradient channels contain a digital signal encoded by the 15th bit. - - `gx`: ADC gate signal - - `gy`: Reference signal for phase correction - - `gz`: RF unblanking signal - The following example shows, how to extract the digital signals - - >>> adc_gate = seq[1::4].astype(np.uint16) >> 15 - >>> reference = seq[2::4].astype(np.uint16) >> 15 - >>> unblanking = seq[3::4].astype(np.uint16) >> 15 - - As the 15th bit is not encoding the sign (as usual for int16), the values are casted to uint16 before shifting. - """ - try: - self._check_parameter(parameter) - self._check_sequence() - except Exception: - self.log.exception("Checks not passed") - raise - - gradient_index: Dimensions = parameter.channel_assignment - - # Get list of all events and list of unique RF and ADC events, since they are frequently reused - events_list = self.block_events - seq_duration, _, _ = self.duration() - seq_samples = round(seq_duration * self.spcm_freq) - - # Calculate the start time (and sample position) and duration of each block - block_durations = np.array( - [self.get_block(block_idx).block_duration for block_idx in list(events_list.keys())], - ) - block_durations = np.round(block_durations * self.spcm_freq).astype(int) - block_pos = np.cumsum(block_durations, dtype=np.int64) - block_pos = np.insert(block_pos, 0, 0) - - if seq_samples != block_pos[-1]: - msg = "Number of sequence samples does not match total number of block samples" - raise IndexError(msg) - - # Setup output arrays - _seq = np.zeros(4 * seq_samples, dtype=np.int16) - _rx_data = [] # list containing rx data objects for each ADC event - - # Count the total number of sample points and gate signals - adc_count: int = 0 - labels = {} - - for event_idx, (event_key, event) in enumerate(events_list.items()): - block = self.get_block(event_key) - # Calculate gradient waveform start and end positions according to block position - waveform_start = block_pos[event_idx] * 4 - - if block.gx is not None: # Gx event - waveform = self._calculate_gradient( - block=block.gx, - # FoV scaling refers to the sequence, block.gx -> x - fov_scaling=parameter.fov_scaling.x, - # Offset value are set independent of the sequence orientation, - # must be considered with respect to the target output channel! - # Offsets mapping: x -> channel 1, y -> channel 2, z -> channel 3 - offset=parameter.gradient_offset.to_list()[int(gradient_index.x-1)], - # Gradient indexing starts at 1 (RF is channel 0) - # -> correct indexing to match tuple index - output_channel=int(gradient_index.x-1), - ) - delay = block.gx.delay - delay_samples = round(delay * self.spcm_freq) - waveform_start_gx = waveform_start + 4 * delay_samples - gx_slice = slice( - waveform_start_gx + gradient_index.x, - waveform_start_gx + 4 * np.size(waveform) + gradient_index.x, - 4, - ) - _seq[gx_slice] = waveform - - if block.gy is not None: # Gy event - waveform = self._calculate_gradient( - block=block.gy, - # FoV scaling refers to the sequence, block.gy -> y - fov_scaling=parameter.fov_scaling.y, - # Offset value are set independent of the sequence orientation, - # must be considered with respect to the target output channel! - # Offsets mapping: x -> channel 1, y -> channel 2, z -> channel 3 - offset=parameter.gradient_offset.to_list()[int(gradient_index.y-1)], - # Gradient indexing starts at 1 (RF is channel 0) - # -> correct indexing to match tuple index - output_channel=int(gradient_index.y-1), - ) - delay = block.gy.delay - delay_samples = round(delay * self.spcm_freq) - waveform_start_gy = waveform_start + 4 * delay_samples - gy_slice = slice( - waveform_start_gy + gradient_index.y, - waveform_start_gy + 4 * np.size(waveform) + gradient_index.y, - 4, - ) - _seq[gy_slice] = waveform - - if block.gz is not None: # Gz event - waveform = self._calculate_gradient( - block=block.gz, - # FoV scaling refers to the sequence, block.gz -> z - fov_scaling=parameter.fov_scaling.z, - # Offset value are set independent of the sequence orientation, - # must be considered with respect to the target output channel! - # Offsets mapping: x -> channel 1, y -> channel 2, z -> channel 3 - offset=parameter.gradient_offset.to_list()[int(gradient_index.z-1)], - # Gradient indexing starts at 1 (RF is channel 0) - # -> correct indexing to match tuple index - output_channel=int(gradient_index.z-1), - ) - delay = block.gz.delay - delay_samples = round(delay * self.spcm_freq) - waveform_start_gz = waveform_start + 4 * delay_samples - gz_slice = slice( - waveform_start_gz + gradient_index.z, - waveform_start_gz + 4 * np.size(waveform) + gradient_index.z, - 4, - ) - _seq[gz_slice] = waveform - - if block.rf is not None: # RF event - # Pre-calculated RF event size can be shorter than the duration of the block since it doesn't - # consider the post-pulse ring-down time. The RF waveform is placed at the start of the block - # and the array is then sliced using the duration of the RF waveform to ensure a good fit - rf_waveform, rf_unblanking = self._calculate_rf( - block=block.rf, - b1_scaling=parameter.b1_scaling, - larmor_frequency=parameter.larmor_frequency, - ) - - rf_size = np.size(rf_waveform) # Get size of the RF waveform - if rf_size > (block_pos[event_idx + 1] - block_pos[event_idx]): - msg = "RF waveform size exceeds block size." - raise IndexError(msg) - - # Calculate RF waveform start and end positions according to block position - rf_start = block_pos[event_idx] * 4 - rf_end = (block_pos[event_idx] + rf_size) * 4 - - # Add RF waveform and unblanking signal to Z gradient - _seq[rf_start:rf_end:4] = rf_waveform - _seq[rf_start + 3:rf_end + 3:4] = _seq[rf_start + 3:rf_end + 3:4] | rf_unblanking - - if block.label is not None: - # Update dictionary with current labels - for label in block.label.values(): - labels[label.label] = label.value - - if block.adc is not None: # ADC event - # Calculate the number of samples to be discarded from the decimated signal - num_samples_discard = floor(block.adc.dead_time / block.adc.dwell) - # Calculate the total gate duration, given by number of samples - # and two times the number of discarded samples for symmetric adc dead time - # Note: The total gate duration is only increased if the dead time is a multiple of the adc dwell time. - total_gate_duration = (block.adc.num_samples + 2 * num_samples_discard) * block.adc.dwell - num_samples_raw = round(total_gate_duration * self.spcm_freq) - - # Remaining delay = dead_time minus pre- and post-sampling fractions - remaining_delay = block.adc.delay - num_samples_discard * block.adc.dwell - num_delay_samples = round(remaining_delay * self.spcm_freq) - - adc_start = (block_pos[event_idx] + num_delay_samples) * 4 - adc_end = adc_start + num_samples_raw * 4 - - # Add ADC gate to 16th bit of output channel 1 (first gradient channel) - _seq[adc_start + 1:adc_end + 1:4] |= np.uint16(2**15) - - # Add phase reference signal - num_samples_reference = min(num_samples_raw, self.phase_reference.size) - phase_ref_end = adc_start + num_samples_reference * 4 - _seq[adc_start + 2:phase_ref_end + 2:4] |= self.phase_reference[:num_samples_reference] - - _rx_data.append( - RxData( - index=adc_count, - num_samples=block.adc.num_samples, - num_samples_raw=num_samples_raw, - num_samples_discard=num_samples_discard, - dwell_time=block.adc.dwell, - dwell_time_raw=self.spcm_dwell_time, - phase_offset=block.adc.phase_offset, - freq_offset=block.adc.freq_offset, - total_averages=parameter.num_averages, - ddc_method=parameter.ddc_method, - phase_ref_frequency=REFERENCE_FREQUENCY, - labels=labels, - ) - ) - adc_count += 1 - labels = {} # Reset labels dict - - self.log.debug( - "Unrolled sequence; Total sample points: %s; Total block events: %s", - seq_samples, - len(block_durations), - ) - - return UnrolledSequence( - seq=_seq, - sample_count=seq_samples, - gpa_gain=self.gpa_gain, - gradient_efficiency=self.grad_eff, - rf_to_mvolt=self.rf_to_mvolt, - dwell_time=self.spcm_dwell_time, - gradient_output_limits=self.gradient_out_limits, - rf_output_limit=self.rf_out_limit, - duration=self.duration()[0], - adc_count=adc_count, - parameter=parameter, - rx_data=_rx_data, - ) - - # -------- Private waveform calculation functions -------- # - - @profile - def _calculate_rf( - self, - block: SimpleNamespace, - larmor_frequency: float, - b1_scaling: float, - ) -> tuple[np.ndarray, np.ndarray]: - """Calculate RF sample points to be played by TX card. - - Parameters - ---------- - block - Pulseq RF block - larmor_frequency - Larmor frequency of RF waveform - b1_scaling - Experiment dependent scaling factor of the RF amplitude - - Returns - ------- - List with the RF pulse in the first element, and the unblanking signal in the second element - - Raises - ------ - ValueError - Invalid RF block - """ - try: - if not block.type == "rf": - raise ValueError("Sequence block event is not a valid RF event.") - if not larmor_frequency > 0.: - raise ValueError(f"Invalid Larmor frequency: {larmor_frequency}") - except ValueError as err: - self.log.exception(err, exc_info=True) - raise err - - # Calculate the number of delay samples before an RF event (and unblanking) - # Note that the RF ring-down time is handled implicitly: the block duration used to place the RF waveform - # already includes the post-pulse dead time, so no additional handling is required. - # Dead-time is automatically set as delay! Delay accounts for start of RF event - num_samples_delay = round(max(block.dead_time, block.delay) * self.spcm_freq) - # Calculate the number of dead-time samples between unblanking and RF event - # Delay - dead-time samples account for start of unblanking - num_samples_dead_time = round(block.dead_time * self.spcm_freq) - # Calculate the number of RF shape sample points - num_samples = round(block.shape_dur * self.spcm_freq) - - # Set unblanking signal: 16th bit set to 1 (high) - rf_unblanking_start = num_samples_delay - num_samples_dead_time - rf_unblanking_end = num_samples_delay + num_samples - rf_unblanking = np.zeros(rf_unblanking_end, dtype=np.uint16) - rf_unblanking[rf_unblanking_start:] = 2**15 - - # Calculate the static phase offset, defined by RF pulse - phase_offset = np.exp(1j * block.phase_offset) - - # Calculate scaled envelope and convert to int16 scale (not datatype, since we use complex numbers) - # Perform this step here to save computation time, num. of envelope samples << num. of resampled signal - try: - # RF scaling according to B1 calibration and "device" (translation from pulseq to output voltage) - rf_scaling = b1_scaling * self.rf_to_mvolt * phase_offset / self.rf_out_limit - if np.abs(np.amax(envelope_scaled := block.signal * rf_scaling)) > 1: - raise ValueError(f"RF magnitude exceeds output limit by {np.amax(envelope_scaled)*100}%.") - except ValueError as err: - self.log.exception(err, exc_info=True) - raise err - - envelope_scaled = envelope_scaled * INT16_MAX - - # Resampling of scaled complex envelope - envelope = resample(envelope_scaled, num=num_samples) - - # Only precalculate carrier time array, calculate carrier here to take into account the - # frequency and phase offsets of an RF block event - carrier_time = np.arange(num_samples) * self.spcm_dwell_time - carrier = np.exp(2j * np.pi * (larmor_frequency + block.freq_offset) * carrier_time) - - try: - rf_waveform = np.concatenate((np.zeros(num_samples_delay, dtype=complex), (envelope * carrier))) - except IndexError as err: - self.log.exception(err, exc_info=True) - - rf_waveform_i16 = rf_waveform.real.astype(np.int16) - if self._rf_gain_lut is not None: - rf_waveform_i16 = self._rf_gain_lut[rf_waveform_i16+INT16_MIN] - - return (rf_waveform_i16, rf_unblanking) - - @profile - def _calculate_gradient( - self, - block: SimpleNamespace, - fov_scaling: float, - offset: float, - output_channel: int, - ) -> np.ndarray: - """Calculate spectrum-card sample points of a pypulseq gradient block event. - - Parameters - ---------- - block - Gradient block from pypulseq sequence, type must be grad or trap - unroll_arr - Section of numpy array which will contain the unrolled gradient event - fov_scaling - Scaling factor to adjust the FoV. - Factor is applied to the whole gradient waveform, exception the amplitude offset. - - Returns - ------- - Array with sample points of RF waveform as int16 values - - Raises - ------ - ValueError - Invalid block type (must be either ``grad`` or ``trap``), - gradient amplitude exceeds channel maximum output level - """ - try: - # Calculate gradient waveform scaling - scaling = fov_scaling / ( - self.system.gamma * 1e-3 * self.gpa_gain[output_channel] * self.grad_eff[output_channel] - ) - - # Calculate the gradient waveform relative to max output (within the interval [0, 1]) - if block.type == "grad": - # Arbitrary gradient waveform, interpolate linearly - # This function requires float input => cast to int16 afterwards - waveform = block.waveform * scaling / self.gradient_out_limits[output_channel] - self._check_gradient_amplitude(output_channel, np.abs(np.amax(waveform))) - # Transfer mV floating point waveform values to int16 if amplitude check passed - waveform_i16 = waveform * INT16_MAX - # Interpolate waveform on spectrum card time raster - gradient = np.interp( - x=np.linspace(block.tt[0], block.tt[-1], round(block.shape_dur / self.spcm_dwell_time)), - xp=block.tt, - fp=waveform_i16, - ) - - elif block.type == "trap": - # Construct trapezoidal gradient from rise, flat and fall sections - flat_amp = block.amplitude * scaling / self.gradient_out_limits[output_channel] - self._check_gradient_amplitude(output_channel, np.abs(np.amax(flat_amp))) - # Transfer relative floating point flat amplitude to int16 if amplitude check passed - flat_amp_i16 = flat_amp * INT16_MAX - # Define rise, flat and fall sections of trapezoidal gradient on spectrum card time raster - rise = np.linspace(0, flat_amp_i16, round(block.rise_time / self.spcm_dwell_time)) - flat = np.full(round(block.flat_time / self.spcm_dwell_time), fill_value=flat_amp_i16) - fall = np.linspace(flat_amp_i16, 0, round(block.fall_time / self.spcm_dwell_time)) - # Combine rise, flat and fall sections to gradient waveform - gradient = np.concatenate((rise, flat, fall)) - - else: - raise ValueError("Block is not a valid gradient block") - - # Calculate gradient offset int16 value from mV - # Gradient offset is used for calculating output limits but is not added to the waveform - offset_i16 = offset * INT16_MAX / self.gradient_out_limits[output_channel] - # This is the combined int16 gradient and offset waveform as float dtype - combined_i16 = gradient + offset_i16 - if (max_strength_i16 := np.amax(combined_i16)) > INT16_MAX: - # Report maximum strength in mV - max_strength = max_strength_i16 * self.gradient_out_limits[output_channel] / INT16_MAX - msg = f"Amplitude of combined gradient and shim waveforms {max_strength} exceed max gradient amplitude" - raise ValueError(msg) - - # Shifting gradient waveform to 15 bits already for adding the gate signals later - return gradient.astype(np.int16).view(np.uint16) >> 1 - - except (ValueError, IndexError): - self.log.exception("Error calculating gradient") - raise - - # -------- Private validation methods -------- # - - def _check_gradient_amplitude(self, idx: int, rel_value: float) -> None: - """Raise error if amplitude exceeds output limit.""" - limit = self.gradient_out_limits[idx] - if np.abs(rel_value) > 1.: - msg = f"Amplitude of gradient channel {idx+1} ({rel_value*limit}) exceeded output limit ({limit}))" - raise ValueError(msg) - - def _check_parameter(self, parameter: AcquisitionParameter) -> None: - """Check acquisition parameter and raise error if invalid.""" - # Check larmor frequency - f0_limit = self.spcm_freq / 2 - if parameter.larmor_frequency >= f0_limit: - msg = f"Larmor frequency too high ({parameter.larmor_frequency * 1e-6} MHz), violating sampling theorem" - raise ValueError(msg) - if parameter.larmor_frequency < 0: - msg = "Larmor frequency invalid (< 0)." - raise ValueError(msg) - - # Validate channel assignment - grad_ch: Dimensions = parameter.channel_assignment - if not all(isinstance(v, int) for v in (grad_ch.x, grad_ch.y, grad_ch.z)): - raise TypeError("All channel_assignment values must be integers.") - if {grad_ch.x, grad_ch.y, grad_ch.z} != {1, 2, 3}: - msg = f"Invalid channel assignment, must contain each of 1, 2, and 3 exactly once, got: {grad_ch}" - raise ValueError(msg) - - def _check_sequence(self) -> None: - """Check sequence.""" - # Check number of block events - if not len(self.block_events) > 0: - raise ValueError("No block events found") - # Perform sequence timing check - check, seq_err = self.check_timing() - if not check: - raise ValueError(f"Sequence timing check failed: {seq_err}") - - def _load_rfpa_lut(self, rf_gain_lut_path: Path | None) -> None: - if rf_gain_lut_path is None or not rf_gain_lut_path.exists(): - self.log.info("No RFPA LUT file.") - return - - _lut = np.load(rf_gain_lut_path) - _required_lut_size = INT16_MAX - INT16_MIN + 1 - - if _lut.size != _required_lut_size: - self.log.warning( - f"Error loading RFPA LUT: Invalid size.\nLoaded: {_lut.size}, required: {_required_lut_size}" - ) - return - if _lut.dtype != np.int16: - self.log.warning(f"Invalid data type of RFPA gain LUT: {_lut.dtype}") - return - - self.log.info("Successfully loaded LUT for RF gain correction.") - self._rf_gain_lut = _lut +"""Sequence provider class.""" +import logging +from collections.abc import Callable +from dataclasses import dataclass +from math import floor +from pathlib import Path +from types import SimpleNamespace +from typing import Any + +import numpy as np +from pypulseq.opts import Opts +from pypulseq.Sequence.sequence import Sequence +from scipy.signal import resample + +from console.interfaces.acquisition_parameter import AcquisitionParameter +from console.interfaces.dimensions import Dimensions +from console.interfaces.rx_data import RxData +from console.interfaces.unrolled_sequence import UnrolledSequence + +try: + from line_profiler import profile +except ImportError: + def profile(func: Callable[..., Any]) -> Callable[..., Any]: + """Define placeholder for profile decorator.""" + return func + + +INT16_MAX = np.iinfo(np.int16).max +INT16_MIN = np.iinfo(np.int16).min + +NUM_REFERENCE_SAMPLES = 1000 +REFERENCE_FREQUENCY = 1.095e6 + + +@dataclass +class ADCGate: + """Define precalculated attributes of an ADC gate.""" + + start: int + num_samples_discard: int + num_samples_raw: int + +class SequenceProvider(Sequence): + """Sequence provider class. + + This object is inherited from pulseq sequence object, so that all methods of the + pypulseq ``Sequence`` object can be accessed. + + The main functionality of the ``SequenceProvider`` is to unroll a given pulseq sequence. + Usually the first step is to read a sequence file. The unrolling step can be achieved using + the ``unroll_sequence()`` function. + + Example + ------- + >>> seq_provider = SequenceProvider() + >>> seq_provider.read("./seq_file.seq") + >>> unrolled = seq_provider.unroll_sequence(acquisition_parameter) + """ + + __name__: str = "SequenceProvider" + + def __init__( + self, + gradient_efficiency: tuple[float, float, float], + gpa_gain: tuple[float, float, float], + gradient_output_limits: tuple[int, int, int], + gradients_50ohms: bool, + rf_output_limit: int, + rf_50ohms: bool, + rf_to_mvolt: float, + spcm_dwell_time: float, + system: Opts, + rf_gain_lut_path: Path | None = None, + ): + """Initialize sequence provider class which is used to unroll a pulseq sequence. + + Parameters + ---------- + gradient_efficiency + Efficiency of the gradient coils in mT/m/A, e.g. [0.4e-3, 0.4e-3, 0.4e-3]. + gpa_gain + Gain factor of the GPA per gradient channel, e.g. [4.7, 4.7, 4.7]. + gradient_output_limits + Integer output limit per gradient channel in mV, e.g. [6000, 6000, 6000]. + gradients_50ohms + Boolean flag which indicates if the gradient output is terminated into 50 ohms or high impedance. + If terminated into high impedance, the card output doubles, + what needs to be considered when calculating the sequence. + rf_output_limit + Integer output limit of the RF channel in mV. + rf_50ohms + Boolean flag which indicates if the rf output is terminated into 50 ohms (see gradients_50ohms). + rf_to_mvolt + Translation of RF waveform from pulseq (Hz) to mV. + spcm_dwell_time + Sampling time raster of the output waveform (depends on spectrum card). + system_limits + Absolute maximum system limits defined in the device configuration. + Used to instantiate the pypulseq `Opts()` class. + """ + if not isinstance(system, Opts): + raise AttributeError("Invalid system: Pypulseq `Opts` definition required.") + super().__init__(system=system) + self.log = logging.getLogger("SeqProv") + + # Set class instance attributes + self.rf_to_mvolt = rf_to_mvolt + self.spcm_dwell_time = spcm_dwell_time + self.spcm_freq = 1 / spcm_dwell_time + self.gpa_gain = gpa_gain + self.grad_eff = gradient_efficiency + + # Scale output limit dependent on high impedance flags: + # If output is terminated into high impedance (flag is true), the channel output is doubled. + # Otherwise, if output is terminated into 50 ohms impedance, output limit remains unchanged. + self.rf_out_limit = rf_output_limit if rf_50ohms else int(2 * rf_output_limit) + # Ensure tuple[int, int, int] for mypy typing + self.gradient_out_limits = ( + gradient_output_limits[0] if gradients_50ohms else int(2 * gradient_output_limits[0]), + gradient_output_limits[1] if gradients_50ohms else int(2 * gradient_output_limits[1]), + gradient_output_limits[2] if gradients_50ohms else int(2 * gradient_output_limits[2]), + ) + + # Setup phase reference signal + time = np.arange(NUM_REFERENCE_SAMPLES) * self.spcm_dwell_time + signal = np.exp(2j * np.pi * REFERENCE_FREQUENCY * time) + self.phase_reference = np.zeros(NUM_REFERENCE_SAMPLES, dtype=np.uint16) + self.phase_reference[signal > 0] = np.uint16(2**15) + + # Load LUT for RF gain correction if path to LUT is provided + self._rf_gain_lut: np.ndarray | None = None + self._load_rfpa_lut(rf_gain_lut_path) + + + # -------- PyPulseq interface -------- # + + def from_pypulseq(self, seq: Sequence) -> None: + """Read a pypulseq sequence to sequence provider. + + Parameters + ---------- + seq + Pypulseq ``Sequence`` instance + + Raises + ------ + AttributeError + seq is not a valid pypulseq ``Sequence`` instance + """ + if not isinstance(seq, Sequence): + raise AttributeError("Invalid sequence.") + # Re-initialize the parent to start from a clean pypulseq sequence + super().__init__(system=self.system) + for block_index, _ in seq.block_events.items(): + block = seq.get_block(block_index) + self.add_block(block) + # Set definitions without overwriting existing ones + for key, value in seq.definitions.items(): + self.definitions.setdefault(key, value) + + def to_pypulseq(self) -> Sequence | None: + """Create a pypulseq sequence from sequence provider.""" + seq = Sequence(system=self.system) + for block_index, _ in self.block_events.items(): + block = self.get_block(block_index) + seq.add_block(block) + seq.definitions = self.definitions + return seq + + # -------- Public interface -------- # + + def dict(self) -> dict: + """Abstract method which returns variables for logging in dictionary.""" + return { + "rf_to_mvolt": self.rf_to_mvolt, + "spcm_freq": self.spcm_freq, + "spcm_dwell_time": self.spcm_dwell_time, + "gpa_gain": self.gpa_gain, + "gradient_efficiency": self.grad_eff, + "output_limits": self.gradient_out_limits, + } + + @profile + def unroll_sequence(self, parameter: AcquisitionParameter) -> UnrolledSequence: + """Unroll the pypulseq sequence description. + + TODO: Reduce complexity. + + Parameters + ---------- + parameter + Instance of AcquisitionParameter containing all necessary parameters to + calculate the sequence waveforms, i.e. larmor frequency, gradient offsets, etc. + + Returns + ------- + UnrolledSequence + Instance of an unrolled sequence object which contains a list of numpy arrays with + the block-wise calculated sample points in correct spectrum card order (Fortran). + + The list of unrolled sequence arrays is returned as uint16 values which contain a digital + signal encoded by 15th bit. Only the RF channel does not contain a digital signal. + In addition, all receive events are described and returned in a list within the unrolled + sequence object. + + Examples + -------- + For channels ch0, ch1, ch2, ch3, data values n = 0, 1, ..., N are ordered the following way. + + >>> data = [ch0_0, ch1_0, ch2_0, ch3_0, ch0_1, ch1_1, ..., ch0_n, ..., ch3_N] + + Per channel data can be extracted by the following code. + + >>> rf = seq[0::4] + >>> gx = (seq[1::4] << 1).astype(np.int16) + >>> gy = (seq[2::4] << 1).astype(np.int16) + >>> gz = (seq[3::4] << 1).astype(np.int16) + + All the gradient channels contain a digital signal encoded by the 15th bit. + - `gx`: ADC gate signal + - `gy`: Reference signal for phase correction + - `gz`: RF unblanking signal + The following example shows, how to extract the digital signals + + >>> adc_gate = seq[1::4].astype(np.uint16) >> 15 + >>> reference = seq[2::4].astype(np.uint16) >> 15 + >>> unblanking = seq[3::4].astype(np.uint16) >> 15 + + As the 15th bit is not encoding the sign (as usual for int16), the values are casted to uint16 before shifting. + """ + try: + self._check_parameter(parameter) + self._check_sequence() + except Exception: + self.log.exception("Checks not passed") + raise + + gradient_index: Dimensions = parameter.channel_assignment + + # Get list of all events and list of unique RF and ADC events, since they are frequently reused + events_list = self.block_events + seq_duration, _, _ = self.duration() + seq_samples = round(seq_duration * self.spcm_freq) + + # Calculate the start time (and sample position) and duration of each block + block_durations = np.array( + [self.get_block(block_idx).block_duration for block_idx in list(events_list.keys())], + ) + block_durations = np.round(block_durations * self.spcm_freq).astype(int) + block_pos = np.cumsum(block_durations, dtype=np.int64) + block_pos = np.insert(block_pos, 0, 0) + + if seq_samples != block_pos[-1]: + msg = "Number of sequence samples does not match total number of block samples" + raise IndexError(msg) + + # Setup output arrays + _seq = np.zeros(4 * seq_samples, dtype=np.int16) + _rx_data = [] # list containing rx data objects for each ADC event + + # Count the total number of sample points and gate signals + adc_count: int = 0 + labels = {} + + for event_idx, (event_key, event) in enumerate(events_list.items()): + block = self.get_block(event_key) + # Calculate gradient waveform start and end positions according to block position + waveform_start = block_pos[event_idx] * 4 + + if block.gx is not None: # Gx event + waveform = self._calculate_gradient( + block=block.gx, + # FoV scaling refers to the sequence, block.gx -> x + fov_scaling=parameter.fov_scaling.x, + # Offset value are set independent of the sequence orientation, + # must be considered with respect to the target output channel! + # Offsets mapping: x -> channel 1, y -> channel 2, z -> channel 3 + offset=parameter.gradient_offset.to_list()[int(gradient_index.x-1)], + # Gradient indexing starts at 1 (RF is channel 0) + # -> correct indexing to match tuple index + output_channel=int(gradient_index.x-1), + ) + delay = block.gx.delay + delay_samples = round(delay * self.spcm_freq) + waveform_start_gx = waveform_start + 4 * delay_samples + gx_slice = slice( + waveform_start_gx + gradient_index.x, + waveform_start_gx + 4 * np.size(waveform) + gradient_index.x, + 4, + ) + _seq[gx_slice] = waveform + + if block.gy is not None: # Gy event + waveform = self._calculate_gradient( + block=block.gy, + # FoV scaling refers to the sequence, block.gy -> y + fov_scaling=parameter.fov_scaling.y, + # Offset value are set independent of the sequence orientation, + # must be considered with respect to the target output channel! + # Offsets mapping: x -> channel 1, y -> channel 2, z -> channel 3 + offset=parameter.gradient_offset.to_list()[int(gradient_index.y-1)], + # Gradient indexing starts at 1 (RF is channel 0) + # -> correct indexing to match tuple index + output_channel=int(gradient_index.y-1), + ) + delay = block.gy.delay + delay_samples = round(delay * self.spcm_freq) + waveform_start_gy = waveform_start + 4 * delay_samples + gy_slice = slice( + waveform_start_gy + gradient_index.y, + waveform_start_gy + 4 * np.size(waveform) + gradient_index.y, + 4, + ) + _seq[gy_slice] = waveform + + if block.gz is not None: # Gz event + waveform = self._calculate_gradient( + block=block.gz, + # FoV scaling refers to the sequence, block.gz -> z + fov_scaling=parameter.fov_scaling.z, + # Offset value are set independent of the sequence orientation, + # must be considered with respect to the target output channel! + # Offsets mapping: x -> channel 1, y -> channel 2, z -> channel 3 + offset=parameter.gradient_offset.to_list()[int(gradient_index.z-1)], + # Gradient indexing starts at 1 (RF is channel 0) + # -> correct indexing to match tuple index + output_channel=int(gradient_index.z-1), + ) + delay = block.gz.delay + delay_samples = round(delay * self.spcm_freq) + waveform_start_gz = waveform_start + 4 * delay_samples + gz_slice = slice( + waveform_start_gz + gradient_index.z, + waveform_start_gz + 4 * np.size(waveform) + gradient_index.z, + 4, + ) + _seq[gz_slice] = waveform + + if block.rf is not None: # RF event + # Pre-calculated RF event size can be shorter than the duration of the block since it doesn't + # consider the post-pulse ring-down time. The RF waveform is placed at the start of the block + # and the array is then sliced using the duration of the RF waveform to ensure a good fit + rf_waveform, rf_unblanking = self._calculate_rf( + block=block.rf, + b1_scaling=parameter.b1_scaling, + larmor_frequency=parameter.larmor_frequency, + ) + + rf_size = np.size(rf_waveform) # Get size of the RF waveform + if rf_size > (block_pos[event_idx + 1] - block_pos[event_idx]): + msg = "RF waveform size exceeds block size." + raise IndexError(msg) + + # Calculate RF waveform start and end positions according to block position + rf_start = block_pos[event_idx] * 4 + rf_end = (block_pos[event_idx] + rf_size) * 4 + + # Add RF waveform and unblanking signal to Z gradient + _seq[rf_start:rf_end:4] = rf_waveform + _seq[rf_start + 3:rf_end + 3:4] = _seq[rf_start + 3:rf_end + 3:4] | rf_unblanking + + if block.label is not None: + # Update dictionary with current labels + for label in block.label.values(): + labels[label.label] = label.value + + if block.adc is not None: # ADC event + # Calculate the number of samples to be discarded from the decimated signal + num_samples_discard = floor(block.adc.dead_time / block.adc.dwell) + # Calculate the total gate duration, given by number of samples + # and two times the number of discarded samples for symmetric adc dead time + # Note: The total gate duration is only increased if the dead time is a multiple of the adc dwell time. + total_gate_duration = (block.adc.num_samples + 2 * num_samples_discard) * block.adc.dwell + num_samples_raw = round(total_gate_duration * self.spcm_freq) + + # Remaining delay = dead_time minus pre- and post-sampling fractions + remaining_delay = block.adc.delay - num_samples_discard * block.adc.dwell + num_delay_samples = round(remaining_delay * self.spcm_freq) + + adc_start = (block_pos[event_idx] + num_delay_samples) * 4 + adc_end = adc_start + num_samples_raw * 4 + + # Add ADC gate to 16th bit of output channel 1 (first gradient channel) + _seq[adc_start + 1:adc_end + 1:4] |= np.uint16(2**15) + + # Add phase reference signal + num_samples_reference = min(num_samples_raw, self.phase_reference.size) + phase_ref_end = adc_start + num_samples_reference * 4 + _seq[adc_start + 2:phase_ref_end + 2:4] |= self.phase_reference[:num_samples_reference] + + _rx_data.append( + RxData( + index=adc_count, + num_samples=block.adc.num_samples, + num_samples_raw=num_samples_raw, + num_samples_discard=num_samples_discard, + dwell_time=block.adc.dwell, + dwell_time_raw=self.spcm_dwell_time, + phase_offset=block.adc.phase_offset, + freq_offset=block.adc.freq_offset, + total_averages=parameter.num_averages, + ddc_method=parameter.ddc_method, + phase_ref_frequency=REFERENCE_FREQUENCY, + labels=labels, + ) + ) + adc_count += 1 + labels = {} # Reset labels dict + + self.log.debug( + "Unrolled sequence; Total sample points: %s; Total block events: %s", + seq_samples, + len(block_durations), + ) + + return UnrolledSequence( + seq=_seq, + sample_count=seq_samples, + gpa_gain=self.gpa_gain, + gradient_efficiency=self.grad_eff, + rf_to_mvolt=self.rf_to_mvolt, + dwell_time=self.spcm_dwell_time, + gradient_output_limits=self.gradient_out_limits, + rf_output_limit=self.rf_out_limit, + duration=self.duration()[0], + adc_count=adc_count, + parameter=parameter, + rx_data=_rx_data, + ) + + # -------- Private waveform calculation functions -------- # + + @profile + def _calculate_rf( + self, + block: SimpleNamespace, + larmor_frequency: float, + b1_scaling: float, + ) -> tuple[np.ndarray, np.ndarray]: + """Calculate RF sample points to be played by TX card. + + Parameters + ---------- + block + Pulseq RF block + larmor_frequency + Larmor frequency of RF waveform + b1_scaling + Experiment dependent scaling factor of the RF amplitude + + Returns + ------- + List with the RF pulse in the first element, and the unblanking signal in the second element + + Raises + ------ + ValueError + Invalid RF block + """ + try: + if not block.type == "rf": + raise ValueError("Sequence block event is not a valid RF event.") + if not larmor_frequency > 0.: + raise ValueError(f"Invalid Larmor frequency: {larmor_frequency}") + except ValueError as err: + self.log.exception(err, exc_info=True) + raise err + + # Calculate the number of delay samples before an RF event (and unblanking) + # Note that the RF ring-down time is handled implicitly: the block duration used to place the RF waveform + # already includes the post-pulse dead time, so no additional handling is required. + # Dead-time is automatically set as delay! Delay accounts for start of RF event + num_samples_delay = round(max(block.dead_time, block.delay) * self.spcm_freq) + # Calculate the number of dead-time samples between unblanking and RF event + # Delay - dead-time samples account for start of unblanking + num_samples_dead_time = round(block.dead_time * self.spcm_freq) + # Calculate the number of RF shape sample points + num_samples = round(block.shape_dur * self.spcm_freq) + + # Set unblanking signal: 16th bit set to 1 (high) + rf_unblanking_start = num_samples_delay - num_samples_dead_time + rf_unblanking_end = num_samples_delay + num_samples + rf_unblanking = np.zeros(rf_unblanking_end, dtype=np.uint16) + rf_unblanking[rf_unblanking_start:] = 2**15 + + # Calculate the static phase offset, defined by RF pulse + phase_offset = np.exp(1j * block.phase_offset) + + # Calculate scaled envelope and convert to int16 scale (not datatype, since we use complex numbers) + # Perform this step here to save computation time, num. of envelope samples << num. of resampled signal + try: + # RF scaling according to B1 calibration and "device" (translation from pulseq to output voltage) + rf_scaling = b1_scaling * self.rf_to_mvolt * phase_offset / self.rf_out_limit + if np.abs(np.amax(envelope_scaled := block.signal * rf_scaling)) > 1: + raise ValueError(f"RF magnitude exceeds output limit by {np.amax(envelope_scaled)*100}%.") + except ValueError as err: + self.log.exception(err, exc_info=True) + raise err + + envelope_scaled = envelope_scaled * INT16_MAX + + # Resampling of scaled complex envelope + envelope = resample(envelope_scaled, num=num_samples) + + # Only precalculate carrier time array, calculate carrier here to take into account the + # frequency and phase offsets of an RF block event + carrier_time = np.arange(num_samples) * self.spcm_dwell_time + carrier = np.exp(2j * np.pi * (larmor_frequency + block.freq_offset) * carrier_time) + + try: + rf_waveform = np.concatenate((np.zeros(num_samples_delay, dtype=complex), (envelope * carrier))) + except IndexError as err: + self.log.exception(err, exc_info=True) + + rf_waveform_i16 = rf_waveform.real.astype(np.int16) + if self._rf_gain_lut is not None: + rf_waveform_i16 = self._rf_gain_lut[rf_waveform_i16+INT16_MIN] + + return (rf_waveform_i16, rf_unblanking) + + @profile + def _calculate_gradient( + self, + block: SimpleNamespace, + fov_scaling: float, + offset: float, + output_channel: int, + ) -> np.ndarray: + """Calculate spectrum-card sample points of a pypulseq gradient block event. + + Parameters + ---------- + block + Gradient block from pypulseq sequence, type must be grad or trap + unroll_arr + Section of numpy array which will contain the unrolled gradient event + fov_scaling + Scaling factor to adjust the FoV. + Factor is applied to the whole gradient waveform, exception the amplitude offset. + + Returns + ------- + Array with sample points of RF waveform as int16 values + + Raises + ------ + ValueError + Invalid block type (must be either ``grad`` or ``trap``), + gradient amplitude exceeds channel maximum output level + """ + try: + # Calculate gradient waveform scaling + scaling = fov_scaling / ( + self.system.gamma * 1e-3 * self.gpa_gain[output_channel] * self.grad_eff[output_channel] + ) + + # Calculate the gradient waveform relative to max output (within the interval [0, 1]) + if block.type == "grad": + # Arbitrary gradient waveform, interpolate linearly + # This function requires float input => cast to int16 afterwards + waveform = block.waveform * scaling / self.gradient_out_limits[output_channel] + self._check_gradient_amplitude(output_channel, np.abs(np.amax(waveform))) + # Transfer mV floating point waveform values to int16 if amplitude check passed + waveform_i16 = waveform * INT16_MAX + # Interpolate waveform on spectrum card time raster + gradient = np.interp( + x=np.linspace(block.tt[0], block.tt[-1], round(block.shape_dur / self.spcm_dwell_time)), + xp=block.tt, + fp=waveform_i16, + ) + + elif block.type == "trap": + # Construct trapezoidal gradient from rise, flat and fall sections + flat_amp = block.amplitude * scaling / self.gradient_out_limits[output_channel] + self._check_gradient_amplitude(output_channel, np.abs(np.amax(flat_amp))) + # Transfer relative floating point flat amplitude to int16 if amplitude check passed + flat_amp_i16 = flat_amp * INT16_MAX + # Define rise, flat and fall sections of trapezoidal gradient on spectrum card time raster + rise = np.linspace(0, flat_amp_i16, round(block.rise_time / self.spcm_dwell_time)) + flat = np.full(round(block.flat_time / self.spcm_dwell_time), fill_value=flat_amp_i16) + fall = np.linspace(flat_amp_i16, 0, round(block.fall_time / self.spcm_dwell_time)) + # Combine rise, flat and fall sections to gradient waveform + gradient = np.concatenate((rise, flat, fall)) + + else: + raise ValueError("Block is not a valid gradient block") + + # Calculate gradient offset int16 value from mV + # Gradient offset is used for calculating output limits but is not added to the waveform + offset_i16 = offset * INT16_MAX / self.gradient_out_limits[output_channel] + # This is the combined int16 gradient and offset waveform as float dtype + combined_i16 = gradient + offset_i16 + if (max_strength_i16 := np.amax(combined_i16)) > INT16_MAX: + # Report maximum strength in mV + max_strength = max_strength_i16 * self.gradient_out_limits[output_channel] / INT16_MAX + msg = f"Amplitude of combined gradient and shim waveforms {max_strength} exceed max gradient amplitude" + raise ValueError(msg) + + # Shifting gradient waveform to 15 bits already for adding the gate signals later + return gradient.astype(np.int16).view(np.uint16) >> 1 + + except (ValueError, IndexError): + self.log.exception("Error calculating gradient") + raise + + # -------- Private validation methods -------- # + + def _check_gradient_amplitude(self, idx: int, rel_value: float) -> None: + """Raise error if amplitude exceeds output limit.""" + limit = self.gradient_out_limits[idx] + if np.abs(rel_value) > 1.: + msg = f"Amplitude of gradient channel {idx+1} ({rel_value*limit}) exceeded output limit ({limit}))" + raise ValueError(msg) + + def _check_parameter(self, parameter: AcquisitionParameter) -> None: + """Check acquisition parameter and raise error if invalid.""" + # Check larmor frequency + f0_limit = self.spcm_freq / 2 + if parameter.larmor_frequency >= f0_limit: + msg = f"Larmor frequency too high ({parameter.larmor_frequency * 1e-6} MHz), violating sampling theorem" + raise ValueError(msg) + if parameter.larmor_frequency < 0: + msg = "Larmor frequency invalid (< 0)." + raise ValueError(msg) + + # Validate channel assignment + grad_ch: Dimensions = parameter.channel_assignment + if not all(isinstance(v, int) for v in (grad_ch.x, grad_ch.y, grad_ch.z)): + raise TypeError("All channel_assignment values must be integers.") + if {grad_ch.x, grad_ch.y, grad_ch.z} != {1, 2, 3}: + msg = f"Invalid channel assignment, must contain each of 1, 2, and 3 exactly once, got: {grad_ch}" + raise ValueError(msg) + + def _check_sequence(self) -> None: + """Check sequence.""" + # Check number of block events + if not len(self.block_events) > 0: + raise ValueError("No block events found") + # Perform sequence timing check + check, seq_err = self.check_timing() + if not check: + raise ValueError(f"Sequence timing check failed: {seq_err}") + + def _load_rfpa_lut(self, rf_gain_lut_path: Path | None) -> None: + if rf_gain_lut_path is None or not rf_gain_lut_path.exists(): + self.log.info("No RFPA LUT file.") + return + + _lut = np.load(rf_gain_lut_path) + _required_lut_size = INT16_MAX - INT16_MIN + 1 + + if _lut.size != _required_lut_size: + self.log.warning( + f"Error loading RFPA LUT: Invalid size.\nLoaded: {_lut.size}, required: {_required_lut_size}" + ) + return + if _lut.dtype != np.int16: + self.log.warning(f"Invalid data type of RFPA gain LUT: {_lut.dtype}") + return + + self.log.info("Successfully loaded LUT for RF gain correction.") + self._rf_gain_lut = _lut diff --git a/src/console/spcm_control/rx_device.py b/src/console/spcm_control/rx_device.py index 527d3f4e..7845b169 100644 --- a/src/console/spcm_control/rx_device.py +++ b/src/console/spcm_control/rx_device.py @@ -1,533 +1,533 @@ -"""Implementation of receive card.""" - -import logging -import threading -import time -from collections.abc import Callable -from ctypes import POINTER, addressof, byref, c_short, cast -from dataclasses import dataclass -from itertools import compress - -import numpy as np - -import console.spcm_control.spcm.pyspcm as sp -from console.interfaces.rx_data import RxData -from console.pulseq_interpreter.sequence_provider import NUM_REFERENCE_SAMPLES -from console.spcm_control.abstract_device import SpectrumDevice -from console.spcm_control.spcm.tools import create_dma_buffer, type_to_name - -# Define registers lists -CH_SELECT = [ - sp.CHANNEL0, - sp.CHANNEL1, - sp.CHANNEL2, - sp.CHANNEL3, - sp.CHANNEL4, - sp.CHANNEL5, - sp.CHANNEL6, - sp.CHANNEL7, -] -AMP_SELECT = [ - sp.SPC_AMP0, - sp.SPC_AMP1, - sp.SPC_AMP2, - sp.SPC_AMP3, - sp.SPC_AMP4, - sp.SPC_AMP5, - sp.SPC_AMP6, - sp.SPC_AMP7, -] -IMP_SELECT = [ - sp.SPC_50OHM0, - sp.SPC_50OHM1, - sp.SPC_50OHM2, - sp.SPC_50OHM3, - sp.SPC_50OHM4, - sp.SPC_50OHM5, - sp.SPC_50OHM6, - sp.SPC_50OHM7, -] - - -@dataclass -class RxCard(SpectrumDevice): - """Implementation of RX device.""" - - __name__: str = "RxCard" - - def __init__( - self, - path: str, - sample_rate: int, - channel_enable: tuple[bool, ...], - max_amplitude: tuple[int, ...], - impedance_50_ohms: tuple[bool, ...], - ) -> None: - """Execute after init function to do further class setup.""" - self.log = logging.getLogger(self.__name__) - super().__init__(path, log=self.log) - - self.sample_rate = sample_rate - self.channel_enable = [int(val) for val in channel_enable] - self.max_amplitude = max_amplitude - self.impedance_50_ohms = [int(val) for val in impedance_50_ohms] - self.rx_data: list[RxData | None] | None = None - - self.num_channels = sp.int32(0) - self.card_type = sp.int32(0) - - self.worker: threading.Thread | None = None - self.is_running = threading.Event() - self.is_receiving = threading.Event() - self._total_gates: int = 0 - - # Pre trigger is set to minimum, post trigger depends on active channel count and is defined later. - self.pre_trigger: int = 8 - self.post_trigger: None | int = None - - self.rx_scaling = [amp / (2**15) for amp in self.max_amplitude] - - self._submit_fn: Callable[[int, RxData], None] | None = None - self._index_offset: int = 0 - - @property - def total_gates(self) -> int: - """Helper function to return the number of gates that have been collected by the Rx Card.""" - return self._total_gates - - def setup_card(self): - """Set up spectrum card in transmit (Rx) mode. - - At the very beginning, a card reset is performed. The clock mode is set according to the sample rate, - defined by the class attribute. - Two receive channels are enables and configured by max. amplitude according to class variables and impedance. - - Raises - ------ - Warning - The actual set sample rate deviates from the corresponding class attribute to be set, - class attribute is overwritten. - """ - # Get the card type and reset card - sp.spcm_dwGetParam_i32(self.card, sp.SPC_PCITYP, byref(self.card_type)) - sp.spcm_dwSetParam_i64(self.card, sp.SPC_M2CMD, sp.M2CMD_CARD_RESET) # Needed? - - try: - if "M2p.59" not in (device_type := type_to_name(self.card_type.value)): - raise ConnectionError("Device with path %s is of type %s, no receive card" % (self.path, device_type)) - except ConnectionError as err: - self.log.exception(err, exc_info=True) - raise err - - # Setup the internal clockmode, clock output enable (use RX clock output to enable anti-alias filter) - sp.spcm_dwSetParam_i32(self.card, sp.SPC_CLOCKMODE, sp.SPC_CM_INTPLL) - sp.spcm_dwSetParam_i32(self.card, sp.SPC_CLOCKOUT, 1) - - # Use external clock: Terminate to 50 Ohms, set threshold to 1.5V, suitable for 3.3V clock - # sp.spcm_dwSetParam_i32(self.card, sp.SPC_CLOCKMODE, sp.SPC_CM_EXTERNAL) - # sp.spcm_dwSetParam_i32(self.card, sp.SPC_CLOCK50OHM, 1) - # sp.spcm_dwSetParam_i32(self.card, sp.SPC_CLOCK_THRESHOLD, 1500) - - # Set card sampling rate in MHz and read the actual sampling rate - sp.spcm_dwSetParam_i64(self.card, sp.SPC_SAMPLERATE, sp.MEGA(self.sample_rate)) - sample_rate = sp.int64(0) - sp.spcm_dwGetParam_i64(self.card, sp.SPC_SAMPLERATE, byref(sample_rate)) - self.log.info("Device sampling rate: %s MHz", sample_rate.value * 1e-6) - - if sample_rate.value != sp.MEGA(self.sample_rate): - self.log.warning( - "Actual device sample rate %s MHz does not match set sample rate of %s MHz; Updating class attribute", - sample_rate.value * 1e-6, - self.sample_rate, - ) - self.sample_rate = int(sample_rate.value * 1e-6) - - # Check channel enable, max. amplitude per channel and impedance values - try: - # Check that the length of the channel enable list is 8 - # this has to be true for cards with fewer channels too - if (num_enable := len(self.channel_enable)) != 8: - raise ValueError("Channel enable list is incomplete: %s/8" % num_enable) - # Impedance and amplitude configuration lists must also be of length 8 - if (num_imp := len(self.impedance_50_ohms)) != 8: - raise ValueError("Channel impedance list is incomplete: %s/8" % num_imp) - if (num_amp := len(self.max_amplitude)) != 8: - raise ValueError("channel max. amplitude list is incomplete: %s/8" % num_amp) - # Number of enabled channels must be either 1, 2, 4 or 8 - if not np.log2(sum(self.channel_enable)).is_integer(): - raise ValueError("Invalid number of enabled channels, must be power of 2.") - except ValueError as err: - self.log.exception(err, exc_info=True) - raise err - - # Enable receive channels, compress list of channel select registers to obtain list of channels to be enabled - # Sum of the compressed list equals logical or operator - # e.g. sp.CHANNEL0 | sp.CHANNEL1 | sp.CHANNEL5 = sum([sp.CHANNEL0, sp.CHANNEL1, sp.CHANNEL5]) = 35 - channel_selection = sum(list(compress(CH_SELECT, map(bool, self.channel_enable)))) - sp.spcm_dwSetParam_i32(self.card, sp.SPC_CHENABLE, channel_selection) - - # Set impedance and amplitude limits for each channel according to device configuration - for k, enable in enumerate(map(bool, self.channel_enable)): - if enable: - self.log.info( - "Channel %s enabled; 50 ohms impedance: %s; Max. amplitude: %s mV", - k, - self.impedance_50_ohms[k], - self.max_amplitude[k], - ) - sp.spcm_dwSetParam_i32(self.card, IMP_SELECT[k], self.impedance_50_ohms[k]) - sp.spcm_dwSetParam_i32(self.card, AMP_SELECT[k], self.max_amplitude[k]) - - # Get the number of actual active channels and compare against provided channel enable list - sp.spcm_dwGetParam_i32(self.card, sp.SPC_CHCOUNT, byref(self.num_channels)) - try: - self.log.info( - "Number of enabled receive channels (read from card): %s", - self.num_channels.value, - ) - if not self.num_channels.value == sum(self.channel_enable): - raise ValueError("Actual number of enabled channels does not match the provided channel enable list") - except ValueError as err: - self.log.exception(err, exc_info=True) - raise err - - # Digital filter setting for receiver, 0 = disable digital bandwidth filter - sp.spcm_dwSetParam_i32(self.card, sp.SPC_DIGITALBWFILTER, 0) - - # Configure X2 as digital input for phase reference signal and sample it in sync with analog channel 0 - sp.spcm_dwSetParam_i32(self.card, sp.SPCM_X2_MODE, sp.SPCM_XMODE_DIGIN) - sp.spcm_dwSetParam_i32(self.card, sp.SPC_DIGMODE0, (sp.DIGMODEMASK_BIT15 & sp.SPCM_DIGMODE_X2)) - - # Calculate trigger size depending on the number of active channels - # Since data can only be gathered in notify size chunks, post_trigger // channel_count should be at least one - # notify size to ensure that we can always access the full gate data. - self.post_trigger = 4096 // self.num_channels.value - - # Set the memory size, pre and post trigger and loop parameters, SPC_LOOPS = 0 => runs infinitely long - sp.spcm_dwSetParam_i32(self.card, sp.SPC_POSTTRIGGER, self.post_trigger) - sp.spcm_dwSetParam_i32(self.card, sp.SPC_PRETRIGGER, self.pre_trigger) - sp.spcm_dwSetParam_i32(self.card, sp.SPC_LOOPS, 0) - - # Setup timestamp mode to read number of samples per gate if available - sp.spcm_dwSetParam_i32( - self.card, - sp.SPC_TIMESTAMP_CMD, - sp.SPC_TSMODE_STARTRESET | sp.SPC_TSCNT_INTERNAL, - ) - # Configure trigger on EXT1 channel; and trigger on positive edge - sp.spcm_dwSetParam_i32(self.card, sp.SPC_TRIG_EXT1_MODE, sp.SPC_TM_POS) - sp.spcm_dwSetParam_i32(self.card, sp.SPC_TRIG_ORMASK, sp.SPC_TMASK_EXT1) - - # Setup gated FIFO mode - sp.spcm_dwSetParam_i32(self.card, sp.SPC_CARDMODE, sp.SPC_REC_FIFO_GATE) - - # Get gate length alignment, number of samples must be integer multiple of this - gate_alignment = sp.int64(0) - sp.spcm_dwGetParam_i64(self.card, sp.SPC_GATE_LEN_ALIGNMENT, byref(gate_alignment)) - self.gate_alignment = gate_alignment.value - self.log.debug("Alignment samples: %d samples" % (self.gate_alignment)) - - # Set timeout used for DMA wait to 10 ms - sp.spcm_dwSetParam_i32(self.card, sp.SPC_TIMEOUT, 10) - - self.log.debug("Device setup completed") - - def start_operation( - self, - submit_fn: Callable[[int, RxData], None] | None = None, - index_offset: int = 0, - ) -> None: - """Start card operation. - - Parameters - ---------- - submit_fn - Optional callback invoked as ``submit_fn(global_index, rx_data)`` - for each completed gate. When provided the gate's slot in - ``rx_data`` is set to ``None`` immediately after the call so the - caller owns the object. When ``None`` the populated items remain - in ``rx_data`` for the caller to collect after ``stop_operation()``. - index_offset - Added to the per-gate index before calling *submit_fn*, allowing - the caller to assign globally unique indices across multiple averages. - """ - self._submit_fn = submit_fn - self._index_offset = index_offset - self.is_running.clear() - self.is_receiving.clear() - self.worker = threading.Thread(target=self._gated_timestamps_stream) - self.worker.start() - - def stop_operation(self): - """Stop card thread.""" - if self.worker is not None: - self.is_running.set() - self.worker.join() - self._submit_fn = None - self._index_offset = 0 - - # Stop card operation with the following steps: - # 1. Stop card acquisition - # 2. Stop data DMA transfer - # 3. Stop timestamp DMA transfer - self.handle_error( - sp.spcm_dwSetParam_i32( - self.card, - sp.SPC_M2CMD, - sp.M2CMD_CARD_STOP | sp.M2CMD_DATA_STOPDMA | sp.M2CMD_EXTRA_STOPDMA, - ) - ) - else: - # No thread is running - self.log.error("No active process found") - - def _gated_timestamps_stream(self): - # Rx buffer size must be a multiple of notify size. Min. notify size is 4096 bytes/4 kBytes. - rx_notify = sp.int32(sp.KILO_B(4)) - - # Buffer size set to maximum. - rx_size = 1024**3 - rx_buffer_size = sp.uint64(rx_size) - - # Create DMA buffer for receive data and tell the card to use it - rx_buffer = create_dma_buffer(rx_buffer_size.value) - sp.spcm_dwDefTransfer_i64( - self.card, - sp.SPCM_BUF_DATA, - sp.SPCM_DIR_CARDTOPC, - rx_notify, - rx_buffer, - sp.uint64(0), - rx_buffer_size, - ) - - # Define the timestamps notify size. Min. notify size is 4096 bytes. - ts_notify = sp.int32(sp.KILO_B(4)) - # Define timestamp buffer size, must be multiple of timestamps notify size - ts_buffer_size = sp.uint64(2 * 4096) - - # Create DMA buffer for timestamp data and tell the card to use it - ts_buffer = create_dma_buffer(ts_buffer_size.value) - sp.spcm_dwDefTransfer_i64( - self.card, - sp.SPCM_BUF_TIMESTAMP, - sp.SPCM_DIR_CARDTOPC, - ts_notify, - ts_buffer, - sp.uint64(0), - ts_buffer_size, - ) - - pll_data = cast(ts_buffer, sp.ptr64) # cast to pointer to 64bit integer - adc_data = cast(rx_buffer, sp.ptr16) # cast to pointer to 16bit integer - - # Setup polling mode for timestamp data - self.handle_error(sp.spcm_dwSetParam_i32(self.card, sp.SPC_M2CMD, sp.M2CMD_EXTRA_POLL)) - - # Start card acquisition and DMA usage - self.handle_error( - sp.spcm_dwSetParam_i32( - self.card, - sp.SPC_M2CMD, - sp.M2CMD_CARD_START | sp.M2CMD_CARD_ENABLETRIGGER | sp.M2CMD_DATA_STARTDMA, - ) - ) - - # Define helpers/buffer to read card parameter - available_timestamp_bytes = sp.int32(0) - available_timestamp_position = sp.int32(0) - available_data_bytes = sp.int32(0) - available_data_position = sp.int32(0) - - # Track bytes from incomplete gate reads for next iteration - remaining_bytes = 0 - # Track the amount of gate events recorded - self._total_gates = 0 - - # Check that the list of RxData objects has been passed - if self.rx_data is None: - self.log.critical("No RxData objects found for storing ADC data") - raise RuntimeError("No RxData objects found for storing ADC data") - - # Signal that acquisition has started - self.log.debug("Starting receive") - self.is_receiving.set() - - while not self.is_running.is_set(): - # Read the available timestamp buffer size - sp.spcm_dwGetParam_i64(self.card, sp.SPC_TS_AVAIL_USER_LEN, byref(available_timestamp_bytes)) - - # Process, if buffer size is greater or equal 32 (corresponds to 2 timestamps) - if available_timestamp_bytes.value >= 32: - # Read timestamp position - sp.spcm_dwGetParam_i32( - self.card, - sp.SPC_TS_AVAIL_USER_POS, - byref(available_timestamp_position), - ) - - # Read exactly two timestamps - timestamp_0 = pll_data[int(available_timestamp_position.value / 8)] - timestamp_1 = pll_data[int(available_timestamp_position.value / 8) + 2] - - # Calculate gate duration and the number of adc gate sample points (per channel) - num_gate_samples = timestamp_1 - timestamp_0 - gate_duration = num_gate_samples / (self.sample_rate * 1e6) - - self.log.info( - "Gate: (%s s, %s s); ADC duration: %s ms ; Samples/gate/channel: %s", - timestamp_0 / (self.sample_rate * 1e6), - timestamp_1 / (self.sample_rate * 1e6), - float(gate_duration) * 1e3, # Can be trimmed. - num_gate_samples, - ) - - # Tell buffer 32 bytes were read from timestamp buffer - try: - self.handle_error(sp.spcm_dwSetParam_i32(self.card, sp.SPC_TS_AVAIL_CARD_LEN, 32)) - except RuntimeError: # Reraise error for traceability - raise RuntimeError - - # Calculate size of relevant data (pre_trigger needed to get position of start of gate) - # This is the minimum amount of data must be available to get full gate data - total_bytes_gate = (num_gate_samples + self.pre_trigger) * 2 * self.num_channels.value - # Get the total data duration, including post trigger, to accurately track buffer position - samples_sequence = num_gate_samples + self.pre_trigger + self.post_trigger - # Ensure data alignment - alignment_samples = samples_sequence % self.gate_alignment - samples_sequence += alignment_samples - bytes_sequence = samples_sequence * 2 * self.num_channels.value - - # Check if total gate data does not exceed buffer size - if bytes_sequence > rx_size: - error_msg = ( - f"ADC gate data ({bytes_sequence} bytes) exceeds " - f"available buffer ({rx_size} bytes). " - f"Reduce adc length, sample rate or channel count" - ) - self.log.critical(error_msg) - raise ValueError(error_msg) - - # Wait for ADC data to arrive in DMA buffer - try: - self.handle_error(sp.spcm_dwSetParam_i32(self.card, sp.SPC_M2CMD, sp.M2CMD_DATA_WAITDMA)) - except RuntimeError as e: # Reraise error for traceability - self.log.error(f"DMA wait failed with error: {e}") - break - - # Read available data length and position - sp.spcm_dwGetParam_i32(self.card, sp.SPC_DATA_AVAIL_USER_POS, byref(available_data_position)) - sp.spcm_dwGetParam_i32(self.card, sp.SPC_DATA_AVAIL_USER_LEN, byref(available_data_bytes)) - - # # Debug log statements - self.log.debug( - "ADC event size: %d bytes, Available data length: %s bytes" - % (total_bytes_gate, available_data_bytes.value) - ) - - # If insufficient data is in buffer wait for more to arrive. - if available_data_bytes.value + remaining_bytes < total_bytes_gate: - # Wait for sufficient data to come in - wait_start = time.time() - while ( - available_data_bytes.value + remaining_bytes < total_bytes_gate - ) and not self.is_running.is_set(): - try: - self.handle_error(sp.spcm_dwSetParam_i32(self.card, sp.SPC_M2CMD, sp.M2CMD_DATA_WAITDMA)) - except RuntimeError as e: # Reraise error for traceability - self.log.error(f"DMA wait failed with error: {e}") - break - sp.spcm_dwGetParam_i32(self.card, sp.SPC_DATA_AVAIL_USER_LEN, byref(available_data_bytes)) - self.log.debug(f"Waited {(time.time() - wait_start) * 1e3:.3f} ms for extra data to enter buffer") - - if remaining_bytes + available_data_bytes.value > rx_size: - error_msg = ( - f"Memory overflow. Sum of remaining bytes ({remaining_bytes} bytes) " - f"and newly available bytes ({available_data_bytes.value} bytes) " - f"exceeds receive buffer size ({rx_size} bytes)" - ) - self.log.critical(error_msg) - raise MemoryError(error_msg) - - # Check if sufficient data is available (while loop doesn't guarantee it since it can be interrupted) - if available_data_bytes.value + remaining_bytes >= total_bytes_gate: - # Adjust memory position to account for bytes remaining after previous acquisition - byte_position = available_data_position.value - remaining_bytes - - # Handle buffer wraparound - if byte_position + total_bytes_gate >= rx_size: - # Calculate number of bytes to end of buffer - bytes_to_end = rx_size - byte_position - # calculates number of samples to end of buffer (2 bytes per sample) - samples_to_end = bytes_to_end // 2 - # Get the remaining number of samples after overflow - samples_leftover = total_bytes_gate // 2 - samples_to_end - - # Get the first part of the data - # Handle edge case when memory position is exactly at end - if samples_to_end == 0: - slice_1 = np.array([], dtype=np.int16) - else: - ptr_to_slice_1 = cast(addressof(adc_data.contents) + byte_position, POINTER(c_short)) - slice_1 = np.ctypeslib.as_array(ptr_to_slice_1, (samples_to_end,)) - - # Get the second part of the numpy slice - ptr_to_slice_2 = cast(addressof(adc_data.contents), POINTER(c_short)) - slice_2 = np.ctypeslib.as_array(ptr_to_slice_2, (samples_leftover,)) - - # Combine the slices - gate_data = np.concatenate((slice_1, slice_2)) - - else: - # If there is no memory position overflow, just get the data. - ptr_to_slice = cast(addressof(adc_data.contents) + byte_position, POINTER(c_short)) - gate_data = np.ctypeslib.as_array(ptr_to_slice, ((total_bytes_gate // 2),)) - - # Cut the pretrigger (we don't need it) and reshape the data to (num_coils, num_samples) - pre_trigger_cut = (self.pre_trigger) * self.num_channels.value - gate_data = gate_data[pre_trigger_cut:].reshape( - (self.num_channels.value, num_gate_samples), - order="F", - ) - - gate = self.rx_data[self._total_gates] - if gate is None: - msg = f"RxData slot {self._total_gates} was already consumed by the processor." - raise RuntimeError(msg) - # Store digital reference signal first (bit 16 of channel 0) - reference_len = min(num_gate_samples, NUM_REFERENCE_SAMPLES) - gate.phase_reference = (gate_data[0, :reference_len].astype(np.uint16) >> 15).copy() - # Modify gate_data by removing digital signal before writing it to RxData instance - gate_data[0] = (gate_data[0].view(np.uint16) << 1).view(np.int16) - gate.write_raw_data(gate_data) - gate.scaling_factor = self.rx_scaling[: self.num_channels.value] - gate.time_stamp = timestamp_0 / (self.sample_rate * 1e6) - - if self._submit_fn is not None: - self._submit_fn(self._index_offset + self._total_gates, gate) - self.rx_data[self._total_gates] = None - - # The accumulation of the leftover bytes is positive, - # if if the post-trigger event was not fully captured (accumulated sum increases), - # or negative if more then the expected data could be read due to lefter bytes - # from a previous acquisition (accumulated sum decreases). - remaining_bytes += available_data_bytes.value - bytes_sequence - - self._total_gates += 1 - - # Tell the card that data has been read and the buffer can be reused. - # Using the size of available data bytes prevents invalid values. - try: - self.handle_error( - sp.spcm_dwSetParam_i32(self.card, sp.SPC_DATA_AVAIL_CARD_LEN, available_data_bytes) - ) - except RuntimeError: # Reraise error for traceability - raise RuntimeError - - else: - self.log.error( - "Needed at least %d bytes but only %d bytes available" - % (total_bytes_gate, available_data_bytes.value) - ) - - self.log.debug("Card operation stopped") +"""Implementation of receive card.""" + +import logging +import threading +import time +from collections.abc import Callable +from ctypes import POINTER, addressof, byref, c_short, cast +from dataclasses import dataclass +from itertools import compress + +import numpy as np + +import console.spcm_control.spcm.pyspcm as sp +from console.interfaces.rx_data import RxData +from console.pulseq_interpreter.sequence_provider import NUM_REFERENCE_SAMPLES +from console.spcm_control.abstract_device import SpectrumDevice +from console.spcm_control.spcm.tools import create_dma_buffer, type_to_name + +# Define registers lists +CH_SELECT = [ + sp.CHANNEL0, + sp.CHANNEL1, + sp.CHANNEL2, + sp.CHANNEL3, + sp.CHANNEL4, + sp.CHANNEL5, + sp.CHANNEL6, + sp.CHANNEL7, +] +AMP_SELECT = [ + sp.SPC_AMP0, + sp.SPC_AMP1, + sp.SPC_AMP2, + sp.SPC_AMP3, + sp.SPC_AMP4, + sp.SPC_AMP5, + sp.SPC_AMP6, + sp.SPC_AMP7, +] +IMP_SELECT = [ + sp.SPC_50OHM0, + sp.SPC_50OHM1, + sp.SPC_50OHM2, + sp.SPC_50OHM3, + sp.SPC_50OHM4, + sp.SPC_50OHM5, + sp.SPC_50OHM6, + sp.SPC_50OHM7, +] + + +@dataclass +class RxCard(SpectrumDevice): + """Implementation of RX device.""" + + __name__: str = "RxCard" + + def __init__( + self, + path: str, + sample_rate: int, + channel_enable: tuple[bool, ...], + max_amplitude: tuple[int, ...], + impedance_50_ohms: tuple[bool, ...], + ) -> None: + """Execute after init function to do further class setup.""" + self.log = logging.getLogger(self.__name__) + super().__init__(path, log=self.log) + + self.sample_rate = sample_rate + self.channel_enable = [int(val) for val in channel_enable] + self.max_amplitude = max_amplitude + self.impedance_50_ohms = [int(val) for val in impedance_50_ohms] + self.rx_data: list[RxData | None] | None = None + + self.num_channels = sp.int32(0) + self.card_type = sp.int32(0) + + self.worker: threading.Thread | None = None + self.is_running = threading.Event() + self.is_receiving = threading.Event() + self._total_gates: int = 0 + + # Pre trigger is set to minimum, post trigger depends on active channel count and is defined later. + self.pre_trigger: int = 8 + self.post_trigger: None | int = None + + self.rx_scaling = [amp / (2**15) for amp in self.max_amplitude] + + self._submit_fn: Callable[[int, RxData], None] | None = None + self._index_offset: int = 0 + + @property + def total_gates(self) -> int: + """Helper function to return the number of gates that have been collected by the Rx Card.""" + return self._total_gates + + def setup_card(self): + """Set up spectrum card in transmit (Rx) mode. + + At the very beginning, a card reset is performed. The clock mode is set according to the sample rate, + defined by the class attribute. + Two receive channels are enables and configured by max. amplitude according to class variables and impedance. + + Raises + ------ + Warning + The actual set sample rate deviates from the corresponding class attribute to be set, + class attribute is overwritten. + """ + # Get the card type and reset card + sp.spcm_dwGetParam_i32(self.card, sp.SPC_PCITYP, byref(self.card_type)) + sp.spcm_dwSetParam_i64(self.card, sp.SPC_M2CMD, sp.M2CMD_CARD_RESET) # Needed? + + try: + if "M2p.59" not in (device_type := type_to_name(self.card_type.value)): + raise ConnectionError("Device with path %s is of type %s, no receive card" % (self.path, device_type)) + except ConnectionError as err: + self.log.exception(err, exc_info=True) + raise err + + # Setup the internal clockmode, clock output enable (use RX clock output to enable anti-alias filter) + sp.spcm_dwSetParam_i32(self.card, sp.SPC_CLOCKMODE, sp.SPC_CM_INTPLL) + sp.spcm_dwSetParam_i32(self.card, sp.SPC_CLOCKOUT, 1) + + # Use external clock: Terminate to 50 Ohms, set threshold to 1.5V, suitable for 3.3V clock + # sp.spcm_dwSetParam_i32(self.card, sp.SPC_CLOCKMODE, sp.SPC_CM_EXTERNAL) + # sp.spcm_dwSetParam_i32(self.card, sp.SPC_CLOCK50OHM, 1) + # sp.spcm_dwSetParam_i32(self.card, sp.SPC_CLOCK_THRESHOLD, 1500) + + # Set card sampling rate in MHz and read the actual sampling rate + sp.spcm_dwSetParam_i64(self.card, sp.SPC_SAMPLERATE, sp.MEGA(self.sample_rate)) + sample_rate = sp.int64(0) + sp.spcm_dwGetParam_i64(self.card, sp.SPC_SAMPLERATE, byref(sample_rate)) + self.log.info("Device sampling rate: %s MHz", sample_rate.value * 1e-6) + + if sample_rate.value != sp.MEGA(self.sample_rate): + self.log.warning( + "Actual device sample rate %s MHz does not match set sample rate of %s MHz; Updating class attribute", + sample_rate.value * 1e-6, + self.sample_rate, + ) + self.sample_rate = int(sample_rate.value * 1e-6) + + # Check channel enable, max. amplitude per channel and impedance values + try: + # Check that the length of the channel enable list is 8 + # this has to be true for cards with fewer channels too + if (num_enable := len(self.channel_enable)) != 8: + raise ValueError("Channel enable list is incomplete: %s/8" % num_enable) + # Impedance and amplitude configuration lists must also be of length 8 + if (num_imp := len(self.impedance_50_ohms)) != 8: + raise ValueError("Channel impedance list is incomplete: %s/8" % num_imp) + if (num_amp := len(self.max_amplitude)) != 8: + raise ValueError("channel max. amplitude list is incomplete: %s/8" % num_amp) + # Number of enabled channels must be either 1, 2, 4 or 8 + if not np.log2(sum(self.channel_enable)).is_integer(): + raise ValueError("Invalid number of enabled channels, must be power of 2.") + except ValueError as err: + self.log.exception(err, exc_info=True) + raise err + + # Enable receive channels, compress list of channel select registers to obtain list of channels to be enabled + # Sum of the compressed list equals logical or operator + # e.g. sp.CHANNEL0 | sp.CHANNEL1 | sp.CHANNEL5 = sum([sp.CHANNEL0, sp.CHANNEL1, sp.CHANNEL5]) = 35 + channel_selection = sum(list(compress(CH_SELECT, map(bool, self.channel_enable)))) + sp.spcm_dwSetParam_i32(self.card, sp.SPC_CHENABLE, channel_selection) + + # Set impedance and amplitude limits for each channel according to device configuration + for k, enable in enumerate(map(bool, self.channel_enable)): + if enable: + self.log.info( + "Channel %s enabled; 50 ohms impedance: %s; Max. amplitude: %s mV", + k, + self.impedance_50_ohms[k], + self.max_amplitude[k], + ) + sp.spcm_dwSetParam_i32(self.card, IMP_SELECT[k], self.impedance_50_ohms[k]) + sp.spcm_dwSetParam_i32(self.card, AMP_SELECT[k], self.max_amplitude[k]) + + # Get the number of actual active channels and compare against provided channel enable list + sp.spcm_dwGetParam_i32(self.card, sp.SPC_CHCOUNT, byref(self.num_channels)) + try: + self.log.info( + "Number of enabled receive channels (read from card): %s", + self.num_channels.value, + ) + if not self.num_channels.value == sum(self.channel_enable): + raise ValueError("Actual number of enabled channels does not match the provided channel enable list") + except ValueError as err: + self.log.exception(err, exc_info=True) + raise err + + # Digital filter setting for receiver, 0 = disable digital bandwidth filter + sp.spcm_dwSetParam_i32(self.card, sp.SPC_DIGITALBWFILTER, 0) + + # Configure X2 as digital input for phase reference signal and sample it in sync with analog channel 0 + sp.spcm_dwSetParam_i32(self.card, sp.SPCM_X2_MODE, sp.SPCM_XMODE_DIGIN) + sp.spcm_dwSetParam_i32(self.card, sp.SPC_DIGMODE0, (sp.DIGMODEMASK_BIT15 & sp.SPCM_DIGMODE_X2)) + + # Calculate trigger size depending on the number of active channels + # Since data can only be gathered in notify size chunks, post_trigger // channel_count should be at least one + # notify size to ensure that we can always access the full gate data. + self.post_trigger = 4096 // self.num_channels.value + + # Set the memory size, pre and post trigger and loop parameters, SPC_LOOPS = 0 => runs infinitely long + sp.spcm_dwSetParam_i32(self.card, sp.SPC_POSTTRIGGER, self.post_trigger) + sp.spcm_dwSetParam_i32(self.card, sp.SPC_PRETRIGGER, self.pre_trigger) + sp.spcm_dwSetParam_i32(self.card, sp.SPC_LOOPS, 0) + + # Setup timestamp mode to read number of samples per gate if available + sp.spcm_dwSetParam_i32( + self.card, + sp.SPC_TIMESTAMP_CMD, + sp.SPC_TSMODE_STARTRESET | sp.SPC_TSCNT_INTERNAL, + ) + # Configure trigger on EXT1 channel; and trigger on positive edge + sp.spcm_dwSetParam_i32(self.card, sp.SPC_TRIG_EXT1_MODE, sp.SPC_TM_POS) + sp.spcm_dwSetParam_i32(self.card, sp.SPC_TRIG_ORMASK, sp.SPC_TMASK_EXT1) + + # Setup gated FIFO mode + sp.spcm_dwSetParam_i32(self.card, sp.SPC_CARDMODE, sp.SPC_REC_FIFO_GATE) + + # Get gate length alignment, number of samples must be integer multiple of this + gate_alignment = sp.int64(0) + sp.spcm_dwGetParam_i64(self.card, sp.SPC_GATE_LEN_ALIGNMENT, byref(gate_alignment)) + self.gate_alignment = gate_alignment.value + self.log.debug("Alignment samples: %d samples" % (self.gate_alignment)) + + # Set timeout used for DMA wait to 10 ms + sp.spcm_dwSetParam_i32(self.card, sp.SPC_TIMEOUT, 10) + + self.log.debug("Device setup completed") + + def start_operation( + self, + submit_fn: Callable[[int, RxData], None] | None = None, + index_offset: int = 0, + ) -> None: + """Start card operation. + + Parameters + ---------- + submit_fn + Optional callback invoked as ``submit_fn(global_index, rx_data)`` + for each completed gate. When provided the gate's slot in + ``rx_data`` is set to ``None`` immediately after the call so the + caller owns the object. When ``None`` the populated items remain + in ``rx_data`` for the caller to collect after ``stop_operation()``. + index_offset + Added to the per-gate index before calling *submit_fn*, allowing + the caller to assign globally unique indices across multiple averages. + """ + self._submit_fn = submit_fn + self._index_offset = index_offset + self.is_running.clear() + self.is_receiving.clear() + self.worker = threading.Thread(target=self._gated_timestamps_stream) + self.worker.start() + + def stop_operation(self): + """Stop card thread.""" + if self.worker is not None: + self.is_running.set() + self.worker.join() + self._submit_fn = None + self._index_offset = 0 + + # Stop card operation with the following steps: + # 1. Stop card acquisition + # 2. Stop data DMA transfer + # 3. Stop timestamp DMA transfer + self.handle_error( + sp.spcm_dwSetParam_i32( + self.card, + sp.SPC_M2CMD, + sp.M2CMD_CARD_STOP | sp.M2CMD_DATA_STOPDMA | sp.M2CMD_EXTRA_STOPDMA, + ) + ) + else: + # No thread is running + self.log.error("No active process found") + + def _gated_timestamps_stream(self): + # Rx buffer size must be a multiple of notify size. Min. notify size is 4096 bytes/4 kBytes. + rx_notify = sp.int32(sp.KILO_B(4)) + + # Buffer size set to maximum. + rx_size = 1024**3 + rx_buffer_size = sp.uint64(rx_size) + + # Create DMA buffer for receive data and tell the card to use it + rx_buffer = create_dma_buffer(rx_buffer_size.value) + sp.spcm_dwDefTransfer_i64( + self.card, + sp.SPCM_BUF_DATA, + sp.SPCM_DIR_CARDTOPC, + rx_notify, + rx_buffer, + sp.uint64(0), + rx_buffer_size, + ) + + # Define the timestamps notify size. Min. notify size is 4096 bytes. + ts_notify = sp.int32(sp.KILO_B(4)) + # Define timestamp buffer size, must be multiple of timestamps notify size + ts_buffer_size = sp.uint64(2 * 4096) + + # Create DMA buffer for timestamp data and tell the card to use it + ts_buffer = create_dma_buffer(ts_buffer_size.value) + sp.spcm_dwDefTransfer_i64( + self.card, + sp.SPCM_BUF_TIMESTAMP, + sp.SPCM_DIR_CARDTOPC, + ts_notify, + ts_buffer, + sp.uint64(0), + ts_buffer_size, + ) + + pll_data = cast(ts_buffer, sp.ptr64) # cast to pointer to 64bit integer + adc_data = cast(rx_buffer, sp.ptr16) # cast to pointer to 16bit integer + + # Setup polling mode for timestamp data + self.handle_error(sp.spcm_dwSetParam_i32(self.card, sp.SPC_M2CMD, sp.M2CMD_EXTRA_POLL)) + + # Start card acquisition and DMA usage + self.handle_error( + sp.spcm_dwSetParam_i32( + self.card, + sp.SPC_M2CMD, + sp.M2CMD_CARD_START | sp.M2CMD_CARD_ENABLETRIGGER | sp.M2CMD_DATA_STARTDMA, + ) + ) + + # Define helpers/buffer to read card parameter + available_timestamp_bytes = sp.int32(0) + available_timestamp_position = sp.int32(0) + available_data_bytes = sp.int32(0) + available_data_position = sp.int32(0) + + # Track bytes from incomplete gate reads for next iteration + remaining_bytes = 0 + # Track the amount of gate events recorded + self._total_gates = 0 + + # Check that the list of RxData objects has been passed + if self.rx_data is None: + self.log.critical("No RxData objects found for storing ADC data") + raise RuntimeError("No RxData objects found for storing ADC data") + + # Signal that acquisition has started + self.log.debug("Starting receive") + self.is_receiving.set() + + while not self.is_running.is_set(): + # Read the available timestamp buffer size + sp.spcm_dwGetParam_i64(self.card, sp.SPC_TS_AVAIL_USER_LEN, byref(available_timestamp_bytes)) + + # Process, if buffer size is greater or equal 32 (corresponds to 2 timestamps) + if available_timestamp_bytes.value >= 32: + # Read timestamp position + sp.spcm_dwGetParam_i32( + self.card, + sp.SPC_TS_AVAIL_USER_POS, + byref(available_timestamp_position), + ) + + # Read exactly two timestamps + timestamp_0 = pll_data[int(available_timestamp_position.value / 8)] + timestamp_1 = pll_data[int(available_timestamp_position.value / 8) + 2] + + # Calculate gate duration and the number of adc gate sample points (per channel) + num_gate_samples = timestamp_1 - timestamp_0 + gate_duration = num_gate_samples / (self.sample_rate * 1e6) + + self.log.info( + "Gate: (%s s, %s s); ADC duration: %s ms ; Samples/gate/channel: %s", + timestamp_0 / (self.sample_rate * 1e6), + timestamp_1 / (self.sample_rate * 1e6), + float(gate_duration) * 1e3, # Can be trimmed. + num_gate_samples, + ) + + # Tell buffer 32 bytes were read from timestamp buffer + try: + self.handle_error(sp.spcm_dwSetParam_i32(self.card, sp.SPC_TS_AVAIL_CARD_LEN, 32)) + except RuntimeError: # Reraise error for traceability + raise RuntimeError + + # Calculate size of relevant data (pre_trigger needed to get position of start of gate) + # This is the minimum amount of data must be available to get full gate data + total_bytes_gate = (num_gate_samples + self.pre_trigger) * 2 * self.num_channels.value + # Get the total data duration, including post trigger, to accurately track buffer position + samples_sequence = num_gate_samples + self.pre_trigger + self.post_trigger + # Ensure data alignment + alignment_samples = samples_sequence % self.gate_alignment + samples_sequence += alignment_samples + bytes_sequence = samples_sequence * 2 * self.num_channels.value + + # Check if total gate data does not exceed buffer size + if bytes_sequence > rx_size: + error_msg = ( + f"ADC gate data ({bytes_sequence} bytes) exceeds " + f"available buffer ({rx_size} bytes). " + f"Reduce adc length, sample rate or channel count" + ) + self.log.critical(error_msg) + raise ValueError(error_msg) + + # Wait for ADC data to arrive in DMA buffer + try: + self.handle_error(sp.spcm_dwSetParam_i32(self.card, sp.SPC_M2CMD, sp.M2CMD_DATA_WAITDMA)) + except RuntimeError as e: # Reraise error for traceability + self.log.error(f"DMA wait failed with error: {e}") + break + + # Read available data length and position + sp.spcm_dwGetParam_i32(self.card, sp.SPC_DATA_AVAIL_USER_POS, byref(available_data_position)) + sp.spcm_dwGetParam_i32(self.card, sp.SPC_DATA_AVAIL_USER_LEN, byref(available_data_bytes)) + + # # Debug log statements + self.log.debug( + "ADC event size: %d bytes, Available data length: %s bytes" + % (total_bytes_gate, available_data_bytes.value) + ) + + # If insufficient data is in buffer wait for more to arrive. + if available_data_bytes.value + remaining_bytes < total_bytes_gate: + # Wait for sufficient data to come in + wait_start = time.time() + while ( + available_data_bytes.value + remaining_bytes < total_bytes_gate + ) and not self.is_running.is_set(): + try: + self.handle_error(sp.spcm_dwSetParam_i32(self.card, sp.SPC_M2CMD, sp.M2CMD_DATA_WAITDMA)) + except RuntimeError as e: # Reraise error for traceability + self.log.error(f"DMA wait failed with error: {e}") + break + sp.spcm_dwGetParam_i32(self.card, sp.SPC_DATA_AVAIL_USER_LEN, byref(available_data_bytes)) + self.log.debug(f"Waited {(time.time() - wait_start) * 1e3:.3f} ms for extra data to enter buffer") + + if remaining_bytes + available_data_bytes.value > rx_size: + error_msg = ( + f"Memory overflow. Sum of remaining bytes ({remaining_bytes} bytes) " + f"and newly available bytes ({available_data_bytes.value} bytes) " + f"exceeds receive buffer size ({rx_size} bytes)" + ) + self.log.critical(error_msg) + raise MemoryError(error_msg) + + # Check if sufficient data is available (while loop doesn't guarantee it since it can be interrupted) + if available_data_bytes.value + remaining_bytes >= total_bytes_gate: + # Adjust memory position to account for bytes remaining after previous acquisition + byte_position = available_data_position.value - remaining_bytes + + # Handle buffer wraparound + if byte_position + total_bytes_gate >= rx_size: + # Calculate number of bytes to end of buffer + bytes_to_end = rx_size - byte_position + # calculates number of samples to end of buffer (2 bytes per sample) + samples_to_end = bytes_to_end // 2 + # Get the remaining number of samples after overflow + samples_leftover = total_bytes_gate // 2 - samples_to_end + + # Get the first part of the data + # Handle edge case when memory position is exactly at end + if samples_to_end == 0: + slice_1 = np.array([], dtype=np.int16) + else: + ptr_to_slice_1 = cast(addressof(adc_data.contents) + byte_position, POINTER(c_short)) + slice_1 = np.ctypeslib.as_array(ptr_to_slice_1, (samples_to_end,)) + + # Get the second part of the numpy slice + ptr_to_slice_2 = cast(addressof(adc_data.contents), POINTER(c_short)) + slice_2 = np.ctypeslib.as_array(ptr_to_slice_2, (samples_leftover,)) + + # Combine the slices + gate_data = np.concatenate((slice_1, slice_2)) + + else: + # If there is no memory position overflow, just get the data. + ptr_to_slice = cast(addressof(adc_data.contents) + byte_position, POINTER(c_short)) + gate_data = np.ctypeslib.as_array(ptr_to_slice, ((total_bytes_gate // 2),)) + + # Cut the pretrigger (we don't need it) and reshape the data to (num_coils, num_samples) + pre_trigger_cut = (self.pre_trigger) * self.num_channels.value + gate_data = gate_data[pre_trigger_cut:].reshape( + (self.num_channels.value, num_gate_samples), + order="F", + ) + + gate = self.rx_data[self._total_gates] + if gate is None: + msg = f"RxData slot {self._total_gates} was already consumed by the processor." + raise RuntimeError(msg) + # Store digital reference signal first (bit 16 of channel 0) + reference_len = min(num_gate_samples, NUM_REFERENCE_SAMPLES) + gate.phase_reference = (gate_data[0, :reference_len].astype(np.uint16) >> 15).copy() + # Modify gate_data by removing digital signal before writing it to RxData instance + gate_data[0] = (gate_data[0].view(np.uint16) << 1).view(np.int16) + gate.write_raw_data(gate_data) + gate.scaling_factor = self.rx_scaling[: self.num_channels.value] + gate.time_stamp = timestamp_0 / (self.sample_rate * 1e6) + + if self._submit_fn is not None: + self._submit_fn(self._index_offset + self._total_gates, gate) + self.rx_data[self._total_gates] = None + + # The accumulation of the leftover bytes is positive, + # if if the post-trigger event was not fully captured (accumulated sum increases), + # or negative if more then the expected data could be read due to lefter bytes + # from a previous acquisition (accumulated sum decreases). + remaining_bytes += available_data_bytes.value - bytes_sequence + + self._total_gates += 1 + + # Tell the card that data has been read and the buffer can be reused. + # Using the size of available data bytes prevents invalid values. + try: + self.handle_error( + sp.spcm_dwSetParam_i32(self.card, sp.SPC_DATA_AVAIL_CARD_LEN, available_data_bytes) + ) + except RuntimeError: # Reraise error for traceability + raise RuntimeError + + else: + self.log.error( + "Needed at least %d bytes but only %d bytes available" + % (total_bytes_gate, available_data_bytes.value) + ) + + self.log.debug("Card operation stopped") diff --git a/src/console/spcm_control/rx_processor.py b/src/console/spcm_control/rx_processor.py index a275c854..310d2f24 100644 --- a/src/console/spcm_control/rx_processor.py +++ b/src/console/spcm_control/rx_processor.py @@ -1,112 +1,112 @@ -"""Processing worker for RxData objects using a persistent process pool.""" - -import logging -import multiprocessing -import signal -from concurrent.futures import Future, ProcessPoolExecutor, wait - -from console.interfaces.rx_data import RxData - -log = logging.getLogger("RxProc") - -# 'spawn' is the only safe start method on Windows and avoids fork-related -# deadlocks in multi-threaded processes on Linux. -_mp_ctx = multiprocessing.get_context("spawn") - - -def _worker_init() -> None: - """Ignore SIGINT in worker processes so Ctrl-C is handled by the main process only.""" - signal.signal(signal.SIGINT, signal.SIG_IGN) - - -def _noop() -> None: - """No-op submitted to pre-warm pool worker processes.""" - - -def _process_one(index: int, rx_data: RxData, store_unprocessed: bool) -> tuple[int, RxData]: - """Process a single RxData item inside a worker process.""" - rx_data.process_data(store_unprocessed=store_unprocessed) - rx_data.materialize(keep=store_unprocessed) - return index, rx_data - - -class RxProcessor: - """Processes RxData items using a persistent pool of worker processes. - - ``submit()`` is thread-safe and can be called directly from the rx_card - streaming thread. ``collect()`` is called after all averages complete - (by which point every streaming thread has been joined) and waits for - all outstanding futures. - - Usage - ----- - >>> processor = RxProcessor(num_workers=2) - >>> processor.start() - >>> # pass processor.submit as submit_fn to rx_card.start_operation() - >>> results = processor.collect(expected_count=N, timeout=T) - >>> # Repeat collect() for each subsequent run. - >>> processor.shutdown() - """ - - def __init__(self, num_workers: int = 1) -> None: - self._num_workers = num_workers - self._executor: ProcessPoolExecutor | None = None - self._futures: dict[int, Future] = {} - - def start(self) -> None: - """Create the process pool and pre-warm all worker processes.""" - self._executor = ProcessPoolExecutor( - max_workers=self._num_workers, mp_context=_mp_ctx, initializer=_worker_init - ) - # Submit one no-op per worker to force all processes to spawn now so - # the first real acquisition doesn't pay the spawn cost. - warm = [self._executor.submit(_noop) for _ in range(self._num_workers)] - for f in warm: - f.result() - log.debug("RxProcessor pool started with %d worker(s)", self._num_workers) - - def submit(self, index: int, rx_data: RxData, store_unprocessed: bool) -> None: - """Submit one RxData item for async processing. - - Thread-safe: may be called from any thread, including the rx_card - streaming thread. - """ - if self._executor is None: - raise RuntimeError("RxProcessor is not started. Call start() first.") - self._futures[index] = self._executor.submit(_process_one, index, rx_data, store_unprocessed) - - def collect(self, expected_count: int, timeout: float = 60.0) -> list[RxData]: - """Wait for all submitted futures and return results ordered by index. - - Parameters - ---------- - expected_count - Number of RxData items expected in this batch. - timeout - Maximum seconds to wait for all futures to complete. - - Returns - ------- - List of processed RxData objects ordered by global index. - """ - done, not_done = wait(list(self._futures.values()), timeout=timeout) - if not_done: - log.warning("%d future(s) did not complete within timeout", len(not_done)) - - results: dict[int, RxData] = {} - for future in done: - idx, rx_data = future.result() - results[idx] = rx_data - - if len(results) != expected_count: - log.warning("Expected %d items but collected %d", expected_count, len(results)) - - self._futures.clear() - return [results[i] for i in range(expected_count) if i in results] - - def shutdown(self) -> None: - """Shut down the process pool.""" - if self._executor is not None: - self._executor.shutdown(wait=False, cancel_futures=True) - self._executor = None - log.debug("RxProcessor pool shut down") +"""Processing worker for RxData objects using a persistent process pool.""" + +import logging +import multiprocessing +import signal +from concurrent.futures import Future, ProcessPoolExecutor, wait + +from console.interfaces.rx_data import RxData + +log = logging.getLogger("RxProc") + +# 'spawn' is the only safe start method on Windows and avoids fork-related +# deadlocks in multi-threaded processes on Linux. +_mp_ctx = multiprocessing.get_context("spawn") + + +def _worker_init() -> None: + """Ignore SIGINT in worker processes so Ctrl-C is handled by the main process only.""" + signal.signal(signal.SIGINT, signal.SIG_IGN) + + +def _noop() -> None: + """No-op submitted to pre-warm pool worker processes.""" + + +def _process_one(index: int, rx_data: RxData, store_unprocessed: bool) -> tuple[int, RxData]: + """Process a single RxData item inside a worker process.""" + rx_data.process_data(store_unprocessed=store_unprocessed) + rx_data.materialize(keep=store_unprocessed) + return index, rx_data + + +class RxProcessor: + """Processes RxData items using a persistent pool of worker processes. + + ``submit()`` is thread-safe and can be called directly from the rx_card + streaming thread. ``collect()`` is called after all averages complete + (by which point every streaming thread has been joined) and waits for + all outstanding futures. + + Usage + ----- + >>> processor = RxProcessor(num_workers=2) + >>> processor.start() + >>> # pass processor.submit as submit_fn to rx_card.start_operation() + >>> results = processor.collect(expected_count=N, timeout=T) + >>> # Repeat collect() for each subsequent run. + >>> processor.shutdown() + """ + + def __init__(self, num_workers: int = 1) -> None: + self._num_workers = num_workers + self._executor: ProcessPoolExecutor | None = None + self._futures: dict[int, Future] = {} + + def start(self) -> None: + """Create the process pool and pre-warm all worker processes.""" + self._executor = ProcessPoolExecutor( + max_workers=self._num_workers, mp_context=_mp_ctx, initializer=_worker_init + ) + # Submit one no-op per worker to force all processes to spawn now so + # the first real acquisition doesn't pay the spawn cost. + warm = [self._executor.submit(_noop) for _ in range(self._num_workers)] + for f in warm: + f.result() + log.debug("RxProcessor pool started with %d worker(s)", self._num_workers) + + def submit(self, index: int, rx_data: RxData, store_unprocessed: bool) -> None: + """Submit one RxData item for async processing. + + Thread-safe: may be called from any thread, including the rx_card + streaming thread. + """ + if self._executor is None: + raise RuntimeError("RxProcessor is not started. Call start() first.") + self._futures[index] = self._executor.submit(_process_one, index, rx_data, store_unprocessed) + + def collect(self, expected_count: int, timeout: float = 60.0) -> list[RxData]: + """Wait for all submitted futures and return results ordered by index. + + Parameters + ---------- + expected_count + Number of RxData items expected in this batch. + timeout + Maximum seconds to wait for all futures to complete. + + Returns + ------- + List of processed RxData objects ordered by global index. + """ + done, not_done = wait(list(self._futures.values()), timeout=timeout) + if not_done: + log.warning("%d future(s) did not complete within timeout", len(not_done)) + + results: dict[int, RxData] = {} + for future in done: + idx, rx_data = future.result() + results[idx] = rx_data + + if len(results) != expected_count: + log.warning("Expected %d items but collected %d", expected_count, len(results)) + + self._futures.clear() + return [results[i] for i in range(expected_count) if i in results] + + def shutdown(self) -> None: + """Shut down the process pool.""" + if self._executor is not None: + self._executor.shutdown(wait=False, cancel_futures=True) + self._executor = None + log.debug("RxProcessor pool shut down") diff --git a/src/console/spcm_control/spcm/tools.py b/src/console/spcm_control/spcm/tools.py index c744b17a..508efda3 100644 --- a/src/console/spcm_control/spcm/tools.py +++ b/src/console/spcm_control/spcm/tools.py @@ -1,127 +1,127 @@ -"""Tools for spectrum card.""" - -from ctypes import * -from typing import Any - -# load registers for easier access -import console.spcm_control.spcm.registers as regs -from console.spcm_control.spcm.errors import ERR_OK, error_reg -from console.spcm_control.spcm.status import status_reg, status_reg_desc - - -def translate_status(status: int, include_desc: bool = False) -> tuple[dict[int, list[Any]], list[str]]: - """Translate integer value to readable status message. - - Parameters - ---------- - status - Status code from spectrum card device - - Returns - ------- - Description from user manual, default is unknown - """ - # Convert status code to 12-digit bit sequence in reversed order - # >> lowest bit comes first => correspondence to order in manual - bit_reg = list(reversed("{:012b}".format(status))) - - # Status codes are defined for card (0x1 ... 0x8) and data (0x100 ... 0x800) - # >> First 4 bits correspond to (0x1 ... 0x8) - # >> Last 4 bits correspond to (0x100 ... 0x800) - status_flags_card = [bool(int(b)) for b in bit_reg[:4]] - status_flags_data = [bool(int(b)) for b in bit_reg[-4:]] - status_flags = status_flags_card + status_flags_data - - # Construct status dictionary, include description depending on function argument - status_dict: dict[int, list] = {} - for k, (val, stat) in enumerate(status_reg.items()): - status_dict[val] = [status_flags[k], stat, status_reg_desc[val]] if include_desc else [status_flags[k], stat] - - return status_dict, bit_reg - - -def translate_error(error: int) -> str | None: - """Translate error code to description string from manual. - - Parameters - ---------- - error - Error code to be translated - - Returns - ------- - Error description string from user manual - """ - if error in error_reg.keys(): - if error_reg[error] is not ERR_OK: - return f"ERROR: {error_reg[error]}" - return "Unknown error" - - -def type_to_name(card_type: int) -> str: - """Name translation for card type. - - Parameters - ---------- - lCardType - Card code - - Returns - ------- - Card name as string - """ - version = card_type & regs.TYP_VERSIONMASK - code = card_type & regs.TYP_SERIESMASK - match code: - case regs.TYP_M2ISERIES: - return "M2i.%04x" % version - case regs.TYP_M2IEXPSERIES: - return "M2i.%04x-Exp" % version - case regs.TYP_M3ISERIES: - return "M3i.%04x" % version - case regs.TYP_M3IEXPSERIES: - return "M3i.%04x-Exp" % version - case regs.TYP_M4IEXPSERIES: - return "M4i.%04x-x8" % version - case regs.TYP_M4XEXPSERIES: - return "M4x.%04x-x4" % version - case regs.TYP_M2PEXPSERIES: - return "M2p.%04x-x4" % version - case regs.TYP_M5IEXPSERIES: - return "M5i.%04x-x16" % version - case _: - return "unknown type" - - -def create_dma_buffer(buffer_size: int): - """Allocate memory for page-aligned DMA buffer. - - Parameters - ---------- - buffer_size - Size of the buffer - - Returns - ------- - Buffer - """ - dwAlignment = 4096 - dwMask = dwAlignment - 1 - - # allocate non-aligned, slightly larger buffer - qwRequiredNonAlignedBytes = buffer_size * sizeof(c_char) + dwMask - pvNonAlignedBuf = (c_char * qwRequiredNonAlignedBytes)() - - # get offset of next aligned address in non-aligned buffer - misalignment = addressof(pvNonAlignedBuf) & dwMask - if misalignment: - dwOffset = dwAlignment - misalignment - else: - dwOffset = 0 - - aligned_buffer = (c_char * buffer_size).from_buffer(pvNonAlignedBuf, dwOffset) - - # zero the aligned buffer explicitly - memset(addressof(aligned_buffer), 0, buffer_size) - - return aligned_buffer +"""Tools for spectrum card.""" + +from ctypes import * +from typing import Any + +# load registers for easier access +import console.spcm_control.spcm.registers as regs +from console.spcm_control.spcm.errors import ERR_OK, error_reg +from console.spcm_control.spcm.status import status_reg, status_reg_desc + + +def translate_status(status: int, include_desc: bool = False) -> tuple[dict[int, list[Any]], list[str]]: + """Translate integer value to readable status message. + + Parameters + ---------- + status + Status code from spectrum card device + + Returns + ------- + Description from user manual, default is unknown + """ + # Convert status code to 12-digit bit sequence in reversed order + # >> lowest bit comes first => correspondence to order in manual + bit_reg = list(reversed("{:012b}".format(status))) + + # Status codes are defined for card (0x1 ... 0x8) and data (0x100 ... 0x800) + # >> First 4 bits correspond to (0x1 ... 0x8) + # >> Last 4 bits correspond to (0x100 ... 0x800) + status_flags_card = [bool(int(b)) for b in bit_reg[:4]] + status_flags_data = [bool(int(b)) for b in bit_reg[-4:]] + status_flags = status_flags_card + status_flags_data + + # Construct status dictionary, include description depending on function argument + status_dict: dict[int, list] = {} + for k, (val, stat) in enumerate(status_reg.items()): + status_dict[val] = [status_flags[k], stat, status_reg_desc[val]] if include_desc else [status_flags[k], stat] + + return status_dict, bit_reg + + +def translate_error(error: int) -> str | None: + """Translate error code to description string from manual. + + Parameters + ---------- + error + Error code to be translated + + Returns + ------- + Error description string from user manual + """ + if error in error_reg.keys(): + if error_reg[error] is not ERR_OK: + return f"ERROR: {error_reg[error]}" + return "Unknown error" + + +def type_to_name(card_type: int) -> str: + """Name translation for card type. + + Parameters + ---------- + lCardType + Card code + + Returns + ------- + Card name as string + """ + version = card_type & regs.TYP_VERSIONMASK + code = card_type & regs.TYP_SERIESMASK + match code: + case regs.TYP_M2ISERIES: + return "M2i.%04x" % version + case regs.TYP_M2IEXPSERIES: + return "M2i.%04x-Exp" % version + case regs.TYP_M3ISERIES: + return "M3i.%04x" % version + case regs.TYP_M3IEXPSERIES: + return "M3i.%04x-Exp" % version + case regs.TYP_M4IEXPSERIES: + return "M4i.%04x-x8" % version + case regs.TYP_M4XEXPSERIES: + return "M4x.%04x-x4" % version + case regs.TYP_M2PEXPSERIES: + return "M2p.%04x-x4" % version + case regs.TYP_M5IEXPSERIES: + return "M5i.%04x-x16" % version + case _: + return "unknown type" + + +def create_dma_buffer(buffer_size: int): + """Allocate memory for page-aligned DMA buffer. + + Parameters + ---------- + buffer_size + Size of the buffer + + Returns + ------- + Buffer + """ + dwAlignment = 4096 + dwMask = dwAlignment - 1 + + # allocate non-aligned, slightly larger buffer + qwRequiredNonAlignedBytes = buffer_size * sizeof(c_char) + dwMask + pvNonAlignedBuf = (c_char * qwRequiredNonAlignedBytes)() + + # get offset of next aligned address in non-aligned buffer + misalignment = addressof(pvNonAlignedBuf) & dwMask + if misalignment: + dwOffset = dwAlignment - misalignment + else: + dwOffset = 0 + + aligned_buffer = (c_char * buffer_size).from_buffer(pvNonAlignedBuf, dwOffset) + + # zero the aligned buffer explicitly + memset(addressof(aligned_buffer), 0, buffer_size) + + return aligned_buffer diff --git a/src/console/utilities/json_encoder.py b/src/console/utilities/json_encoder.py index f256e2e6..26e5baa7 100644 --- a/src/console/utilities/json_encoder.py +++ b/src/console/utilities/json_encoder.py @@ -1,26 +1,26 @@ -"""Implementation of custom JSON encoder.""" -import dataclasses -import json -from pathlib import Path - - -class JSONEncoder(json.JSONEncoder): - """JSON Encoder class.""" - - def default(self, obj) -> object: - """Encode object default method. - - Parameters - ---------- - o - Object to encode - - Returns - ------- - JSON encoded object - """ - if bool(dataclasses.is_dataclass(obj)) and not isinstance(obj, type): - return dataclasses.asdict(obj) - if isinstance(obj, Path): - return str(obj) - return super().default(obj) +"""Implementation of custom JSON encoder.""" +import dataclasses +import json +from pathlib import Path + + +class JSONEncoder(json.JSONEncoder): + """JSON Encoder class.""" + + def default(self, obj) -> object: + """Encode object default method. + + Parameters + ---------- + o + Object to encode + + Returns + ------- + JSON encoded object + """ + if bool(dataclasses.is_dataclass(obj)) and not isinstance(obj, type): + return dataclasses.asdict(obj) + if isinstance(obj, Path): + return str(obj) + return super().default(obj) diff --git a/tests/acquisition/test_ddc.py b/tests/acquisition/test_ddc.py index a7f38654..de2536e6 100644 --- a/tests/acquisition/test_ddc.py +++ b/tests/acquisition/test_ddc.py @@ -1,56 +1,56 @@ -"""Test digital down converter (DDC) functions.""" -import numpy as np -import pytest -from scipy import signal - -from console.utilities.ddc import filter_cic_fir_comp, filter_moving_average - - -@pytest.mark.parametrize("num_coils", [1, 2, 4]) -@pytest.mark.parametrize("num_samples", [8000, 9511, 80640]) -@pytest.mark.parametrize("decimation", [100, 200, 400]) -@pytest.mark.parametrize("overlap", [2, 4]) -def test_moving_average_filter(num_coils, num_samples, decimation, overlap, random_complex_data): - """Test moving average filter for various parameter configurations with FIR as reference.""" - input_data = random_complex_data(shape=(num_coils, num_samples)) - processed_avg = filter_moving_average(input_data, decimation=decimation, overlap=overlap) - processed_fir = signal.decimate(input_data, q=decimation, ftype="fir") - - assert processed_avg.shape == processed_fir.shape - assert np.iscomplex(processed_avg).all() - assert np.iscomplex(processed_fir).all() - - -@pytest.mark.parametrize("num_coils", [1, 2, 4]) -@pytest.mark.parametrize("num_samples", [8000, 9511, 80640]) -@pytest.mark.parametrize("decimation", [100, 200, 400]) -@pytest.mark.parametrize("filter_stages", [2, 3, 5]) -def test_cic_fir_comp(num_coils, num_samples, decimation, filter_stages, random_complex_data): - """Test CIC FIR filter composition with FIR as reference.""" - input_data = random_complex_data(shape=(num_coils, num_samples)) - processed_cic = filter_cic_fir_comp(input_data, decimation=decimation, number_of_stages=filter_stages) - processed_fir = signal.decimate(input_data, q=decimation, ftype="fir") - - assert processed_cic.shape == processed_fir.shape - assert np.iscomplex(processed_cic).all() - - -@pytest.mark.parametrize("fixture_name", ["rx_data_fid", "rx_data_trapezoid"]) -def test_rx_data_processing(fixture_name, request): - """Test RxData processing.""" - rx_data = request.getfixturevalue(fixture_name) - rx_data.process_data(store_unprocessed=True) - - # Get the envelope of the raw signal using the Hilbert transform and discard samples - raw_envelope = signal.hilbert(rx_data.raw_data[0]) - processed = rx_data.processed_data[0] - - if rx_data.num_samples_discard > 0: - processed = processed[rx_data.num_samples_discard:-rx_data.num_samples_discard] - discard_raw = int(rx_data.num_samples_discard * rx_data.decimation_factor) - raw_envelope = raw_envelope[discard_raw:-discard_raw] - - # Compare the 99th percentile of the absolute signals - val_raw = np.percentile(np.abs(raw_envelope), 99) - val_proc = np.percentile(np.abs(processed), 99) - np.testing.assert_allclose(val_proc, val_raw, rtol=0.01) +"""Test digital down converter (DDC) functions.""" +import numpy as np +import pytest +from scipy import signal + +from console.utilities.ddc import filter_cic_fir_comp, filter_moving_average + + +@pytest.mark.parametrize("num_coils", [1, 2, 4]) +@pytest.mark.parametrize("num_samples", [8000, 9511, 80640]) +@pytest.mark.parametrize("decimation", [100, 200, 400]) +@pytest.mark.parametrize("overlap", [2, 4]) +def test_moving_average_filter(num_coils, num_samples, decimation, overlap, random_complex_data): + """Test moving average filter for various parameter configurations with FIR as reference.""" + input_data = random_complex_data(shape=(num_coils, num_samples)) + processed_avg = filter_moving_average(input_data, decimation=decimation, overlap=overlap) + processed_fir = signal.decimate(input_data, q=decimation, ftype="fir") + + assert processed_avg.shape == processed_fir.shape + assert np.iscomplex(processed_avg).all() + assert np.iscomplex(processed_fir).all() + + +@pytest.mark.parametrize("num_coils", [1, 2, 4]) +@pytest.mark.parametrize("num_samples", [8000, 9511, 80640]) +@pytest.mark.parametrize("decimation", [100, 200, 400]) +@pytest.mark.parametrize("filter_stages", [2, 3, 5]) +def test_cic_fir_comp(num_coils, num_samples, decimation, filter_stages, random_complex_data): + """Test CIC FIR filter composition with FIR as reference.""" + input_data = random_complex_data(shape=(num_coils, num_samples)) + processed_cic = filter_cic_fir_comp(input_data, decimation=decimation, number_of_stages=filter_stages) + processed_fir = signal.decimate(input_data, q=decimation, ftype="fir") + + assert processed_cic.shape == processed_fir.shape + assert np.iscomplex(processed_cic).all() + + +@pytest.mark.parametrize("fixture_name", ["rx_data_fid", "rx_data_trapezoid"]) +def test_rx_data_processing(fixture_name, request): + """Test RxData processing.""" + rx_data = request.getfixturevalue(fixture_name) + rx_data.process_data(store_unprocessed=True) + + # Get the envelope of the raw signal using the Hilbert transform and discard samples + raw_envelope = signal.hilbert(rx_data.raw_data[0]) + processed = rx_data.processed_data[0] + + if rx_data.num_samples_discard > 0: + processed = processed[rx_data.num_samples_discard:-rx_data.num_samples_discard] + discard_raw = int(rx_data.num_samples_discard * rx_data.decimation_factor) + raw_envelope = raw_envelope[discard_raw:-discard_raw] + + # Compare the 99th percentile of the absolute signals + val_raw = np.percentile(np.abs(raw_envelope), 99) + val_proc = np.percentile(np.abs(processed), 99) + np.testing.assert_allclose(val_proc, val_raw, rtol=0.01)