unit EngineeringLib.Signal;

{-----------------------------------------------------------------------------
 EngineeringLib.Signal

 Digital signal processing toolkit.

 Provides:
   Basic Filtering
     MovingAverage    — simple sliding-window average

   Windowing
     GenerateWindow   — Rectangular, Hamming, Hann, Blackman
     ApplyWindow      — element-wise multiply signal × window

   FFT / Spectral Analysis
     FFT              — in-place Cooley-Tukey radix-2 DIT (N must be power of 2)
     IFFT             — inverse FFT
     CalculateFFT     — convenience wrapper: real input → complex output arrays
     CalculateIFFT    — inverse: complex arrays → real output
     CalculateFFTMagnitudePhase — complete magnitude and phase spectra

   FIR Filter Design (windowed-sinc)
     DesignFIRLowPass   — low-pass FIR coefficients
     DesignFIRHighPass  — high-pass FIR coefficients
     DesignFIRBandPass  — band-pass FIR coefficients
     DesignFIRBandStop  — band-stop (notch) FIR coefficients
     ApplyFIRFilter     — convolve signal with FIR coefficients

   Signal Properties
     SignalPower      — mean square
     SignalEnergy     — sum of squares
     RootMeanSquare   — RMS value
-----------------------------------------------------------------------------}

{$mode ObjFPC}{$H+}

interface

uses
  Classes, SysUtils, Math, MathBase.SharedTypes, MathBase.Complex,
  EngineeringLib.Common;

type
  { Keep the familiar EngineeringLib.Signal name while sharing the project-wide
    dynamic-array type from MathBase. }
  TDoubleArray = MathBase.SharedTypes.TDoubleArray;

  TWindowType = (wtRectangular, wtHamming, wtHann, wtBlackman);

  TSignalKit = class
  public
    { --- Basic Filtering --- }
    class function MovingAverage(const InputSignal: TDoubleArray; WindowSize: Integer): TDoubleArray; static;

    { --- Windowing --- }
    class function GenerateWindow(WindowType: TWindowType; Size: Integer): TDoubleArray; static;
    class function ApplyWindow(const InputSignal, Window: TDoubleArray): TDoubleArray; static;

    { --- FFT / Spectral Analysis ---
      FFT and IFFT operate on complex data split into separate Real/Imag arrays.
      N must be a power of 2. }

    { In-place Cooley-Tukey radix-2 DIT FFT.
      RealPart and ImagPart must both have length N (a power of 2).
      Set Inverse=True for IFFT (includes the 1/N scaling). }
    class procedure FFT(var RealPart, ImagPart: TDoubleArray; Inverse: Boolean = False); static;

    { In-place complex-array FFT. This is a source-compatible companion to
      the split-array overload above and has the same power-of-two contract. }
    class procedure FFT(var Data: TComplexArray; Inverse: Boolean = False); overload; static;

    { Convenience wrapper: real-valued input signal → complete complex spectrum.
      Zero-pads to the next power-of-2 length automatically; it never truncates.
      Empty input produces two empty output arrays. }
    class procedure CalculateFFT(
      const InputSignal: TDoubleArray;
      out OutRealPart, OutImagPart: TDoubleArray); overload; static;

    { Convenience wrapper: real-valued input signal → complete complex
      spectrum. It has the same zero-padding and empty-input behavior as the
      split-array overload. }
    class procedure CalculateFFT(const InputSignal: TDoubleArray;
      out OutputSpectrum: TComplexArray); overload; static;

    { Inverse FFT: complete complex spectrum → real-valued signal.
      The real and imaginary arrays must have equal, power-of-2 lengths.
      Two empty input arrays produce an empty output array. }
    class procedure CalculateIFFT(const InRealPart, InImagPart: TDoubleArray; out OutputSignal: TDoubleArray); static;

    { Inverse FFT from a complete complex spectrum → real-valued signal. }
    class procedure CalculateIFFT(const InputSpectrum: TComplexArray;
      out OutputSignal: TDoubleArray); overload; static;

    { Complete N-bin magnitude and phase spectra. }
    class procedure CalculateFFTMagnitudePhase(
      const InputSignal: TDoubleArray;
      out Magnitude, Phase: TDoubleArray); static;

    { --- FIR Filter Design (windowed-sinc) ---
      CutoffFreq is normalised: 0 < fc < 0.5 (where 0.5 = Nyquist).
      Order is the filter order (number of coefficients = Order+1). It must be
      at least 2; odd orders are incremented to the next even order so every
      design remains a symmetric, odd-length linear-phase FIR.
      WindowType selects the window used to taper the ideal sinc kernel. }

    class function DesignFIRLowPass(
      CutoffFreq: Double;
      Order: Integer;
      WindowType: TWindowType = wtHamming): TDoubleArray; static;

    class function DesignFIRHighPass(
      CutoffFreq: Double;
      Order: Integer;
      WindowType: TWindowType = wtHamming): TDoubleArray; static;

    class function DesignFIRBandPass(
      LowCutoff, HighCutoff: Double;
      Order: Integer;
      WindowType: TWindowType = wtHamming): TDoubleArray; static;

    class function DesignFIRBandStop(
      LowCutoff, HighCutoff: Double;
      Order: Integer;
      WindowType: TWindowType = wtHamming): TDoubleArray; static;

    { Convolve Signal with FIR Coefficients (linear/direct-form convolution).
      Output length = Length(Signal) + Length(Coeffs) - 1. }
    class function ApplyFIRFilter(const Signal, Coeffs: TDoubleArray): TDoubleArray; static;

    { --- Signal Properties --- }
    class function SignalPower(const InputSignal: TDoubleArray): Double; static;
    class function SignalEnergy(const InputSignal: TDoubleArray): Double; static;
    class function RootMeanSquare(const InputSignal: TDoubleArray): Double; static;
  end;

