TextureHelper, 3Dcuts, Angular Sp, STFT_SSA, Wavelet,

.clangd

@@ -0,0 +1,2 @@
+CompileFlags:
+  Add: [-std=c++20]

.vscode/tasks.json

@@ -31,7 +31,7 @@
         {
             "type": "shell",
             "label": "CMake: build",
-            "command": "build\\generators\\conanbuild.bat && cmake --build --preset conan-relwithdebinfo",
+            "command": "build\\generators\\conanbuild.bat && cmake --build --preset conan-release",
             "group": {
                 "kind": "build",
                 "isDefault": true
@@ -56,4 +56,4 @@
             }
         }
     ]
-}
+}

Backend/cuda_includes/morlet_wavelet.cuh

@@ -0,0 +1,8 @@
+#pragma once
+#include <cuda_runtime.h>
+#include <cuComplex.h>
+#include <math.h>
+#include "complex_utils.cuh"
+
+void createMorletKernel(cuComplex* d_kernel, int N, float dt,
+                        float scale, float omega0, const cudaStream_t stream);

Backend/cuda_includes/wavelet_transform.cuh

@@ -0,0 +1,6 @@
+#pragma once
+#include "complex_utils.cuh"
+#include "morlet_wavelet.cuh"
+#include "frame_desc.hh"
+
+void wavelet_transform(cuComplex* output, cuComplex* input, const cufftHandle plan1d,  const camera::FrameDescriptor& fd, int tranformation_size, const cudaStream_t stream, float target_freq);

Backend/cuda_sources/angular_spectrum.cu

