Use WDSP analyzer for wideband

WidebandView.pas

@@ -4,10 +4,10 @@ unit WidebandView;
   WidebandView.pas - raw ADC wideband spectrum pane.
 
   The network layer feeds 16-bit ADC samples collected from Protocol V4
-  wideband packets. This view windows the frame, runs an in-process real FFT
-  and renders a compact Thetis-style wideband panadapter. CPU paint uses a
-  bitmap; when MainForm creates a TOpenGLControl the same spectrum is rendered
-  as GL primitives.
+  wideband packets. This view feeds them to a dedicated WDSP analyzer and
+  renders a compact Thetis-style wideband panadapter. CPU paint uses a bitmap;
+  when MainForm creates a TOpenGLControl the same spectrum is rendered as GL
+  primitives.
 }
 
 {$IFDEF FPC}
@@ -18,7 +18,7 @@ interface
 
 uses
   Classes, SysUtils, Math, Graphics, Controls, ExtCtrls,
-  OpenGLContext, GL, AppTheme;
+  OpenGLContext, GL, AppTheme, WDSP;
 
 type
   TWidebandView = class
@@ -39,13 +39,22 @@ type
     FSourceStartHz: Double;
     FSourceEndHz: Double;
     FMarkerHz: Double;
+    FWDSPAnalyzerOpen: Boolean;
+    FWDSPAnalyzerConfigured: Boolean;
+    FWDSPAnalyzerID: Integer;
+    FWDSPFlp: array[0..0] of Integer;
+    FWDSPIn: array of Double;
+    FWDSPPixels: array of Single;
     procedure EnsureBitmap(W, H: Integer);
-    procedure ComputeSpectrum(const Samples: array of SmallInt; Count: Integer);
+    function EnsureWDSPAnalyzer: Boolean;
+    function ComputeSpectrumWDSP(const Samples: array of SmallInt;
+      Count: Integer): Boolean;
+    procedure CloseWDSPAnalyzer;
     procedure DrawCPU(W, H: Integer);
     procedure PaintCPU(C: TCanvas);
     procedure PaintGL(C: TOpenGLControl);
     procedure ColorToGL(AColor: TColor; out R, G, B: GLFloat);
-    function BinToMHz(Bin, BinCount: Integer): Double;
+    function SourceFreqToBin(FreqHz: Double; BinCount: Integer): Integer;
     function FreqToX(FreqHz: Double; W: Integer): Integer;
     function GridStepHz(W: Integer): Double;
     function RulerGridStepHz(C: TCanvas; PlotW: Integer): Double;
@@ -80,6 +89,10 @@ implementation
 
 const
   WB_DB_SCALE_W = 34;
+  WB_WDSP_ID = 32;
+  WB_WDSP_FFT = 16384;
+  WB_WDSP_BLOCK = 512;
+  WB_WDSP_PIXELS = 4096;
 
 constructor TWidebandView.Create;
 begin
@@ -98,12 +111,19 @@ begin
   FSourceStartHz := 0.0;
   FSourceEndHz := FSampleRateHz * 0.5;
   FMarkerHz := 0.0;
+  FWDSPAnalyzerOpen := False;
+  FWDSPAnalyzerConfigured := False;
+  FWDSPAnalyzerID := WB_WDSP_ID;
+  FWDSPFlp[0] := 0;
+  SetLength(FWDSPIn, WB_WDSP_BLOCK * 2);
+  SetLength(FWDSPPixels, WB_WDSP_PIXELS);
   FGLTex := 0;
   FDirty := True;
 end;
 
 destructor TWidebandView.Destroy;
 begin
+  CloseWDSPAnalyzer;
   if FGLTex <> 0 then
     glDeleteTextures(1, @FGLTex);
   FRulerBitmap.Free;
@@ -122,6 +142,7 @@ begin
   if AHz > 0 then
   begin
     FSampleRateHz := AHz;
+    FWDSPAnalyzerConfigured := False;
     SetFrequencyView(0.0, FSampleRateHz * 0.5, 0.0, FSampleRateHz * 0.5);
   end;
 end;
@@ -203,100 +224,109 @@ begin
     FRulerBitmap.SetSize(W, H);
 end;
 
-procedure TWidebandView.ComputeSpectrum(const Samples: array of SmallInt;
-  Count: Integer);
-var
-  N, Half, I, J, K, M, Step: Integer;
-  Wr, Wi, Ur, Ui, Tr, Ti, Ang, Re, Im, Mag, Win, WinSum: Double;
-  RealBuf, ImagBuf: array of Double;
+procedure TWidebandView.CloseWDSPAnalyzer;
 begin
