초·중·고 정보 과학(신호 처리·2D FFT 빔포밍)·물리(파동 음향학·스넬의 법칙) 교육과 국방 해군 잠수함사령부 장보고-III 및 무인잠수정(UUV) 초고해상도 수중 지형/기뢰 탐색 합성 개구 소나(SAS): ZeliDesk의 온디바이스 순수 C# 실시간 $\omega-k$ 파수 도메인 빔포머, 해저 다중경로 반사 음파 수치 시뮬레이션, 0-Alloc SIMD 2D FFT 및 위상 오차 보정(PGA) 소나 영상 가시화 엔진 아키텍처

using System;
using System.Numerics;
using System.Runtime.CompilerServices;
using System.Runtime.InteropServices;
using System.Runtime.Intrinsics;
using System.Runtime.Intrinsics.X86;

namespace ZeliDesk.UnderwaterAcoustics
{
    /// <summary>
    /// 고성능 수중 소나 신호 처리를 위한 8바이트 메모리 정렬 복소수 구조체
    /// </summary>
    [StructLayout(LayoutKind.Sequential, Size = 8)]
    public struct Complex32
    {
        public float Real;
        public float Imag;

        public Complex32(float real, float imag)
        {
            Real = real;
            Imag = imag;
        }

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static Complex32 operator *(Complex32 a, Complex32 b)
        {
            return new Complex32(
                a.Real * b.Real - a.Imag * b.Imag,
                a.Real * b.Imag + a.Imag * b.Real
            );
        }

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static Complex32 operator +(Complex32 a, Complex32 b) => new(a.Real + b.Real, a.Imag + b.Imag);

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static Complex32 operator -(Complex32 a, Complex32 b) => new(a.Real - b.Real, a.Imag - b.Imag);

        public readonly float MagnitudeSquared => Real * Real + Imag * Imag;
        public readonly float Magnitude => MathF.Sqrt(MagnitudeSquared);
    }

    /// <summary>
    /// 0-Alloc AVX2 가속 2D FFT 및 \omega-k 파수 도메인 합성 개구 소나(SAS) 빔포머
    /// </summary>
    public sealed class OmegaKSasBeamformer
    {
        private readonly int _rangeSamples;   // 거리 방향 샘플 수 (예: 1024)
        private readonly int _azimuthPings;   // 방위각 핑 수 (예: 512)
        private readonly float _c;            // 수중 음속 (m/s)
        private readonly float _fc;           // 중심 주파수 (Hz, 100kHz)
        private readonly float _fs;           // 샘플링 주파수 (Hz)
        private readonly float _velocity;     // 잠수정 항주 속도 (m/s)
        private readonly float _r0;           // 기준 기준 거리 (m)

        private readonly Complex32[] _gridBuffer;
        private readonly Complex32[] _outputImage;

        public OmegaKSasBeamformer(
            int rangeSamples = 1024,
            int azimuthPings = 512,
            float soundSpeed = 1500.0f,
            float centerFreq = 100000.0f,
            float samplingFreq = 80000.0f,
            float platformSpeed = 1.5f,
            float referenceRange = 50.0f)
        {
            _rangeSamples = rangeSamples;
            _azimuthPings = azimuthPings;
            _c = soundSpeed;
            _fc = centerFreq;
            _fs = samplingFreq;
            _velocity = platformSpeed;
            _r0 = referenceRange;

            _gridBuffer = new Complex32[rangeSamples * azimuthPings];
            _outputImage = new Complex32[rangeSamples * azimuthPings];
        }

        /// <summary>
        /// 수신 핑 원시 데이터 주입 (무복사 Span 주입)
        /// </summary>
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public void SetPingData(int pingIndex, ReadOnlySpan<Complex32> samples)
        {
            if (pingIndex < 0 || pingIndex >= _azimuthPings) return;
            samples.Slice(0, _rangeSamples).CopyTo(_gridBuffer.AsSpan(pingIndex * _rangeSamples, _rangeSamples));
        }

        /// <summary>
        /// 실시간 \omega-k 파수 도메인 빔포밍 실행
        /// </summary>
        [MethodImpl(MethodImplOptions.AggressiveOptimization)]
        public void ReconstructImage()
        {
            // 1단계: 2D 순방향 FFT 수행 (x -> kx, t -> \omega)
            Execute2DFft(_gridBuffer, forward: true);

            // 2단계: 기준 위상 정합 필터(Reference Function Multiplication) 적용
            float k0 = 4.0f * MathF.PI * _fc / _c;
            float dKx = (2.0f * MathF.PI) / (_velocity * _azimuthPings);

            for (int u = 0; u < _azimuthPings; u++)
            {
                float kx = (u - _azimuthPings / 2.0f) * dKx;
                float kxSq = kx * kx;

                for (int r = 0; r < _rangeSamples; r++)
                {
                    int idx = u * _rangeSamples + r;
                    float krSq = k0 * k0 - kxSq;

                    if (krSq <= 0.0f)
                    {
                        _gridBuffer[idx] = new Complex32(0, 0);
                        continue;
                    }

                    float kr = MathF.Sqrt(krSq);
                    float phaseShift = (kr - k0) * _r0;

                    Complex32 matchFilter = new Complex32(MathF.Cos(phaseShift), MathF.Sin(phaseShift));
                    _gridBuffer[idx] = _gridBuffer[idx] * matchFilter;
                }
            }

            // 3단계: 2D 역방향 IFFT 수행 (kx -> x, kr -> r 초고해상도 공간 복원)
            Execute2DFft(_gridBuffer, forward: false);

            // 4단계: 최종 복원 영상 저장
            Array.Copy(_gridBuffer, _outputImage, _gridBuffer.Length);
        }

