Overmodulation

// WHAT: Simulate visible picture and colour damage from analogue-TV
// overmodulation.
// HOW: Encode each RGB pixel into eight representative composite-video
// samples, pass every sample through an idealized RF-envelope fold or limit,
// approximately decode NTSC/PAL/SECAM colour, and convert the result to RGB.
// WHY: Luminance and encoded chrominance share one composite excursion, so
// applying abs() or clipping directly to RGB would miss their interaction.

// OVERMODULATION — VISUAL SIMULATION
// ----------------------------------
// This is not a calibrated television transmitter or receiver.
// It keeps the important order of events:
//
//   RGB picture
//      -> luminance + two colour-difference signals
//      -> one composite waveform
//      -> RF-envelope fold or limit
//      -> approximate colour decoder
//      -> displayed RGB
//
// We do not generate millions of RF cycles. For each picture pixel, we build
// only eight representative samples of the composite waveform.


// 1. CONSTANTS
// ------------

// Eight samples represent one turn around the virtual colour subcarrier.
const int QAM_SAMPLE_COUNT = 8;

// A convenient visual boundary for the composite signal. This is not a voltage
// from PAL, NTSC, or SECAM. It keeps ordinary colours in range at 100% drive.
const float COMPOSITE_LIMIT = 1.4;

const float TAU = 6.28318530718;       // 2*pi: one complete cycle in radians
const float SQRT_HALF = 0.70710678;    // sqrt(1/2), also cos(45 degrees)


// 2. RGB <-> VIDEO COMPONENTS
// ---------------------------
// GLSL's vec3 stores three numbers. Here they mean:
//
//   video.x = Y: luminance, or approximate black-to-white information
//   video.y = U: blue minus luminance, scaled for video
//   video.z = V: red minus luminance, scaled for video
//
// Analogue colour television combines Y, U, and V before transmission. That is
// why applying abs() directly to R, G, B, or Y would model the wrong operation.

vec3 rgbToVideoComponents(vec3 rgb) {
  float luminance = dot(rgb, vec3(0.299, 0.587, 0.114));
  float blueDifference = 0.492 * (rgb.b - luminance);
  float redDifference = 0.877 * (rgb.r - luminance);
  return vec3(luminance, blueDifference, redDifference);
}

vec3 videoComponentsToRgb(vec3 video) {
  float red = video.x + video.z / 0.877;
  float blue = video.x + video.y / 0.492;
  float green = (video.x - 0.299 * red - 0.114 * blue) / 0.587;
  return vec3(red, green, blue);
}


// 3. ONE IDEALIZED RF NONLINEARITY
// --------------------------------

// qamBasis() returns (cos(phi), sin(phi)) for one of eight angles:
// 0, 45, 90, ... 315 degrees. These are the momentary weights of U and V in:
//
//   composite = Y + U*cos(phi) + V*sin(phi)
//
// A lookup table is cheaper and easier to compare than calling cos() and sin()
// for every NTSC/PAL sample.

vec2 qamBasis(int sampleIndex) {
  if (sampleIndex == 0) return vec2(1.0, 0.0);              //   0 degrees
  if (sampleIndex == 1) return vec2(SQRT_HALF, SQRT_HALF);  //  45 degrees
  if (sampleIndex == 2) return vec2(0.0, 1.0);              //  90 degrees
  if (sampleIndex == 3) return vec2(-SQRT_HALF, SQRT_HALF); // 135 degrees
  if (sampleIndex == 4) return vec2(-1.0, 0.0);             // 180 degrees
  if (sampleIndex == 5) return vec2(-SQRT_HALF, -SQRT_HALF);// 225 degrees
  if (sampleIndex == 6) return vec2(0.0, -1.0);             // 270 degrees
  return vec2(SQRT_HALF, -SQRT_HALF);                       // 315 degrees
}

// Pass one composite sample through the transmitter/receiver approximation.
// normalizedComposite=1 means that this sample has reached the chosen boundary.
// depth=1 means 100% requested modulation; depth=1.4 means 140%.
// visionPolarity=0 selects common negative vision AM; 1 selects positive AM.