-  if Count < 256 then Exit;
-  N := 1;
-  while (N shl 1 <= Count) and (N shl 1 <= 16384) do
-    N := N shl 1;
-  Half := N div 2;
-  SetLength(RealBuf, N);
-  SetLength(ImagBuf, N);
-  WinSum := 0.0;
-  for I := 0 to N - 1 do
-  begin
-    Win := 0.35875 - 0.48829 * Cos(2 * Pi * I / (N - 1)) +
-           0.14128 * Cos(4 * Pi * I / (N - 1)) -
-           0.01168 * Cos(6 * Pi * I / (N - 1));
-    WinSum := WinSum + Win;
-    RealBuf[I] := Samples[I] / 32768.0 * Win;
-    ImagBuf[I] := 0.0;
-  end;
-  if WinSum <= 0.0 then WinSum := N;
+  if FWDSPAnalyzerOpen and Assigned(@DestroyAnalyzer) then
+    DestroyAnalyzer(FWDSPAnalyzerID);
+  FWDSPAnalyzerOpen := False;
+  FWDSPAnalyzerConfigured := False;
+end;
 
-  J := 0;
-  for I := 1 to N - 2 do
+function TWidebandView.EnsureWDSPAnalyzer: Boolean;
+var
+  Success: Integer;
+  SampleRate: Integer;
+  AvBackmult: Double;
+begin
+  Result := False;
+  if (not Assigned(@XCreateAnalyzer)) or (not Assigned(@SetAnalyzer)) or
+     (not Assigned(@Spectrum0)) or (not Assigned(@GetPixels)) then Exit;
+
+  if not FWDSPAnalyzerOpen then
   begin
-    K := N shr 1;
-    while J >= K do
-    begin
-      Dec(J, K);
-      K := K shr 1;
-    end;
-    Inc(J, K);
-    if I < J then
-    begin
-      Re := RealBuf[I]; RealBuf[I] := RealBuf[J]; RealBuf[J] := Re;
-      Im := ImagBuf[I]; ImagBuf[I] := ImagBuf[J]; ImagBuf[J] := Im;
-    end;
+    Success := -1;
+    XCreateAnalyzer(FWDSPAnalyzerID, @Success, WB_WDSP_FFT, 1, 1, nil);
+    if Success <> 0 then Exit;
+    FWDSPAnalyzerOpen := True;
+    FWDSPAnalyzerConfigured := False;
   end;
 
-  M := 2;
-  while M <= N do
+  if not FWDSPAnalyzerConfigured then
   begin
-    Ang := -2 * Pi / M;
-    Wr := Cos(Ang);
-    Wi := Sin(Ang);
-    Step := M div 2;
-    K := 0;
-    while K < N do
-    begin
-      Ur := 1.0;
-      Ui := 0.0;
-      for J := 0 to Step - 1 do
-      begin
-        I := K + J;
-        Tr := Ur * RealBuf[I + Step] - Ui * ImagBuf[I + Step];
-        Ti := Ur * ImagBuf[I + Step] + Ui * RealBuf[I + Step];
-        RealBuf[I + Step] := RealBuf[I] - Tr;
-        ImagBuf[I + Step] := ImagBuf[I] - Ti;
-        RealBuf[I] := RealBuf[I] + Tr;
-        ImagBuf[I] := ImagBuf[I] + Ti;
-        Re := Ur * Wr - Ui * Wi;
-        Ui := Ur * Wi + Ui * Wr;
-        Ur := Re;
-      end;
-      Inc(K, M);
-    end;
-    M := M shl 1;
+    SampleRate := Max(1, Round(FSampleRateHz));
+    AvBackmult := Exp(-1.0 / (15.0 * 0.120)); // Thetis wideband default.
+    SetAnalyzer(
+      FWDSPAnalyzerID,
+      2,
+      1,
+      1,                 // Feed like piHPSDR: interleaved I/Q via Spectrum0.
+      @FWDSPFlp[0],
+      WB_WDSP_FFT,
+      WB_WDSP_BLOCK,
+      6,                 // Thetis wideband default window.
+      14.0,
+      0,                 // no overlap for discontinuous WB packets.
+      0,
+      0.0,
+      0.0,
+      WB_WDSP_PIXELS,
+      1,
+      0,
+      0.0,
+      0.0,
+      2 * WB_WDSP_FFT);
+    SetDisplayAverageMode(FWDSPAnalyzerID, 0, AVERAGE_MODE_LOG_RECURSIVE);
+    SetDisplayAvBackmult(FWDSPAnalyzerID, 0, AvBackmult);
+    SetDisplaySampleRate(FWDSPAnalyzerID, SampleRate);
+    ResetPixelBuffers(FWDSPAnalyzerID);
+    FWDSPAnalyzerConfigured := True;
   end;
 
