analytics-api/signal_processing_convolution.ts

// signal-processing.ts - Convolution-based signal processing functions
import { convolve1D, convolve2D, ConvolutionKernels, ConvolutionOptions } from './convolution';

export interface SmoothingOptions {
  method?: 'gaussian' | 'moving_average';
  windowSize?: number;
  sigma?: number;
}

export interface EdgeDetectionOptions {
  method?: 'sobel' | 'laplacian' | 'canny';
  threshold?: number;
}

export interface FilterOptions {
  type: 'lowpass' | 'highpass' | 'bandpass' | 'bandstop';
  cutoffLow?: number;
  cutoffHigh?: number;
  order?: number;
}

export interface DerivativeOptions {
  order?: 1 | 2;
  method?: 'gradient' | 'laplacian';
}

/**
 * Convolution-Based Signal Processing Library
 * Functions that leverage convolution operations for signal processing
 */
export class SignalProcessor {
  
  /**
   * Smooth a 1D signal using convolution-based methods
   */
  static smooth(signal: number[], options: SmoothingOptions = {}): number[] {
    const { method = 'gaussian', windowSize = 5, sigma = 1.0 } = options;

    if (signal.length === 0) {
      throw new Error('Signal cannot be empty');
    }
    
    let kernel: number[];
    
    switch (method) {
      case 'gaussian':
        kernel = ConvolutionKernels.gaussian1D(windowSize, sigma);
        break;
        
      case 'moving_average':
        kernel = ConvolutionKernels.average1D(windowSize);
        break;
           
      default:
        throw new Error(`Unsupported smoothing method: ${method}`);
    }
    
    return convolve1D(signal, kernel, { mode: 'same' }).values;
  }

  /**
   * Detect edges in 2D image data using convolution-based methods
   */
  static detectEdges2D(image: number[][], options: EdgeDetectionOptions = {}): number[][] {
    const { method = 'sobel', threshold = 0.1 } = options;
    
    let kernelX: number[][];
    let kernelY: number[][];
    
    switch (method) {
      case 'sobel':
        kernelX = ConvolutionKernels.sobel("x");
        kernelY = ConvolutionKernels.sobel("y");
        break;
        
      case 'laplacian':
        const laplacianKernel = ConvolutionKernels.laplacian();
        return convolve2D(image, laplacianKernel, { mode: 'same' }).matrix.map(row =>
          row.map(val => Math.abs(val) > threshold ? Math.abs(val) : 0)
        );
        
      default:
        throw new Error(`Unsupported edge detection method: ${method}`);
    }
    
    // Apply both kernels and combine results
    const edgesX = convolve2D(image, kernelX, { mode: 'same' }).matrix;
    const edgesY = convolve2D(image, kernelY, { mode: 'same' }).matrix;
    
    // Calculate gradient magnitude
    const result: number[][] = [];
    for (let i = 0; i < edgesX.length; i++) {
      result[i] = [];
      for (let j = 0; j < edgesX[i].length; j++) {
        const magnitude = Math.sqrt(edgesX[i][j] ** 2 + edgesY[i][j] ** 2);
        result[i][j] = magnitude > threshold ? magnitude : 0;
      }
    }
    
    return result;
  }

  /**
   * Apply digital filters using convolution
   */
  static filter(signal: number[], options: FilterOptions): number[] {
    const { type, cutoffLow = 0.1, cutoffHigh = 0.5, order = 4 } = options;
    
    let kernel: number[];
    
    switch (type) {
      case 'lowpass':
        // Low-pass filter using Gaussian kernel
        kernel = ConvolutionKernels.gaussian1D(order * 4 + 1, order / 2);
        return convolve1D(signal, kernel, { mode: 'same' }).values;
        
      case 'highpass':
        // High-pass filter using difference of Gaussians
        const lpKernel = ConvolutionKernels.gaussian1D(order * 4 + 1, order / 2);
        const smoothed = convolve1D(signal, lpKernel, { mode: 'same' }).values;
        return signal.map((val, i) => val - smoothed[i]);
        
      case 'bandpass':
        // Band-pass as combination of high-pass and low-pass
        const hp = this.filter(signal, { type: 'highpass', cutoffLow, order });
        return this.filter(hp, { type: 'lowpass', cutoffLow: cutoffHigh, order });
        
      case 'bandstop':
        // Band-stop as original minus band-pass
        const bp = this.filter(signal, { type: 'bandpass', cutoffLow, cutoffHigh, order });
        return signal.map((val, i) => val - bp[i]);
        
      default:
        throw new Error(`Unsupported filter type: ${type}`);
    }
  }