implementation

{ ---------------------------------------------------------------------------
  Helpers
  --------------------------------------------------------------------------- }

{ Return smallest power of 2 >= N }
function NextPow2(N: Integer): Integer;
begin
  Result := 1;
  while Result < N do Result := Result shl 1;
end;

{ Bit-reversal permutation for Cooley-Tukey }
procedure BitReverse(var Re, Im: TDoubleArray);
var
  N, I, J, K: Integer;
  Tmp: Double;
begin
  N := Length(Re);
  J := 0;
  for I := 1 to N - 1 do
  begin
    K := N shr 1;
    while J >= K do begin J := J - K; K := K shr 1; end;
    J := J + K;
    if I < J then
    begin
      Tmp := Re[I]; Re[I] := Re[J]; Re[J] := Tmp;
      Tmp := Im[I]; Im[I] := Im[J]; Im[J] := Tmp;
    end;
  end;
end;

{ ---------------------------------------------------------------------------
  TSignalKit — Moving Average
  --------------------------------------------------------------------------- }

class function TSignalKit.MovingAverage(const InputSignal: TDoubleArray; WindowSize: Integer): TDoubleArray;
var
  N, I: Integer;
  Sum: Double;
begin
  if InputSignal = nil then
    Exit(nil);

  N := Length(InputSignal);
  if (WindowSize <= 0) or (WindowSize > N) then
    raise ESignalError.Create('Invalid window size for moving average.');

  SetLength(Result, N);

  Sum := 0.0;
  for I := 0 to WindowSize - 1 do
    Sum := Sum + InputSignal[I];
  Result[WindowSize - 1] := Sum / WindowSize;

  for I := WindowSize to N - 1 do
  begin
    Sum := Sum - InputSignal[I - WindowSize] + InputSignal[I];
    Result[I] := Sum / WindowSize;
  end;

  for I := 0 to WindowSize - 2 do
    Result[I] := Result[WindowSize - 1];
end;

{ ---------------------------------------------------------------------------
  TSignalKit — Windowing
  --------------------------------------------------------------------------- }

class function TSignalKit.GenerateWindow(WindowType: TWindowType; Size: Integer): TDoubleArray;
var
  I: Integer;
  N_1: Double;
const
  Blackman_a0 = 0.42;
  Blackman_a1 = 0.5;
  Blackman_a2 = 0.08;
begin
  Result := nil;
  if Size <= 0 then
    raise ESignalError.Create('Window size must be positive.');
  SetLength(Result, Size);
  N_1 := Size - 1;
  if N_1 = 0 then begin Result[0] := 1.0; Exit; end;

  case WindowType of
    wtRectangular:
      for I := 0 to Size - 1 do Result[I] := 1.0;
    wtHamming:
      for I := 0 to Size - 1 do
        Result[I] := 0.54 - 0.46 * Cos(2 * Pi * I / N_1);
    wtHann:
      for I := 0 to Size - 1 do
        Result[I] := 0.5 * (1 - Cos(2 * Pi * I / N_1));
    wtBlackman:
      for I := 0 to Size - 1 do
        Result[I] := Blackman_a0
                   - Blackman_a1 * Cos(2 * Pi * I / N_1)
                   + Blackman_a2 * Cos(4 * Pi * I / N_1);
  else
    raise ESignalError.Create('Unknown window type specified.');
  end;
