초·중·고 정보 과학(적응형 신호 처리·LMS 필터·역위상 간섭)·물리(음파 파동 방정식·음향 임피던스 정합·공진 흡수) 교육과 국방 해군 잠수함사령부 장보고-III(KSS-III) 도산안창호급 잠수함 수중 음향 스텔스 무반사 타일(Acoustic Anechoic Tile) 및 능동 소음 제어(ANC, Active Noise Cancellation): ZeliDesk의 온디바이스 순수 C# 실시간 알베르티히(Alberich) 공진기 점탄성체 음향 반사 계수 수치 모델링, 0-Alloc SIMD 가속 다중 채널 FxLMS(Filtered-x Least Mean Squares) 적응형 역위상 상쇄 제어기 및 능동 소나 반사파 무력화 엔진 아키텍처

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

namespace ZeliDesk.Core.Acoustics
{
    /// <summary>
    /// 알베르티히 점탄성 무반사 타일 층상 음향 임피던스 전송 행렬 수치 모델
    /// </summary>
    [StructLayout(LayoutKind.Sequential, Pack = 16)]
    public struct TileAcousticProperty
    {
        public float Thickness;         // 타일 두께 (m, 통상 0.04m)
        public float Density;           // 매질 밀도 (kg/m^3, 약 1,150)
        public float RealSoundSpeed;    // 종파 음속 실수부 (m/s, 약 1,600)
        public float LossFactorEta;     // 점탄성 손실 계수 eta (약 0.35)
        public float CavityVoidRatio;   // 체적 대비 공동 공극률 (0.0 ~ 0.3)

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public readonly float ComputeEchoReduction(float frequencyHz)
        {
            // 각주파수 omega 및 복소 파수 k* 계산
            float omega = 2.0f * MathF.PI * frequencyHz;
            float complexSoundSpeedSquared = (RealSoundSpeed * RealSoundSpeed) * (1.0f + LossFactorEta * LossFactorEta);
            float effectiveC = MathF.Sqrt(complexSoundSpeedSquared);
            
            // 공동 공진에 따른 음속 저하 및 흡음 피크 모델
            float resonancePeakFreq = 18000.0f; // 18kHz 알베르티히 주 공진점
            float freqDelta = (frequencyHz - resonancePeakFreq) / 5000.0f;
            float resonanceAbsorption = CavityVoidRatio * 2.5f / (1.0f + freqDelta * freqDelta);

            float totalAlpha = (omega * LossFactorEta / effectiveC) + resonanceAbsorption;
            float waveNumberK = omega / effectiveC;

            // 해수 임피던스 (1.54e6) 및 강철 내압 선체 임피던스 (40.0e6)
            float zw = 1.54e6f;
            float zs = 40.0e6f;
            float zTile = Density * effectiveC;

            // 전송선 모델 기반 입력 임피던스 Zin 계산
            float kd = waveNumberK * Thickness;
            float cosKd = MathF.Cos(kd);
            float sinKd = MathF.Sin(kd);

            float numeratorReal = zs * cosKd;
            float numeratorImag = zTile * sinKd;
            float denominatorReal = zTile * cosKd;
            float denominatorImag = zs * sinKd;

            // 음향 반사 계수 크기 |R| 계산
            float zinApprox = zTile * (zs + zTile * MathF.Tan(kd)) / (zTile + zs * MathF.Tan(kd) + 1e-6f);
            float rNumerator = zinApprox - zw;
            float rDenominator = zinApprox + zw;
            float reflectionCoeff = MathF.Abs(rNumerator / rDenominator);

            // 감쇄율 적용
            reflectionCoeff *= MathF.Exp(-totalAlpha * Thickness * 2.0f);
            if (reflectionCoeff < 0.001f) reflectionCoeff = 0.001f;

            // 반사 손실 Echo Reduction (dB)
            return -10.0f * MathF.Log10(reflectionCoeff * reflectionCoeff);
        }
    }

    /// <summary>
    /// 0-Alloc AVX2 SIMD 가속 다중 채널 FxLMS 능동 소음 제어기
    /// 2,000Hz 실시간 선체 진동 역위상 상쇄 엔진
    /// </summary>
    public sealed class MultiChannelFxLmsController
    {
        private const int FilterLength = 64;           // FIR 필터 차수 L
        private const int SecondaryPathLength = 32;    // 2차 경로 추정 필터 차수 M

        private readonly float[] _w;                   // 적응형 제어 필터 가중치 (FilterLength)
        private readonly float[] _sHat;                // 2차 경로 모델 가중치 (SecondaryPathLength)
        private readonly float[] _xBuffer;             // 참조 신호 링 버퍼 (FilterLength * 2)
        private readonly float[] _filteredXBuffer;     // Filtered-x 신호 버퍼 (FilterLength * 2)
        
        private int _bufferIndex;
        private readonly float _muStepSize;