  /**
   * Calculate derivatives using convolution with derivative kernels
   */
  static derivative(signal: number[], options: DerivativeOptions = {}): number[] {
    const { order = 1, method = 'gradient' } = options;
    
    let kernel: number[];
    
    if (method === 'gradient') {
      switch (order) {
        case 1:
          // First derivative using gradient kernel
          kernel = [-0.5, 0, 0.5]; // Simple gradient
          break;
        case 2:
          // Second derivative using Laplacian-like kernel
          kernel = [1, -2, 1]; // Simple second derivative
          break;
        default:
          throw new Error(`Unsupported derivative order: ${order}`);
      }
    } else if (method === 'laplacian' && order === 2) {
      // 1D Laplacian
      kernel = [1, -2, 1];
    } else {
      throw new Error(`Unsupported derivative method: ${method}`);
    }
    
    return convolve1D(signal, kernel, { mode: 'same' }).values;
  }

  /**
   * Blur 2D image using Gaussian convolution
   */
  static blur2D(image: number[][], sigma: number = 1.0, kernelSize?: number): number[][] {
    const size = kernelSize || Math.ceil(sigma * 6) | 1; // Ensure odd size
    const kernel = ConvolutionKernels.gaussian(size, sigma);

    return convolve2D(image, kernel, { mode: 'same' }).matrix;
  }

  /**
   * Sharpen 2D image using unsharp masking (convolution-based)
   */
  static sharpen2D(image: number[][], strength: number = 1.0): number[][] {
    const sharpenKernel = [
      [0, -strength, 0],
      [-strength, 1 + 4 * strength, -strength],
      [0, -strength, 0]
    ];
    
    return convolve2D(image, sharpenKernel, { mode: 'same' }).matrix;
  }

  /**
   * Apply emboss effect using convolution
   */
  static emboss2D(image: number[][], direction: 'ne' | 'nw' | 'se' | 'sw' = 'ne'): number[][] {
    const embossKernels = {
      ne: [[-2, -1, 0], [-1, 1, 1], [0, 1, 2]],
      nw: [[0, -1, -2], [1, 1, -1], [2, 1, 0]],
      se: [[0, 1, 2], [-1, 1, 1], [-2, -1, 0]],
      sw: [[2, 1, 0], [1, 1, -1], [0, -1, -2]]
    };
    
    const kernel = embossKernels[direction];
    return convolve2D(image, kernel, { mode: 'same' }).matrix;
  }

  /**
   * Apply motion blur using directional convolution kernel
   */
  static motionBlur(signal: number[], direction: number, length: number = 9): number[] {
    // Create motion blur kernel
    const kernel = new Array(length).fill(1 / length);
    
    return convolve1D(signal, kernel, { mode: 'same' }).values;
  }

  /**
   * Detect impulse response using convolution with known impulse
   */
  static matchedFilter(signal: number[], template: number[]): number[] {
    // Matched filter using cross-correlation (convolution with reversed template)
    const reversedTemplate = [...template].reverse();
    return convolve1D(signal, reversedTemplate, { mode: 'same' }).values;
  }

  /**
   * Apply median filtering (note: not convolution-based, but commonly used with other filters)
   */
  static medianFilter(signal: number[], windowSize: number = 3): number[] {
    const result: number[] = [];
    const halfWindow = Math.floor(windowSize / 2);
    
    for (let i = 0; i < signal.length; i++) {
      const window: number[] = [];
      
      for (let j = Math.max(0, i - halfWindow); j <= Math.min(signal.length - 1, i + halfWindow); j++) {
        window.push(signal[j]);
      }
      
      window.sort((a, b) => a - b);
      const median = window[Math.floor(window.length / 2)];
      result.push(median);
    }
    
    return result;
  }

  /**
   * Cross-correlation using convolution
   */
  static crossCorrelate(signal1: number[], signal2: number[]): number[] {
    // Cross-correlation is convolution with one signal reversed
    const reversedSignal2 = [...signal2].reverse();
    return convolve1D(signal1, reversedSignal2, { mode: 'full' }).values;
  }

  /**
   * Auto-correlation using convolution
   */
  static autoCorrelate(signal: number[]): number[] {
    return this.crossCorrelate(signal, signal);
  }