end;

class function TSignalKit.ApplyWindow(const InputSignal, Window: TDoubleArray): TDoubleArray;
var
  N, M, I: Integer;
begin
  N := Length(InputSignal);
  M := Length(Window);
  if N <> M then
    raise ESignalError.Create('Input signal and window must have the same length.');
  if N = 0 then Exit(nil);
  SetLength(Result, N);
  for I := 0 to N - 1 do
    Result[I] := InputSignal[I] * Window[I];
end;

{ ---------------------------------------------------------------------------
  TSignalKit — FFT  (Cooley-Tukey radix-2 DIT)
  --------------------------------------------------------------------------- }

class procedure TSignalKit.FFT(var RealPart, ImagPart: TDoubleArray; Inverse: Boolean);
var
  N, I: Integer;
  Data: TComplexArray;
begin
  N := Length(RealPart);
  if N <> Length(ImagPart) then
    raise ESignalError.Create('FFT: RealPart and ImagPart must have the same length.');
  Data := nil;
  SetLength(Data, N);
  for I := 0 to N - 1 do
    Data[I] := TComplex.Create(RealPart[I], ImagPart[I]);
  FFT(Data, Inverse);
  for I := 0 to N - 1 do
  begin
    RealPart[I] := Data[I].Re;
    ImagPart[I] := Data[I].Im;
  end;
end;

procedure BitReverseComplex(var Data: TComplexArray);
var
  N, I, J, K: Integer;
  Temp: TComplex;
begin
  N := Length(Data);
  J := 0;
  for I := 1 to N - 1 do
  begin
    K := N shr 1;
    while J >= K do begin J := J - K; K := K shr 1; end;
    J := J + K;
    if I < J then
    begin
      Temp := Data[I]; Data[I] := Data[J]; Data[J] := Temp;
    end;
  end;
end;

class procedure TSignalKit.FFT(var Data: TComplexArray; Inverse: Boolean);
var
  N, Len, I, J, K: Integer;
  Angle, WR, WI, Ur, Ui, TR, TI: Double;
  Sign: Double;
begin
  N := Length(Data);
  if N <= 1 then Exit;
  if (N and (N - 1)) <> 0 then
    raise ESignalError.Create('FFT: length must be a power of 2.');
  BitReverseComplex(Data);
  Sign := IfThen(Inverse, 1.0, -1.0);

  Len := 1;
  while Len < N do
  begin
    Len := Len shl 1;
    Angle := Sign * 2 * Pi / Len;
    WR := Cos(Angle);
    WI := Sin(Angle);
    I := 0;
    while I < N do
    begin
      Ur := 1.0; Ui := 0.0;
      for J := 0 to (Len shr 1) - 1 do
      begin
        K := I + J + (Len shr 1);
        TR := Ur * Data[K].Re - Ui * Data[K].Im;
        TI := Ur * Data[K].Im + Ui * Data[K].Re;
        Data[K].Re := Data[I + J].Re - TR;
        Data[K].Im := Data[I + J].Im - TI;
        Data[I + J].Re := Data[I + J].Re + TR;
        Data[I + J].Im := Data[I + J].Im + TI;
        TR := Ur * WR - Ui * WI;
        Ui := Ur * WI + Ui * WR;
        Ur := TR;
      end;
      I := I + Len;
    end;
  end;
  if Inverse then
    for I := 0 to N - 1 do
    begin
      Data[I].Re := Data[I].Re / N;
      Data[I].Im := Data[I].Im / N;
    end;
end;

class procedure TSignalKit.CalculateFFT(const InputSignal: TDoubleArray; out OutRealPart, OutImagPart: TDoubleArray);
var
  I: Integer;
  Spectrum: TComplexArray;
begin
  Spectrum := nil;
  CalculateFFT(InputSignal, Spectrum);
  OutRealPart := nil;
  OutImagPart := nil;
  SetLength(OutRealPart, Length(Spectrum));
  SetLength(OutImagPart, Length(Spectrum));
  for I := 0 to High(Spectrum) do
  begin
    OutRealPart[I] := Spectrum[I].Re;
    OutImagPart[I] := Spectrum[I].Im;
  end;