float recoverCompositeSample(
  float normalizedComposite,
  float depth,
  float visionPolarity
) {
  float requestedExcursion = depth * normalizedComposite;

  // NEGATIVE VISION AM
  // The signed carrier scale is 1-requestedExcursion. Above the boundary it
  // crosses zero and reverses RF phase. An ideal envelope detector forgets the
  // sign, hence abs(). Decoding that magnitude folds the sample backwards.
  float signedCarrierScale = 1.0 - requestedExcursion;
  float detectedEnvelope = abs(signedCarrierScale);
  float negativeAmRecovered = (1.0 - detectedEnvelope) / depth;

  // POSITIVE VISION AM
  // There is no zero-carrier fold on the white side. This visual model instead
  // stops the requested excursion at a normalized upper transmitter limit.
  float positiveAmRecovered = min(requestedExcursion, 1.0) / depth;

  // step() returns 0 for negative AM and 1 for positive AM. mix() then chooses
  // the corresponding result without a dynamic if statement.
  float usePositiveAm = step(0.5, visionPolarity);
  return mix(negativeAmRecovered, positiveAmRecovered, usePositiveAm);
}


// 4A. NTSC / PAL: QAM COLOUR
// --------------------------
// NTSC and PAL place two colour-difference signals on cosine and sine versions
// of the same colour subcarrier. PAL additionally reverses V on alternate lines.

vec3 decodeQamProxy(
  vec3 originalVideo,
  float depth,
  float colourSystem,
  float visionPolarity,
  float sourceLine
) {
  // colourSystem=0 is NTSC; colourSystem=1 is PAL.
  bool isPal = colourSystem > 0.5;
  bool isOddSourceLine = mod(floor(sourceLine), 2.0) > 0.5;
  float palLineSign = (isPal && isOddSourceLine) ? -1.0 : 1.0;

  // PAL changes the sign before transmission and changes it back after decoding.
  float encodedRedDifference = originalVideo.z * palLineSign;

  // These sums are the simple decoder. The mean will recover luminance.
  // Correlation with cosine and sine will recover the two colour differences.
  float luminanceSum = 0.0;
  float blueDifferenceSum = 0.0;
  float redDifferenceSum = 0.0;

  for (int sampleIndex = 0; sampleIndex < QAM_SAMPLE_COUNT; sampleIndex++) {
    vec2 phaseWeights = qamBasis(sampleIndex);

    // ENCODE: all three components now share one composite sample.
    float compositeSample = originalVideo.x
      + originalVideo.y * phaseWeights.x
      + encodedRedDifference * phaseWeights.y;

    // TRANSMIT + RECEIVE: overload acts on the combined sample, not separately
    // on luminance and colour. This creates their visible interaction.
    float normalizedSample = compositeSample / COMPOSITE_LIMIT;
    float recoveredSample = COMPOSITE_LIMIT * recoverCompositeSample(
      normalizedSample,
      depth,
      visionPolarity
    );

    // DECODE: collect the DC, cosine, and sine parts again.
    luminanceSum += recoveredSample;
    blueDifferenceSum += recoveredSample * phaseWeights.x;
    redDifferenceSum += recoveredSample * phaseWeights.y;
  }

  // Eight-sample mean gives Y. The factor 2/N = 1/4 gives each QAM component.
  float recoveredLuminance = luminanceSum / 8.0;
  float recoveredBlueDifference = blueDifferenceSum / 4.0;
  float recoveredRedDifference = palLineSign * redDifferenceSum / 4.0;

  return vec3(
    recoveredLuminance,
    recoveredBlueDifference,
    recoveredRedDifference
  );
}


// 4B. SECAM: SIMPLIFIED FM COLOUR
// --------------------------------
// Real SECAM sends one frequency-modulated colour difference per line and uses
// a one-line delay to reconstruct the other. Implementing a complete FM decoder
// would be expensive and misleading here, so this branch is explicitly a proxy.