  /**
   * Detect peaks using convolution-based edge detection
   */
  /**
   * [REWRITTEN] Detects peaks (local maxima) in a 1D signal.
   * This is a more robust method that directly finds local maxima.
   */
  static detectPeaksConvolution(signal: number[], options: {
    smoothWindow?: number;
    threshold?: number;
    minDistance?: number;
  } = {}): { index: number; value: number }[] {
    const { smoothWindow = 0, threshold = -Infinity, minDistance = 1 } = options;

    let processedSignal = signal;
    // Optionally smooth the signal first to reduce noise
    if (smoothWindow > 1) {
      processedSignal = this.smooth(signal, { method: 'gaussian', windowSize: smoothWindow });
    }
    
    const peaks: { index: number; value: number }[] = [];

    // Find all points that are higher than their immediate neighbors
    for (let i = 1; i < processedSignal.length - 1; i++) {
      const prev = processedSignal[i - 1];
      const curr = processedSignal[i];
      const next = processedSignal[i + 1];

      if (curr > prev && curr > next && curr > threshold) {
        peaks.push({ index: i, value: signal[i] }); // Store index and ORIGINAL value
      }
    }

    // Check boundaries: Is the first or last point a peak?
    if (processedSignal[0] > processedSignal[1] && processedSignal[0] > threshold) {
        peaks.unshift({ index: 0, value: signal[0] });
    }
    const last = processedSignal.length - 1;
    if (processedSignal[last] > processedSignal[last - 1] && processedSignal[last] > threshold) {
        peaks.push({ index: last, value: signal[last] });
    }
    
    // [CORRECTED LOGIC] Enforce minimum distance between peaks
    if (minDistance < 2 || peaks.length <= 1) {
        return peaks;
    }
    
    // Sort peaks by value, highest first
    peaks.sort((a, b) => b.value - a.value);

    const finalPeaks: { index: number; value: number }[] = [];
    const removed = new Array(peaks.length).fill(false);

    for (let i = 0; i < peaks.length; i++) {
        if (!removed[i]) {
            finalPeaks.push(peaks[i]);
            // Remove other peaks within the minimum distance
            for (let j = i + 1; j < peaks.length; j++) {
                if (!removed[j] && Math.abs(peaks[i].index - peaks[j].index) < minDistance) {
                    removed[j] = true;
                }
            }
        }
    }
    
    return finalPeaks.sort((a, b) => a.index - b.index);
  }


  /**
   * [REWRITTEN] Detects valleys (local minima) in a 1D signal.
   */
  static detectValleysConvolution(signal: number[], options: {
    smoothWindow?: number;
    threshold?: number;
    minDistance?: number;
  } = {}): { index: number; value: number }[] {
    const invertedSignal = signal.map(x => -x);
    const invertedThreshold = options.threshold !== undefined ? -options.threshold : undefined;

    const invertedPeaks = this.detectPeaksConvolution(invertedSignal, { ...options, threshold: invertedThreshold });
    
    return invertedPeaks.map(peak => ({
      index: peak.index,
      value: -peak.value,
    }));
  }

  /**
   * Detect outliers using convolution-based methods
   */
  /**
   * [REWRITTEN] Detects outliers using more reliable and statistically sound methods.
   */
  static detectOutliersConvolution(signal: number[], options: {
    method?: 'local_deviation' | 'mean_diff';
    windowSize?: number;
    threshold?: number;
  } = {}): { index: number; value: number; outlierScore: number }[] {
    const { method = 'local_deviation', windowSize = 7, threshold = 3.0 } = options;
    
    let outlierScores: number[];
    
    switch (method) {
      case 'mean_diff':
        // Detects outliers by their difference from the local mean.
        const meanKernel = ConvolutionKernels.average1D(windowSize);
        const localMean = convolve1D(signal, meanKernel, { mode: 'same' }).values;
        outlierScores = signal.map((val, i) => Math.abs(val - localMean[i]));
        break;
        
      case 'local_deviation':
        // A robust method using Z-score: how many local standard deviations away a point is.
        const avgKernel = ConvolutionKernels.average1D(windowSize);
        const localMeanValues = convolve1D(signal, avgKernel, { mode: 'same' }).values;
        const squaredDiffs = signal.map((val, i) => (val - localMeanValues[i]) ** 2);
        const localVar = convolve1D(squaredDiffs, avgKernel, { mode: 'same' }).values;
        outlierScores = signal.map((val, i) => {
          const std = Math.sqrt(localVar[i]);
          return std > 1e-6 ? Math.abs(val - localMeanValues[i]) / std : 0;
        });
        break;
        
      default:
        throw new Error(`Unsupported outlier detection method: ${method}`);
    }
    
    // Find points exceeding the threshold
    const outliers: { index: number; value: number; outlierScore: number }[] = [];
    outlierScores.forEach((score, i) => {
      if (score > threshold) {
        outliers.push({
          index: i,
          value: signal[i],
          outlierScore: score
        });
      }
    });
    
    return outliers;
  }