end;

class procedure TSignalKit.CalculateFFT(const InputSignal: TDoubleArray;
  out OutputSpectrum: TComplexArray);
var
  N, I: Integer;
begin
  if Length(InputSignal) = 0 then
  begin
    OutputSpectrum := nil;
    Exit;
  end;
  N := NextPow2(Length(InputSignal));
  OutputSpectrum := nil;
  SetLength(OutputSpectrum, N);
  for I := 0 to High(InputSignal) do
    OutputSpectrum[I] := TComplex.Create(InputSignal[I], 0.0);
  for I := Length(InputSignal) to N - 1 do
    OutputSpectrum[I] := TComplex.Zero;
  FFT(OutputSpectrum);
end;

class procedure TSignalKit.CalculateIFFT(const InRealPart, InImagPart: TDoubleArray; out OutputSignal: TDoubleArray);
var
  N, I: Integer;
  Spectrum: TComplexArray;
begin
  N := Length(InRealPart);
  if N <> Length(InImagPart) then
    raise ESignalError.Create(
      'CalculateIFFT: real and imaginary parts must have the same length.');
  if N = 0 then
  begin
    OutputSignal := nil;
    Exit;
  end;

  Spectrum := nil;
  SetLength(Spectrum, N);
  for I := 0 to N - 1 do
    Spectrum[I] := TComplex.Create(InRealPart[I], InImagPart[I]);
  CalculateIFFT(Spectrum, OutputSignal);
end;

class procedure TSignalKit.CalculateIFFT(const InputSpectrum: TComplexArray;
  out OutputSignal: TDoubleArray);
var
  I: Integer;
  Work: TComplexArray;
begin
  if Length(InputSpectrum) = 0 then
  begin
    OutputSignal := nil;
    Exit;
  end;
  Work := nil;
  SetLength(Work, Length(InputSpectrum));
  for I := 0 to High(InputSpectrum) do
    Work[I] := InputSpectrum[I];
  FFT(Work, True);
  SetLength(OutputSignal, Length(Work));
  for I := 0 to High(Work) do
    OutputSignal[I] := Work[I].Re;
end;

class procedure TSignalKit.CalculateFFTMagnitudePhase(
  const InputSignal: TDoubleArray;
  out Magnitude, Phase: TDoubleArray);
var
  Re, Im: TDoubleArray;
  N, I: Integer;
begin
  CalculateFFT(InputSignal, Re, Im);
  N := Length(Re);
  SetLength(Magnitude, N);
  SetLength(Phase, N);
  for I := 0 to N - 1 do
  begin
    Magnitude[I] := Sqrt(Re[I]*Re[I] + Im[I]*Im[I]);
    Phase[I]     := ArcTan2(Im[I], Re[I]);
  end;
end;

{ ---------------------------------------------------------------------------
  TSignalKit — FIR Filter Design (windowed-sinc)
  --------------------------------------------------------------------------- }

{ Ideal low-pass sinc kernel, centred at tap M/2, normalised cutoff fc }
function SincKernel(Order: Integer; CutoffFreq: Double; WinType: TWindowType): TDoubleArray;
var
  M, I: Integer;
  Fc2Pi, N0, Sinc: Double;
  Win: TDoubleArray;
begin
  Result := nil;
  M     := Order;
  Fc2Pi := 2 * Pi * CutoffFreq;
  SetLength(Result, M + 1);
  Win := TSignalKit.GenerateWindow(WinType, M + 1);
  for I := 0 to M do
  begin
    N0 := I - M / 2;
    if N0 = 0 then
      Sinc := Fc2Pi / Pi      { limit: sin(x)/x → 1, so 2*fc }
    else
      Sinc := Sin(Fc2Pi * N0) / (Pi * N0);
    Result[I] := Sinc * Win[I];
  end;
end;

{ Normalise coefficients so DC gain = 1 }
procedure NormaliseCoeffs(var C: TDoubleArray);
var
  Sum: Double;
  I: Integer;
begin
  Sum := 0;
  for I := 0 to High(C) do Sum := Sum + C[I];
  if Abs(Sum) > 1E-300 then
    for I := 0 to High(C) do C[I] := C[I] / Sum;
end;