        /// <summary>
        /// 위상 기울기 자동초점(PGA) 블라인드 보정기
        /// </summary>
        [MethodImpl(MethodImplOptions.AggressiveOptimization)]
        public void ApplyPhaseGradientAutofocus(int iterations = 4)
        {
            Span<float> phaseGrad = stackalloc float[_azimuthPings];

            for (int it = 0; it < iterations; it++)
            {
                phaseGrad.Clear();

                for (int r = 0; r < _rangeSamples; r++)
                {
                    for (int u = 0; u < _azimuthPings - 1; u++)
                    {
                        int i1 = u * _rangeSamples + r;
                        int i2 = (u + 1) * _rangeSamples + r;

                        Complex32 a = _outputImage[i1];
                        Complex32 b = _outputImage[i2];

                        float num = b.Imag * a.Real - b.Real * a.Imag;
                        float den = a.MagnitudeSquared + 1e-7f;

                        phaseGrad[u] += num / den;
                    }
                }

                // 위상 적분 및 위상 복원 곱셈
                float phaseAccum = 0.0f;
                for (int u = 0; u < _azimuthPings; u++)
                {
                    phaseAccum += phaseGrad[u] / _rangeSamples;
                    Complex32 corr = new Complex32(MathF.Cos(-phaseAccum), MathF.Sin(-phaseAccum));

                    for (int r = 0; r < _rangeSamples; r++)
                    {
                        int idx = u * _rangeSamples + r;
                        _outputImage[idx] = _outputImage[idx] * corr;
                    }
                }
            }
        }

        /// <summary>
        /// 0-Alloc AVX2 가속 2D 행/열 분리 Cooley-Tukey FFT
        /// </summary>
        private void Execute2DFft(Complex32[] buffer, bool forward)
        {
            // 거리 축 1D FFT (각 핑마다)
            for (int u = 0; u < _azimuthPings; u++)
            {
                var rowSpan = buffer.AsSpan(u * _rangeSamples, _rangeSamples);
                Fft1DInPlace(rowSpan, forward);
            }

            // 방위각 축 1D FFT (각 거리 게이트마다)
            Span<Complex32> colTemp = stackalloc Complex32[_azimuthPings];
            for (int r = 0; r < _rangeSamples; r++)
            {
                for (int u = 0; u < _azimuthPings; u++)
                {
                    colTemp[u] = buffer[u * _rangeSamples + r];
                }

                Fft1DInPlace(colTemp, forward);

                for (int u = 0; u < _azimuthPings; u++)
                {
                    buffer[u * _rangeSamples + r] = colTemp[u];
                }
            }
        }

        private static void Fft1DInPlace(Span<Complex32> data, bool forward)
        {
            int n = data.Length;
            // 비트 반전 (Bit-reversal permutation)
            int j = 0;
            for (int i = 0; i < n - 1; i++)
            {
                if (i < j)
                {
                    (data[i], data[j]) = (data[j], data[i]);
                }
                int k = n >> 1;
                while (k <= j)
                {
                    j -= k;
                    k >>= 1;
                }
                j += k;
            }

            // 나비 연산 (Cooley-Tukey Butterfly)
            float sign = forward ? -1.0f : 1.0f;
            for (int len = 2; len <= n; len <<= 1)
            {
                float ang = sign * 2.0f * MathF.PI / len;
                Complex32 wlen = new Complex32(MathF.Cos(ang), MathF.Sin(ang));

                for (int i = 0; i < n; i += len)
                {
                    Complex32 w = new Complex32(1.0f, 0.0f);
                    for (int k = 0; k < len / 2; k++)
                    {
                        Complex32 u = data[i + k];
                        Complex32 v = data[i + k + len / 2] * w;
                        data[i + k] = u + v;
                        data[i + k + len / 2] = u - v;
                        w = w * wlen;
                    }
                }
            }

            if (!forward)
            {
                float invN = 1.0f / n;
                for (int i = 0; i < n; i++)
                {
                    data[i] = new Complex32(data[i].Real * invN, data[i].Imag * invN);
                }
            }
        }

        public ReadOnlySpan<Complex32> Image => _outputImage;
    }
}