  /**
   * Detect trend vertices (turning points) using convolution
   */
  /**
   * [CORRECTED] Detects trend vertices (turning points) by finding all peaks and valleys.
   * This version fixes a bug that prevented valleys from being detected.
   */
  static detectTrendVertices(signal: number[], options: {
    smoothingWindow?: number;
    threshold?: number;
    minDistance?: number;
  } = {}): { index: number; value: number; type: 'peak' | 'valley' }[] {
    const { 
      smoothingWindow = 5, 
      threshold = 0, // CORRECTED: Changed default from -Infinity to a sensible 0
      minDistance = 3
    } = options;
    
    // Create the options object to pass down. The valley function will handle inverting the threshold itself.
    const detectionOptions = { smoothingWindow, threshold, minDistance };

    const peaks = this.detectPeaksConvolution(signal, detectionOptions).map(p => ({ ...p, type: 'peak' as const }));
    const valleys = this.detectValleysConvolution(signal, detectionOptions).map(v => ({ ...v, type: 'valley' as const }));
    
    // Combine peaks and valleys and sort them by their index to get the sequence of trend changes
    const vertices = [...peaks, ...valleys];
    vertices.sort((a, b) => a.index - b.index);
    
    return vertices;
  }

  /**
   * Detect vertices using curvature (second derivative)
   */
  private static detectCurvatureVertices(
    signal: number[], 
    threshold: number
  ): { index: number; value: number; type: 'peak' | 'valley'; curvature: number }[] {
    // Use second derivative kernel for curvature
    const curvatureKernel = [1, -2, 1]; // Discrete Laplacian
    const curvature = convolve1D(signal, curvatureKernel, { mode: 'same' }).values;
    
    const vertices: { index: number; value: number; type: 'peak' | 'valley'; curvature: number }[] = [];
    
    // Find zero crossings in curvature with sufficient magnitude
    for (let i = 1; i < curvature.length - 1; i++) {
      const prev = curvature[i - 1];
      const curr = curvature[i];
      const next = curvature[i + 1];
      
      // Zero crossing detection
      if ((prev > 0 && next < 0) || (prev < 0 && next > 0)) {
        const curvatureMagnitude = Math.abs(curr);
        
        if (curvatureMagnitude > threshold) {
          const type: 'peak' | 'valley' = prev > 0 ? 'peak' : 'valley';
          
          vertices.push({
            index: i,
            value: signal[i],
            type,
            curvature: curr
          });
        }
      }
    }
    
    return vertices;
  }

  /**
   * Detect vertices using gradient sign changes
   */
  private static detectSignChangeVertices(
    signal: number[], 
    threshold: number
  ): { index: number; value: number; type: 'peak' | 'valley'; curvature: number }[] {
    // First derivative for gradient
    const gradientKernel = [-0.5, 0, 0.5]; // Central difference
    const gradient = convolve1D(signal, gradientKernel, { mode: 'same' }).values;
    
    // Second derivative for curvature
    const curvatureKernel = [1, -2, 1];
    const curvature = convolve1D(signal, curvatureKernel, { mode: 'same' }).values;
    
    const vertices: { index: number; value: number; type: 'peak' | 'valley'; curvature: number }[] = [];
    
    // Find gradient sign changes
    for (let i = 1; i < gradient.length - 1; i++) {
      const prevGrad = gradient[i - 1];
      const nextGrad = gradient[i + 1];
      
      // Check for sign change with sufficient gradient magnitude
      if (Math.abs(prevGrad) > threshold && Math.abs(nextGrad) > threshold) {
        if ((prevGrad > 0 && nextGrad < 0) || (prevGrad < 0 && nextGrad > 0)) {
          const type: 'peak' | 'valley' = prevGrad > 0 ? 'peak' : 'valley';
          
          vertices.push({
            index: i,
            value: signal[i],
            type,
            curvature: curvature[i]
          });
        }
      }
    }
    
    return vertices;
  }