@@ -28,7 +28,7 @@ void angular_spectrum_lens(cuFloatComplex* output,
     kernel_spectral_lens<<<lblocks, lthreads, 0, stream>>>(output, Nx, Ny, z, lambda, x_step, y_step);
     cudaXStreamSynchronize(stream);
     cudaCheckError();
-}
+} 
 
 void angular_spectrum(cuComplex* input,
                       cuComplex* output,
@@ -48,6 +48,7 @@ void angular_spectrum(cuComplex* input,
 
     // Lens and Mask already shifted
     // thus we don't have to shift the 'input' buffer each time
+
     apply_mask(input, lens, output, frame_res, batch_size, stream);
     if (store_frame)
     {

Backend/cuda_sources/masks.cu

@@ -4,6 +4,9 @@
 
 using camera::FrameDescriptor;
 
+#define PI_F 3.14159265358979323846f
+#include <math_constants.h>  // for M_PI
+
 __global__ void kernel_quadratic_lens(
     cuComplex* output, const uint lens_side_size, const float lambda, const float dist, const float pixel_size)
 {
@@ -39,6 +42,46 @@ __global__ void kernel_spectral_lens(cuFloatComplex* output,
     int x = blockIdx.x * blockDim.x + threadIdx.x;
     int y = blockIdx.y * blockDim.y + threadIdx.y;
 
+    if (x >= Nx || y >= Ny)
+        return;
+
+    // Physical frequencies in cycles per meter
+    float u_step = 1.0f / (Nx * x_step);
+    float v_step = 1.0f / (Ny * y_step);
+
+    // Center the frequency grid
+    float u = (x - Nx / 2.0f) * u_step;
+    float v = (y - Ny / 2.0f) * v_step;
+
+    // Convert to angular spatial frequencies
+    float kx = 2.0f * M_PI * u;
+    float ky = 2.0f * M_PI * v;
+
+    // Wave number
+    float k = 2.0f * M_PI / lambda;
+
+    // Compute kz (propagation along z)
+    float kz2 = k * k - kx * kx - ky * ky;
+    float kz = (kz2 > 0.0f) ? sqrtf(kz2) : 0.0f;
+
+    // Compute phase
+    float phase = kz * z;
+
+    // Store as complex exponential
+    output[y * Nx + x] = make_cuFloatComplex(cosf(phase), sinf(phase));
+}
+
+/* __global__ void kernel_spectral_lens(cuFloatComplex* output,
+                                     const int Nx,
+                                     const int Ny,
+                                     const float z,
+                                     const float lambda,
+                                     const float x_step,
+                                     const float y_step)
+{
+    int x = blockIdx.x * blockDim.x + threadIdx.x;
+    int y = blockIdx.y * blockDim.y + threadIdx.y;
+
     if (x < Nx && y < Ny)
     {
         float u_step = 1.0f / (Nx * x_step);
@@ -56,8 +99,7 @@ __global__ void kernel_spectral_lens(cuFloatComplex* output,
         // Store result as complex exponential.
         output[y * Nx + x] = make_cuFloatComplex(cosf(phase), sinf(phase));
     }
-}
-
+} */
 __global__ void
 kernel_circular_mask(float* output, short width, short height, float center_X, float center_Y, float radius)
 {

Backend/cuda_sources/morlet_wavelet.cu

@@ -0,0 +1,44 @@
+#include <cuda_runtime.h>
+#include <cuComplex.h>
+#include <math.h>
+#include "complex_utils.cuh"
+#include "morlet_wavelet.cuh"
+
+#ifndef M_PI
+#define M_PI 3.14159265358979323846
+#endif
+
+
+
+// Kernel to compute Morlet wavelet in frequency domain (simplified version correction terms omitted)
+__global__ void buildMorletKernel(cuComplex* kernel,
+                                  int N, float dt,
+                                  float scale, float omega0)
+{
+    int k = blockIdx.x * blockDim.x + threadIdx.x;
+    if (k >= N) return;
+
+    float freq_hz;
+    if (k < N/2 + 1) {
+        // Positive frequencies (0 to N/2)
+        freq_hz = k / (N * dt);
+    } else {
+        // Negative frequencies (-N/2+1 to -1)
+        freq_hz = (k - N) / (N * dt);
+    }
+
+    float omega = 2.0f * M_PI * freq_hz;
+    float arg = scale * omega - omega0;
+    float value = sqrtf(scale) * expf(-0.5f * arg * arg);
+
+    kernel[k] = make_cuComplex(value, 0.0f);
+}
+
+void createMorletKernel(cuComplex* d_kernel, int N, float dt,
+                        float scale, float omega0, const cudaStream_t stream)
+{
+    uint threads = get_max_threads_1d();
+    uint blocks  = map_blocks_to_problem(N, threads);
+
+    buildMorletKernel<<<blocks, threads, 0, stream>>>(d_kernel, N, dt, scale, omega0);
+}

Backend/cuda_sources/wavelet_transform.cu

@@ -0,0 +1,65 @@
+#include "hardware_limits.hh"
+#include "wavelet_transform.cuh"
+#include "complex_utils.cuh"
+#include "morlet_wavelet.cuh"
+#include "frame_desc.hh"
+#include <iostream>
+using camera::FrameDescriptor;
+
+
+__global__ void mul_conj_kernel(const cuComplex* X,
+                                const cuComplex* Psi,
+                                cuComplex* Y,
+                                int N,
+                                int total)    
+{
+    int idx = blockIdx.x * blockDim.x + threadIdx.x;
+    if (idx >= total) return;
+
+    int frame_res = total / N;
+    int k = idx / frame_res; // freq bin within the transform
+    float xr = X[idx].x, xi = X[idx].y;
+    float pr = Psi[k].x, pi = Psi[k].y;
+
+    // multiply X[idx] by conjugate(Psi[k])
+    Y[idx].x = xr * pr + xi * pi;
+    Y[idx].y = xi * pr - xr * pi; 
+}
+
+__global__ void scale_kernel(cuComplex* data, int total_elements)
+{
+    int k = threadIdx.x + blockIdx.x * blockDim.x;
+    if (k >= total_elements) return;
+    float s = 1.0f / float(total_elements);
+    data[k].x *= s;
+    data[k].y *= s;
+}
+
+void wavelet_transform(cuComplex* output, cuComplex* input, const cufftHandle plan1d,  const FrameDescriptor& fd, int tranformation_size, const cudaStream_t stream, float target_freq)
+{
+    int N =  tranformation_size;
+    float dt = 1.0f; // time step can be set to 1.0 as we work in normalized units (we dont care about the actual time scale here)
+    float omega0 = 6.0f; // central frequency (standard)
+    float scale = omega0 / (2.0f * M_PI * target_freq);
+
+    cuComplex* d_kernel;
+    cudaMalloc(&d_kernel, N * sizeof(cuComplex));
+
+    createMorletKernel(d_kernel, N, dt, scale, omega0, stream);
+
+    cufftExecC2C(plan1d, input, input, CUFFT_FORWARD);
+
+    // Multiply by wavelet in frequency domain
+    int threads = get_max_threads_1d();
+    int total = N * fd.get_frame_res();
+
+    int blocks = map_blocks_to_problem(total, threads);
+
+    mul_conj_kernel<<<blocks, threads, 0, stream>>>(input, d_kernel, output, N, total);
+
+    // Scale by 1/N to normalize FFT
+    scale_kernel<<<blocks, threads, 0,stream>>>(output, N);
+
+    cudaFree(d_kernel);
+}
+

Backend/includes/api/transform_api.hh

@@ -138,6 +138,20 @@ class TransformApi : public IApi
      */
     inline uint get_time_transformation_size() const { return GET_SETTING(TimeTransformationSize); }
 
+
+     /*! \brief Sets the wavelet transform target frequency.
+     *
+     * \return ApiCode the status of the modification: OK, NO_CHANGE or WRONG_COMP_MODE (if in raw mode).
+     */
+    ApiCode set_target_frequency_wavelet(float value) const;
+
+    /*! \brief Returns the wavelet target frequency.
+     *
+     * \return float the wavelet target frequency
+     */
+    inline float get_target_frequency_wavelet() const { return GET_SETTING(TargetFrequencyWavelet); }
+
+
     /*! \brief Modifies the time transformation size. It's the number of frames used for one time transformation. Must
      * be greater than 0.
      *

Backend/includes/compute/fourier_transform.hh

@@ -49,6 +49,7 @@
     holovibes::settings::TimeTransformation,       \
     holovibes::settings::Lambda,                   \
     holovibes::settings::ZDistance,                \
+    holovibes::settings::TargetFrequencyWavelet,   \
     holovibes::settings::PixelSize,                \
     holovibes::settings::SpaceTransformation
 
@@ -167,13 +168,17 @@ class FourierTransform
     /*! \brief Enqueue stft time filtering. */
     void insert_stft();
 
+    /*! \brief Enqueue wavelet transform time filtering. */
+    void insert_wavelet_transform();
+
     /*! \brief Enqueue functions relative to filtering using diagonalization and eigen values.
      *
      * This should eventually replace stft
      */
     void insert_pca();
 
     void insert_ssa_stft();
+    void insert_stft_ssa();
 
     /*!
      * \brief Helper function to get a settings value.

Backend/includes/core/holovibes.hh

@@ -102,6 +102,7 @@
     holovibes::settings::TimeTransformation,                     \
     holovibes::settings::Lambda,                                 \
     holovibes::settings::ZDistance,                              \
+    holovibes::settings::TargetFrequencyWavelet,                 \
     holovibes::settings::ConvolutionMatrix,                      \
     holovibes::settings::DivideConvolutionEnabled,               \
     holovibes::settings::ConvolutionFileName,                    \
@@ -423,6 +424,7 @@ class Holovibes
                                              settings::TimeTransformation{TimeTransformation::NONE},
                                              settings::Lambda{852e-9f},
                                              settings::ZDistance{0.0f},
+                                             settings::TargetFrequencyWavelet{1.5f},
                                              settings::ConvolutionMatrix{std::vector<float>{}},
                                              settings::DivideConvolutionEnabled{false},
                                              settings::ConvolutionFileName{std::string("")},

Backend/includes/enum/enum_time_transformation.hh

@@ -18,15 +18,19 @@ enum class TimeTransformation
     NONE = 0, /*!< No transformation */
     STFT,     /*!< Short-time Fourier transformation */
     PCA,      /*!< Principal component analysis */
-    SSA_STFT  /*!< Self-adaptive Spectrum Analysis Short-time Fourier transformation */
+    SSA_STFT, /*!< Self-adaptive Spectrum Analysis Short-time Fourier transformation */
+    STFT_SSA, /*!< Short-time Fourier transformation Self-adaptive Spectrum Analysis */
+    WAVELET   /*!< Wavelet Transform */
 };
 
 // clang-format off
 SERIALIZE_JSON_ENUM(TimeTransformation, {
     {TimeTransformation::STFT, "STFT"},
     {TimeTransformation::PCA, "PCA"},
     {TimeTransformation::NONE, "NONE"},
-    {TimeTransformation::SSA_STFT, "SSA_STFT"},
+    {TimeTransformation::SSA_STFT, "SSA+STFT"},
+    {TimeTransformation::STFT_SSA, "STFT+SSA"},
+    {TimeTransformation::WAVELET, "WAVELET"},
     {TimeTransformation::NONE, "None"}, // Compat
 
 })

Backend/includes/settings/settings.hh

@@ -171,6 +171,7 @@ DECLARE_SETTING(SpaceTransformation, holovibes::SpaceTransformation);
 DECLARE_SETTING(TimeTransformation, holovibes::TimeTransformation);
 DECLARE_SETTING(Lambda, float);
 DECLARE_SETTING(ZDistance, float);
+DECLARE_SETTING(TargetFrequencyWavelet, float);
 
 DECLARE_SETTING(ConvolutionMatrix, std::vector<float>);
 DECLARE_SETTING(DivideConvolutionEnabled, bool);

Backend/sources/api/transform_api.cc

@@ -117,6 +117,21 @@ ApiCode TransformApi::set_z_distance(float value) const
 
 #pragma region Time Tr.
 
+ApiCode TransformApi::set_target_frequency_wavelet(float value) const
+{
+    NOT_SAME_AND_NOT_RAW(get_target_frequency_wavelet(), value);
+
+    if (value < 0)
+    {
+        LOG_WARN("Target frequency for wavelet cannot be negative. Setting it to 0");
+        value = 0;
+    }
+
+    UPDATE_SETTING(TargetFrequencyWavelet, value);
+
+    return ApiCode::OK;
+}
+
 ApiCode TransformApi::set_time_transformation_size(uint time_transformation_size) const
 {
     NOT_SAME_AND_NOT_RAW(get_time_transformation_size(), time_transformation_size);