-  SetLength(FData, Half);
-  for I := 0 to Half - 1 do
+  Result := True;
+end;
+
+function TWidebandView.ComputeSpectrumWDSP(const Samples: array of SmallInt;
+  Count: Integer): Boolean;
+var
+  Offset, I, Flag: Integer;
+  MaxPix: Single;
+begin
+  Result := False;
+  if Count < WB_WDSP_BLOCK then Exit;
+  if not EnsureWDSPAnalyzer then Exit;
+
+  Offset := 0;
+  while Offset + WB_WDSP_BLOCK <= Count do
   begin
-    Mag := Sqr(RealBuf[I]) + Sqr(ImagBuf[I]);
-    // Convert FFT bin amplitude to dBFS, then apply Thetis-style display
-    // calibration offset. The previous unnormalised FFT power depended on N
-    // and could not line up with the dB scale.
-    FData[I] := 20.0 * Log10((2.0 * Sqrt(Mag) / WinSum) + 1.0E-20) + FCalOffset;
+    for I := 0 to WB_WDSP_BLOCK - 1 do
+    begin
+      FWDSPIn[I * 2] := Samples[Offset + I] / 32768.0;
+      FWDSPIn[I * 2 + 1] := 0.0;
+    end;
+    Spectrum0(1, FWDSPAnalyzerID, 0, 0, @FWDSPIn[0]);
+    Inc(Offset, WB_WDSP_BLOCK);
   end;
+
+  Flag := 0;
+  GetPixels(FWDSPAnalyzerID, 0, @FWDSPPixels[0], @Flag);
+  if Flag = 0 then Exit;
+
+  MaxPix := -1.0E30;
+  for I := 0 to WB_WDSP_PIXELS - 1 do
+    if FWDSPPixels[I] > MaxPix then
+      MaxPix := FWDSPPixels[I];
+  if MaxPix < -250.0 then Exit;
+
+  SetLength(FData, WB_WDSP_PIXELS);
+  for I := 0 to WB_WDSP_PIXELS - 1 do
+    FData[I] := FWDSPPixels[I] + FCalOffset;
+  Result := True;
 end;
 
 procedure TWidebandView.SetSamples(const Samples: array of SmallInt; Count: Integer);
 begin
-  ComputeSpectrum(Samples, Count);
-  FDirty := True;
-end;
-
-function TWidebandView.BinToMHz(Bin, BinCount: Integer): Double;
-begin
-  if BinCount <= 1 then Result := 0
-  else Result := (Bin / (BinCount - 1)) * (FSampleRateHz * 0.5) / 1000000.0;
+  if ComputeSpectrumWDSP(Samples, Count) then
+    FDirty := True;
 end;
 
 function TWidebandView.FreqToX(FreqHz: Double; W: Integer): Integer;
@@ -308,6 +338,19 @@ begin
   Result := Round((FreqHz - FViewStartHz) / SpanHz * Max(1, W - 1));
 end;
 
+function TWidebandView.SourceFreqToBin(FreqHz: Double; BinCount: Integer): Integer;
+var
+  Fs: Double;
+begin
+  if BinCount <= 1 then Exit(0);
+  Fs := FSampleRateHz;
+  if Fs <= 0.0 then Fs := 122880000.0;
+  // We feed WDSP through Spectrum0 as a complex stream with Q=0. GetPixels is
+  // then laid out as -Fs/2..+Fs/2, so raw ADC frequencies 0..Fs/2 live in the
+  // right half of the pixel array. Do not map 0..Nyquist across the whole array.
+  Result := EnsureRange(Round((0.5 + FreqHz / Fs) * (BinCount - 1)), 0, BinCount - 1);
+end;
+
 function TWidebandView.GridStepHz(W: Integer): Double;
 var
   PixPerMHz, SpanMHz: Double;
@@ -536,7 +579,7 @@ begin
         DB := -200.0
       else
       begin
-        Bin := EnsureRange(Round(SrcHz / Nyq * (N - 1)), 0, N - 1);
+        Bin := SourceFreqToBin(SrcHz, N);
         DB := FData[Bin];
       end;
       Y := Round((FRefLevel - DB) * InvRange * TopH);