  /**
   * Detect vertices using momentum changes
   */
  private static detectMomentumVertices(
    signal: number[], 
    threshold: number
  ): { index: number; value: number; type: 'peak' | 'valley'; curvature: number }[] {
    // Create momentum kernel (difference over larger window)
    const momentumKernel = [-1, 0, 0, 0, 1]; // 4-point difference
    const momentum = convolve1D(signal, momentumKernel, { mode: 'same' }).values;
    
    // Detect momentum reversals
    const momentumGradient = convolve1D(momentum, [-0.5, 0, 0.5], { mode: 'same' }).values;
    const curvature = convolve1D(signal, [1, -2, 1], { mode: 'same' }).values;
    
    const vertices: { index: number; value: number; type: 'peak' | 'valley'; curvature: number }[] = [];
    
    for (let i = 2; i < momentum.length - 2; i++) {
      const prevMomentum = momentum[i - 1];
      const currMomentum = momentum[i];
      const nextMomentum = momentum[i + 1];
      
      // Check for momentum reversal
      if (Math.abs(momentumGradient[i]) > threshold) {
        if ((prevMomentum > 0 && nextMomentum < 0) || (prevMomentum < 0 && nextMomentum > 0)) {
          const type: 'peak' | 'valley' = prevMomentum > 0 ? 'peak' : 'valley';
          
          vertices.push({
            index: i,
            value: signal[i],
            type,
            curvature: curvature[i]
          });
        }
      }
    }
    
    return vertices;
  }

  /**
   * Detect trend direction changes using convolution
   */
  static detectTrendChanges(signal: number[], options: {
    windowSize?: number;
    threshold?: number;
    minTrendLength?: number;
  } = {}): { index: number; fromTrend: 'up' | 'down' | 'flat'; toTrend: 'up' | 'down' | 'flat'; strength: number }[] {
    const { windowSize = 10, threshold = 0.01, minTrendLength = 5 } = options;
    
    // Calculate local trends using convolution with trend-detecting kernel
    const trendKernel = new Array(windowSize).fill(0).map((_, i) => {
      const center = windowSize / 2;
      return (i - center) / (windowSize * windowSize / 12); // Linear trend kernel
    });
    
    const trends = convolve1D(signal, trendKernel, { mode: 'same' }).values;
    
    // Classify trends
    const trendDirection = trends.map(t => {
      if (t > threshold) return 'up';
      if (t < -threshold) return 'down';
      return 'flat';
    });
    
    // Find trend changes
    const changes: { index: number; fromTrend: 'up' | 'down' | 'flat'; toTrend: 'up' | 'down' | 'flat'; strength: number }[] = [];
    
    let currentTrend: 'up' | 'down' | 'flat' = trendDirection[0];
    let trendStart = 0;
    
    for (let i = 1; i < trendDirection.length; i++) {
      if (trendDirection[i] !== currentTrend) {
        const trendLength = i - trendStart;
        
        if (trendLength >= minTrendLength) {
          changes.push({
            index: i,
            fromTrend: currentTrend,
            toTrend: trendDirection[i] as 'up' | 'down' | 'flat',
            strength: Math.abs(trends[i] - trends[trendStart])
          });
        }
        
        currentTrend = trendDirection[i] as 'up' | 'down' | 'flat';
        trendStart = i;
      }
    }
    
    return changes;
  }

  /**
   * Enforce minimum distance between vertices
   */
  private static enforceMinDistanceVertices(
    vertices: { index: number; value: number; type: 'peak' | 'valley'; curvature: number }[],
    minDistance: number
  ): { index: number; value: number; type: 'peak' | 'valley'; curvature: number }[] {
    if (vertices.length <= 1) return vertices;
    
    // Sort by curvature magnitude (stronger vertices first)
    const sorted = [...vertices].sort((a, b) => Math.abs(b.curvature) - Math.abs(a.curvature));
    const result: { index: number; value: number; type: 'peak' | 'valley'; curvature: number }[] = [];
    
    for (const vertex of sorted) {
      let tooClose = false;
      
      for (const accepted of result) {
        if (Math.abs(vertex.index - accepted.index) < minDistance) {
          tooClose = true;
          break;
        }
      }
      
      if (!tooClose) {
        result.push(vertex);
      }
    }
    
    // Sort result by index
    return result.sort((a, b) => a.index - b.index);
  }

  /**
   * Enforce minimum distance between detected features
   */
  private static enforceMinDistanceConv(
    features: { index: number; value: number; strength: number }[], 
    minDistance: number
  ): { index: number; value: number; strength: number }[] {
    if (features.length <= 1) return features;
    
    // Sort by strength (descending)
    const sorted = [...features].sort((a, b) => b.strength - a.strength);
    const result: { index: number; value: number; strength: number }[] = [];
    
    for (const feature of sorted) {
      let tooClose = false;
      
      for (const accepted of result) {
        if (Math.abs(feature.index - accepted.index) < minDistance) {
          tooClose = true;
          break;
        }
      }
      
      if (!tooClose) {
        result.push(feature);
      }
    }
    
    // Sort result by index
    return result.sort((a, b) => a.index - b.index);
  }
}