vec3 decodeSecamProxy(
  vec3 originalVideo,
  float depth,
  float visionPolarity,
  float sourceLine
) {
  bool transmitRedDifference = mod(floor(sourceLine), 2.0) > 0.5;
  float transmittedDifference = transmitRedDifference
    ? originalVideo.z
    : originalVideo.y;

  // Colour changes the virtual tone frequency, not its amplitude. The numbers
  // below are visual normalization choices, not SECAM frequencies in MHz.
  float virtualCycles = 1.25
    + 0.38 * clamp(transmittedDifference, -0.65, 0.65);
  float virtualChromaAmplitude = 0.34;

  // FIRST PASS: measure the average level before and after overload.
  // Fractional cycles do not average to exactly zero in eight samples, so we
  // calculate both means explicitly instead of accidentally changing luma.
  float cleanSampleSum = 0.0;
  float recoveredSampleSum = 0.0;

  for (int sampleIndex = 0; sampleIndex < QAM_SAMPLE_COUNT; sampleIndex++) {
    float phase = TAU * virtualCycles * (float(sampleIndex) + 0.5) / 8.0;
    float fmTone = cos(phase);
    float compositeSample = originalVideo.x + virtualChromaAmplitude * fmTone;
    float recoveredSample = COMPOSITE_LIMIT * recoverCompositeSample(
      compositeSample / COMPOSITE_LIMIT,
      depth,
      visionPolarity
    );

    cleanSampleSum += compositeSample;
    recoveredSampleSum += recoveredSample;
  }

  float cleanMean = cleanSampleSum / 8.0;
  float recoveredMean = recoveredSampleSum / 8.0;

  // SECOND PASS: ask how much of the expected FM-like tone remains. Correlation
  // means multiplying by the known cosine/sine shape and adding the results.
  float cleanCosineCorrelation = 0.0;
  float cleanSineCorrelation = 0.0;
  float recoveredCosineCorrelation = 0.0;
  float recoveredSineCorrelation = 0.0;

  for (int sampleIndex = 0; sampleIndex < QAM_SAMPLE_COUNT; sampleIndex++) {
    float phase = TAU * virtualCycles * (float(sampleIndex) + 0.5) / 8.0;
    float cosineReference = cos(phase);
    float sineReference = sin(phase);
    float compositeSample = originalVideo.x
      + virtualChromaAmplitude * cosineReference;
    float recoveredSample = COMPOSITE_LIMIT * recoverCompositeSample(
      compositeSample / COMPOSITE_LIMIT,
      depth,
      visionPolarity
    );

    cleanCosineCorrelation += (compositeSample - cleanMean) * cosineReference;
    cleanSineCorrelation += (compositeSample - cleanMean) * sineReference;
    recoveredCosineCorrelation +=
      (recoveredSample - recoveredMean) * cosineReference;
    recoveredSineCorrelation +=
      (recoveredSample - recoveredMean) * sineReference;
  }

  float cleanToneStrength = max(length(vec2(
    cleanCosineCorrelation,
    cleanSineCorrelation
  )), 0.0001);
  float recoveredToneStrength = length(vec2(
    recoveredCosineCorrelation,
    recoveredSineCorrelation
  ));

  // 1 means the expected tone survived; 0 means it became unusable. This is a
  // visual stand-in for limiter/discriminator confidence, not measured hardware.
  float colourConfidence = clamp(
    recoveredToneStrength / cleanToneStrength,
    0.0,
    1.0
  );

  vec3 recoveredVideo = originalVideo;

  // Preserve clean luma below overload. Apply only the change introduced by
  // the nonlinear transfer, not the non-zero mean of our short virtual tone.
  recoveredVideo.x += recoveredMean - cleanMean;

  // Only one colour difference is transmitted on this line. The untouched one
  // represents the value that a real receiver would obtain from its delay line.
  if (transmitRedDifference) {
    recoveredVideo.z *= colourConfidence;
  } else {
    recoveredVideo.y *= colourConfidence;
  }

  return recoveredVideo;
}


// 5. PUBLIC SHADER FUNCTION
// -------------------------
// colourSystem:   0=NTSC, 1=PAL, 2=SECAM proxy
// visionPolarity: 0=negative AM, 1=positive AM
// sourcePixel.y:  original image-row number, needed for PAL/SECAM alternation

vec3 applyOvermodulation(
  vec3 sourceRgb,
  float modulationDepth,
  float colourSystem,
  float visionPolarity,
  vec2 sourcePixel
) {
  vec3 originalVideo = rgbToVideoComponents(sourceRgb);
  float safeDepth = max(modulationDepth, 0.01); // avoid division by zero

  bool useQamColour = colourSystem < 1.5;
  vec3 recoveredVideo = useQamColour
    ? decodeQamProxy(
        originalVideo,
        safeDepth,
        colourSystem,
        visionPolarity,
        sourcePixel.y
      )
    : decodeSecamProxy(
        originalVideo,
        safeDepth,
        visionPolarity,
        sourcePixel.y
      );

  // Convert to display RGB only after the receiver approximation. Clamping any
  // earlier would hide the interaction between luminance and encoded colour.
  vec3 recoveredRgb = videoComponentsToRgb(recoveredVideo);
  return clamp(recoveredRgb, 0.0, 1.0);
}