class function TSignalKit.DesignFIRLowPass(CutoffFreq: Double; Order: Integer; WindowType: TWindowType): TDoubleArray;
begin
  if (CutoffFreq <= 0) or (CutoffFreq >= 0.5) then
    raise ESignalError.Create('DesignFIRLowPass: CutoffFreq must be in (0, 0.5).');
  if Order < 2 then
    raise ESignalError.Create('DesignFIRLowPass: Order must be at least 2.');
  if Odd(Order) then Inc(Order);
  Result := SincKernel(Order, CutoffFreq, WindowType);
  NormaliseCoeffs(Result);
end;

class function TSignalKit.DesignFIRHighPass(CutoffFreq: Double; Order: Integer; WindowType: TWindowType): TDoubleArray;
var
  LP: TDoubleArray;
  I: Integer;
begin
  Result := nil;
  { High-pass = spectral inversion of low-pass }
  LP := DesignFIRLowPass(CutoffFreq, Order, WindowType);
  SetLength(Result, Length(LP));
  for I := 0 to High(LP) do
    Result[I] := -LP[I];
  Result[High(Result) div 2] := Result[High(Result) div 2] + 1.0;
end;

class function TSignalKit.DesignFIRBandPass(
  LowCutoff, HighCutoff: Double;
  Order: Integer;
  WindowType: TWindowType): TDoubleArray;
var
  LP_Lo, LP_Hi: TDoubleArray;
  I: Integer;
begin
  Result := nil;
  { Band-pass = low-pass(fc_hi) − low-pass(fc_lo) }
  if LowCutoff >= HighCutoff then
    raise ESignalError.Create('DesignFIRBandPass: LowCutoff must be < HighCutoff.');
  LP_Lo := DesignFIRLowPass(LowCutoff,  Order, WindowType);
  LP_Hi := DesignFIRLowPass(HighCutoff, Order, WindowType);
  SetLength(Result, Length(LP_Lo));
  for I := 0 to High(Result) do
    Result[I] := LP_Hi[I] - LP_Lo[I];
end;

class function TSignalKit.DesignFIRBandStop(
  LowCutoff, HighCutoff: Double;
  Order: Integer;
  WindowType: TWindowType): TDoubleArray;
var
  BP: TDoubleArray;
  I: Integer;
begin
  Result := nil;
  { Band-stop = 1 − band-pass (spectral inversion) }
  BP := DesignFIRBandPass(LowCutoff, HighCutoff, Order, WindowType);
  SetLength(Result, Length(BP));
  for I := 0 to High(BP) do
    Result[I] := -BP[I];
  Result[High(Result) div 2] := Result[High(Result) div 2] + 1.0;
end;

class function TSignalKit.ApplyFIRFilter(const Signal, Coeffs: TDoubleArray): TDoubleArray;
var
  NS, NC, OutLen, I, J, K: Integer;
begin
  NS     := Length(Signal);
  NC     := Length(Coeffs);
  if (NS = 0) or (NC = 0) then Exit(nil);
  OutLen := NS + NC - 1;
  SetLength(Result, OutLen);
  for I := 0 to OutLen - 1 do
  begin
    Result[I] := 0;
    for J := 0 to NC - 1 do
    begin
      K := I - J;
      if (K >= 0) and (K < NS) then
        Result[I] := Result[I] + Coeffs[J] * Signal[K];
    end;
  end;
end;

{ ---------------------------------------------------------------------------
  TSignalKit — Signal Properties
  --------------------------------------------------------------------------- }

class function TSignalKit.SignalPower(const InputSignal: TDoubleArray): Double;
var
  SumSq: Double;
  I, N: Integer;
begin
  N := Length(InputSignal);
  if N = 0 then Exit(0.0);
  SumSq := 0.0;
  for I := 0 to N - 1 do SumSq := SumSq + Sqr(InputSignal[I]);
  Result := SumSq / N;
end;

class function TSignalKit.SignalEnergy(const InputSignal: TDoubleArray): Double;
var
  SumSq: Double;
  I, N: Integer;
begin
  N := Length(InputSignal);
  if N = 0 then Exit(0.0);
  SumSq := 0.0;
  for I := 0 to N - 1 do SumSq := SumSq + Sqr(InputSignal[I]);
  Result := SumSq;
end;

class function TSignalKit.RootMeanSquare(const InputSignal: TDoubleArray): Double;
begin
  Result := Sqrt(SignalPower(InputSignal));
end;

end.