        public MultiChannelFxLmsController(float mu = 0.0025f)
        {
            _muStepSize = mu;
            _w = new float[FilterLength];
            _sHat = new float[SecondaryPathLength];
            _xBuffer = new float[FilterLength * 2];
            _filteredXBuffer = new float[FilterLength * 2];
            _bufferIndex = 0;

            // 2차 경로(DAC -> 액추에이터 -> 선체 전달) 기본 전달 특성 초기화 (지연 감쇄 모델)
            for (int i = 0; i < SecondaryPathLength; i++)
            {
                _sHat[i] = MathF.Exp(-0.15f * i) * MathF.Sin(0.4f * MathF.PI * i);
            }
        }

        /// <summary>
        /// 1 샘플 단위 실시간 FxLMS 갱신 (0-Alloc 실행, 2,000Hz 주기 호출)
        /// </summary>
        /// <param name="primaryNoiseRef">선체 내부 추진기 진동 기준 센서 신호</param>
        /// <param name="errorHydrophoneMic">선체 외판 잔여 음압 에러 마이크 신호</param>
        /// <param name="outAntiNoiseControl">선체 프레임 액추에이터 역위상 구동 사령</param>
        /// <returns>소음 감쇄량 (dB)</returns>
        [MethodImpl(MethodImplOptions.AggressiveOptimization)]
        public float ProcessSample(float primaryNoiseRef, float errorHydrophoneMic, out float outAntiNoiseControl)
        {
            // 1. 참조 신호 링 버퍼 주입 (더블 버퍼 미러링으로 분기 제거)
            _bufferIndex = (_bufferIndex - 1) & (FilterLength - 1);
            _xBuffer[_bufferIndex] = primaryNoiseRef;
            _xBuffer[_bufferIndex + FilterLength] = primaryNoiseRef;

            ReadOnlySpan<float> xSpan = _xBuffer.AsSpan(_bufferIndex, FilterLength);

            // 2. Filtered-x 신호 생성: x'(n) = s_hat * x(n)
            float filteredX = 0.0f;
            for (int j = 0; j < SecondaryPathLength; j++)
            {
                filteredX += _sHat[j] * _xBuffer[_bufferIndex + j];
            }

            _filteredXBuffer[_bufferIndex] = filteredX;
            _filteredXBuffer[_bufferIndex + FilterLength] = filteredX;
            ReadOnlySpan<float> fXSpan = _filteredXBuffer.AsSpan(_bufferIndex, FilterLength);

            // 3. 역위상 제어 출력 합성: y(n) = w^T * x(n) (AVX2 SIMD 내적)
            float y = SimdDotProduct(xSpan, _w.AsSpan());
            outAntiNoiseControl = -y; // 역위상 상쇄 제어 신호

            // 4. 필터 가중치 적응 갱신: w(n+1) = w(n) + mu * e(n) * x'(n) (AVX2 FMA 가속)
            float error = errorHydrophoneMic;
            float muError = _muStepSize * error;

            SimdWeightUpdate(_w.AsSpan(), fXSpan, muError);

            // 5. 실시간 소음 감쇄 지표 산출
            float primaryMagnitude = MathF.Abs(primaryNoiseRef) + 1e-6f;
            float residualMagnitude = MathF.Abs(error) + 1e-6f;
            float attenuationDb = 20.0f * MathF.Log10(primaryMagnitude / residualMagnitude);
            return MathF.Max(0.0f, attenuationDb);
        }

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        private static float SimdDotProduct(ReadOnlySpan<float> a, ReadOnlySpan<float> b)
        {
            if (Avx2.IsSupported)
            {
                var sumVector = Vector256<float>.Zero;
                int i = 0;
                for (; i <= FilterLength - 8; i += 8)
                {
                    var va = Vector256.Create(a.Slice(i, 8));
                    var vb = Vector256.Create(b.Slice(i, 8));
                    sumVector = Fma.MultiplyAdd(va, vb, sumVector);
                }

                // 수평 합산 (Horizontal Add)
                var v128 = Sse.Add(sumVector.GetLower(), sumVector.GetUpper());
                v128 = Sse3.HorizontalAdd(v128, v128);
                v128 = Sse3.HorizontalAdd(v128, v128);
                return v128.ToScalar();
            }
            else
            {
                float sum = 0.0f;
                for (int i = 0; i < FilterLength; i++)
                {
                    sum += a[i] * b[i];
                }
                return sum;
            }
        }

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        private static void SimdWeightUpdate(Span<float> w, ReadOnlySpan<float> fx, float muE)
        {
            if (Avx2.IsSupported)
            {
                var vMuE = Vector256.Create(muE);
                for (int i = 0; i <= FilterLength - 8; i += 8)
                {
                    var vw = Vector256.Create(w.Slice(i, 8));
                    var vfx = Vector256.Create(fx.Slice(i, 8));
                    var vNewW = Fma.MultiplyAdd(vMuE, vfx, vw);
                    vNewW.CopyTo(w.Slice(i, 8));
                }
            }
            else
            {
                for (int i = 0; i < FilterLength; i++)
                {
                    w[i] += muE * fx[i];
                }
            }
        }
    }
}