Moonlight Sonata Waterfall

audiofft.py

#!python
# Spectrogram preprocessor
# 
#
# Python code and Mathematics by Waldeck Schutzer (@infinitymathart)
# See the license information at the end of this file.
#
# This app uses the Fast Fourier Transform to separate the frequencies
# present in an audio file. It creates a binary file with each individual
# frequency, typically resampled to 60Hz to drive an animation (other
# sample rates are possible). The file format is easily inferred from
# the code.
#
import numpy as np
import argparse
import matplotlib.pyplot as plt
from scipy.signal import resample, get_window, detrend
from scipy.fft import rfft, rfftfreq
from scipy.io import wavfile
import struct

def old_compute_spectrogram(data, sample_rate, NFFT, noverlap, window_func, detrend_func):
    hop_size = NFFT - noverlap
    num_frames = (len(data) - NFFT) // hop_size + 1
    window = get_window(window_func, NFFT)
    spectrogram = np.zeros((NFFT // 2 + 1, num_frames), dtype=np.float32)
    print("Audio data is",len(data),"long")
    print("Have",num_frames,"to process")

    for i in range(num_frames):
        start_index = i * hop_size
        end_index = start_index + NFFT
        segment = data[start_index:end_index]

        if detrend_func:
            segment = detrend_func(segment)

        segment = segment * window
        fft_result = rfft(segment, n=NFFT)
        power_density = np.abs(fft_result) ** 2
        window_power = np.sum(window ** 2)
        power_density /= (window_power * sample_rate)

        spectrogram[:, i] = power_density

    times = np.arange(num_frames) * hop_size / sample_rate
    freqs = rfftfreq(NFFT, 1 / sample_rate)
    print("Computed",spectrogram.shape[1],"frames")

    return spectrogram, freqs, times

def compute_spectrogram(data, sample_rate, NFFT, target_fps, window_func, detrend_func):
    # Calculate hop size based on the desired output frame rate
    hop_size = int(sample_rate / target_fps)
    noverlap = NFFT - hop_size
    
    # Calculate the total number of frames based on hop size
    num_frames = (len(data) - NFFT) // hop_size + 1
    window = get_window(window_func, NFFT)
    spectrogram = np.zeros((NFFT // 2 + 1, num_frames), dtype=np.float32)

    print("Audio data length:", len(data))
    print("Window size (NFFT):", NFFT)
    print("Hop size:", hop_size)
    print("Number of frames to process:", num_frames)

    for i in range(num_frames):
        start_index = i * hop_size
        end_index = start_index + NFFT
        segment = data[start_index:end_index]

        if detrend_func:
            segment = detrend_func(segment)

        segment = segment * window
        fft_result = rfft(segment)
        power_density = np.abs(fft_result) ** 2
        window_power = np.sum(window ** 2)
        power_density /= (window_power * sample_rate)

        spectrogram[:, i] = power_density

    times = np.arange(num_frames) * hop_size / sample_rate
    freqs = rfftfreq(NFFT, 1 / sample_rate)
    print("Computed", spectrogram.shape[1], "frames")

    return spectrogram, freqs, times

def resample_spectrogram(spectrogram, original_times, target_rate):
    num_target_frames = int(original_times[-1] * target_rate)
    resampled_spectrogram = np.zeros((spectrogram.shape[0], num_target_frames), dtype=np.float32)

    # Resample each frequency bin independently
    for i in range(spectrogram.shape[0]):
        resampled_spectrogram[i] = resample(spectrogram[i], num_target_frames)

    return resampled_spectrogram

def filter_frequencies(spectrogram, freqs, lower_freq, higher_freq):
    # Find the indices of the frequencies within the specified range
    freq_indices = np.where((freqs >= lower_freq) & (freqs <= higher_freq))[0]
    
    # Filter the spectrogram and frequencies based on these indices
    filtered_spectrogram = spectrogram[freq_indices, :]
    filtered_freqs = freqs[freq_indices]

    return filtered_spectrogram, filtered_freqs

def compute_min_max(spectrogram):
    min_values = np.min(spectrogram, axis=1)
    max_values = np.max(spectrogram, axis=1)
    return min_values, max_values

def clean_values(arr):
    # Replace any NaN, -Inf, or +Inf values with 0.0
    arr = np.nan_to_num(arr, nan=0.0, posinf=0.0, neginf=0.0)
    return arr

def write_spectrogram_to_file(output_file, sample_rate, freqs, resampled_spectrogram):
    with open(output_file, 'wb') as f:
        # Prepare header information
        num_frequencies = len(freqs)
        total_packets = resampled_spectrogram.shape[1]

        # Create the header packet
        header = struct.pack('<i', sample_rate)                 # Sample Rate as 32-bit integer
        header += struct.pack('<i', num_frequencies)            # Number of frequencies as 32-bit integer
        header += struct.pack('<i', total_packets)              # Total packets as 32-bit integer

        # Pad the header with zeros to reach 128 bytes
        header = header.ljust(128, b'\x00')
        f.write(header)

        # Write the frequency values (first packet)
        f.write(freqs.astype(np.float32).tobytes())

        # Clean the spectrogram values before computing min and max
        cleaned_spectrogram = clean_values(resampled_spectrogram)

        # Compute min and max values for each frequency bin
        min_values = np.min(cleaned_spectrogram, axis=1)
        max_values = np.max(cleaned_spectrogram, axis=1)

        # Write the min and max values
        f.write(min_values.astype(np.float32).tobytes())
        f.write(max_values.astype(np.float32).tobytes())

        # Write the spectrum values (subsequent packets)
        for frame in cleaned_spectrogram.T:
            f.write(frame.astype(np.float32).tobytes())
        
        print("\n\nHeader and", total_packets, "packets written.")
        print("Sample rate:", sample_rate)
        print("Number of frequencies:", num_frequencies)

def plot_spectrogram(spectrogram, freqs, times):
    # Convert the power spectrum to dB scale
    dB_spectrogram = 10 * np.log10(spectrogram + 1e-10)
    
    # Replace any -inf or NaN values with -100 dB
    dB_spectrogram = np.where(np.isfinite(dB_spectrogram), dB_spectrogram, -100)

    plt.figure(figsize=(10, 6))
    plt.imshow(dB_spectrogram, aspect='auto', origin='lower',
               extent=[times.min(), times.max(), freqs.min(), freqs.max()])
    plt.colorbar(label='Power Spectral Density (dB)')
    plt.xlabel('Time (s)')
    plt.ylabel('Frequency (Hz)')
    plt.title('Spectrogram')
    plt.show()

if __name__ == "__main__":
    parser = argparse.ArgumentParser(description="Compute and save a spectrogram resampled to 60Hz.")
    parser.add_argument("input_file", type=str, help="Path to the input WAV file")
    parser.add_argument("output_file", type=str, help="Path to the output binary file")
    parser.add_argument("-n", "--nfft", type=int, default=1024, help="Number of data points per segment (default: 1024)")
    parser.add_argument("-r", "--frame_rate", type=float, default=60.0, help="Target output frame rate in FPS (default: 60)")
    parser.add_argument("-w", "--window", type=str, default='hann', help="Window function to use (default: 'hann')")
    parser.add_argument("-d", "--detrend", type=str, choices=['none', 'mean', 'linear'], default='none', help="Detrending method (default: 'none')")
    parser.add_argument("-lf", "--lower_frequency", type=float, default=0.0, help="Lower frequency limit for filtering (default: 0.0 Hz)")
    parser.add_argument("-hf", "--higher_frequency", type=float, default=float('inf'), help="Higher frequency limit for filtering (default: no upper limit)")

    args = parser.parse_args()

    # Read the WAV file
    sample_rate, audio_data = wavfile.read(args.input_file)
    if len(audio_data.shape) == 2:  # Stereo
        audio_data = np.mean(audio_data, axis=1)  # Convert to mono by averaging channels

    if args.detrend == 'none':
        detrend_func = None
    elif args.detrend == 'mean':
        detrend_func = lambda x: detrend(x, type='constant')
    elif args.detrend == 'linear':
        detrend_func = lambda x: detrend(x, type='linear')

    spectrogram, freqs, times = compute_spectrogram(audio_data, sample_rate, args.nfft, args.frame_rate, args.window, detrend_func)

    # Filter the spectrogram to keep only frequencies within the specified range
    filtered_spectrogram, filtered_freqs = filter_frequencies(spectrogram, freqs, args.lower_frequency, args.higher_frequency)

    # Write the resampled and filtered spectrogram to a file
    write_spectrogram_to_file(args.output_file, sample_rate, filtered_freqs, filtered_spectrogram)

    # Optional: plot the spectrogram
    plot_spectrogram(filtered_spectrogram, filtered_freqs, np.linspace(0, times[-1], filtered_spectrogram.shape[1]))

# License:
#
# Copyright (c) 2024 Waldeck Schutzer
#
# All rights reserved.
#
# This code after the template and the related animations are the property of the
# copyright holder. Any reproduction, distribution, or use of this material,
# in whole or in part, without the express written permission of the copyright
# holder is strictly prohibited.

sketch_241026a_moonlight.pde

// Sound visualizer animation
//
// Processing code and Mathematics by Waldeck Schutzer (@infinitymathart)
// See the license information at the end of this file.
//
// The sound is first processed with a separate python app available here:
// which creates a binary spectrogram file with a sample rate of 60Hz with
// each individual frequency, thanks to the Fast Fourier Transform.
// This binary data is read into this app by the SpectrogramReader class.
//
import peasy.*;
PeasyCam cam;
boolean recording = false;
int frameWidth = recording ? 1080 : 600;
int frameHeight = recording ? 1920 : 1067;
SpectrogramReader reader;
int numFrames = 0; // Total number of frames in the spectrogram
int currentFrame = 0; // Current frame index for display
PImage spectrogramImage; // Image to display the spectrogram
boolean reading = true;
boolean stamping = true;
track[] tracks;
float fac = 1.0;
float track_dist = 20.0;
float hook_size = 0.1;
final int padding = 380;
int padding_count = 0;
float h_stretch = 1.4;
final int max_samples = int(h_stretch*600);
final float filterq = 0.7;
double h = 4*frameWidth; // height of waterfall
double d = h_stretch*frameWidth; // horizontal distance until water hits the ground
double d2 = d*d;
double alpha = Math.sqrt(2*h)/d;
double alpha2 = 2*h/d2;
float t = 0;

final float[] track_freqs = { // Notes of a Grand Piano
   27.5000, 29.1352, 30.8677, // A0, A#0, B0
            32.7032, 34.6478, 36.7081, 38.8909, 41.2034, 43.6535, 46.2493, 48.9994, 51.9131, // C1 to B1
            55.0000, 58.2705, 61.7354, 65.4064, 69.2957, 73.4162, 77.7817, 82.4069, 87.3071, 92.4986, 97.9989, 103.826, // A2 to B2
            110.000, 116.541, 123.471, 130.813, 138.591, 146.832, 155.563, 164.814, 174.614, 184.997, 195.998, 207.652, // A3 to B3
            220.000, 233.082, 246.942, 261.626, 277.183, 293.665, 311.127, 329.628, 349.228, 369.994, 391.995, 415.305, // A4 to B4
            440.000, 466.164, 493.883, 523.251, 554.365, 587.330, 622.254, 659.255, 698.456, 739.989, 783.991, 830.609, // A5 to B5
            880.000, 932.328, 987.767, 1046.50, 1108.73, 1174.66, 1244.51, 1318.51, 1396.91, 1479.98, 1567.98, 1661.22, // A6 to B6
            1760.00, 1864.66, 1975.53, 2093.00, 2217.46, 2349.32, 2489.02, 2637.02, 2793.83, 2959.96, 3135.96, 3322.44, // A7 to B7
            3520.00, 3729.31, 3951.07, 4186.01
};
PFont courier;

class track
{
  float freq; // Track frequency
  int num_samples = 0; 
  int max_samples;
  float[] samples = null;
  int index = 0;
  float sbuf = 0;
  float ampf = 1.0;
  float fq = filterq;
  
  track(float _freq, int _max_samples)
  {
    freq = _freq;
    max_samples = _max_samples;
    samples = new float[max_samples];
  }
  
  void store_sample(float s)
  {
    if (num_samples < max_samples)
      num_samples++;
    sbuf = fq*sbuf + (1-fq)*s;
    samples[index] = sbuf;
    index = (index + 1) % max_samples;
  }
  
  int get_num_samples()
  {
    return num_samples;
  }
    
  float[] get_samples()
  {
    if (num_samples == 0)
      return null;
    float[] res = new float[num_samples];
    if (num_samples < max_samples)
      for(int i=0; i<num_samples; i++)
        res[i] = constrain(ampf*samples[(index+i)%max_samples],-100,160);
    else
      for(int i=0; i<num_samples; i++)
        res[i] = constrain(ampf*samples[(index+i)%max_samples],-100,160);
    return res;
  }
}

float ease(float p, float g) {
    if (p < 0.5) 
        return 0.5 * pow(2*p, g);
    else
        return 1 - 0.5 * pow(2*(1 - p), g);
}

/**
 * A function to emulate a "jet" colormap similar to Matplotlib's.
 * It maps a normalized value (0 to 1) to a color ranging from blue to red, 
 * passing through cyan, green, and yellow.
 */
 
int jetColormap(float value) {
    value = constrain(value, 0, 1);
    float hue, saturation = 250, brightness = 255;

    if (value < 0.25f) {
        hue = map(value, 0.0f, 0.25f, 240, 180); // Blue to Cyan
    } else if (value < 0.5f) {
        hue = map(value, 0.25f, 0.5f, 180, 120); // Cyan to Green
    } else if (value < 0.75f) {
        hue = map(value, 0.5f, 0.75f, 120, 60);  // Green to Yellow
    } else {
        hue = map(value, 0.75f, 1.0f, 60, 0);    // Yellow to Red
    }

    colorMode(HSB, 360, 255, 255); // Set HSB mode with hue range [0, 360]
    int col = color(hue, saturation, brightness);
    colorMode(RGB, 255); // Reset color mode to default RGB for other parts of the program
    return col;
}

int jetColormap_old(float value) {
    float r, g, b;
    if (value < 0.25f) {
        r = 0;
        g = 4 * value;
        b = 1;
    } else if (value < 0.5f) {
        r = 0;
        g = 1;
        b = 1 - 4 * (value - 0.25f);
    } else if (value < 0.75f) {
        r = 4 * (value - 0.5f);
        g = 1;
        b = 0;
    } else {
        r = 1;
        g = 1 - 4 * (value - 0.75f);
        b = 0;
    }
    
    return color(r * 255, g * 255, b * 255);
}

public static double asinh(double a) {
    final double sign;
    // check the sign bit of the raw representation to handle -0
    if (Double.doubleToRawLongBits(a) < 0) {
        a = Math.abs(a);
        sign = -1.0d;
    } else {
        sign = 1.0d;
    }

    return sign * Math.log(Math.sqrt(a * a + 1.0d) + a);
}

// Function s(t) as defined
double s(double t) {
   return (t*Math.sqrt(1+alpha2*t*t) + asinh(alpha*t)/alpha)/2;
}

// Derivative of s(t) with respect to t: s'(t)
double ds(double t) {
    return Math.sqrt(1+alpha2*t*t);
}

// Function to approximate gamma(s) by numerically solving for t such that s(t) ≈ input sValue
PVector gamma(float t) 
{
    return new PVector(t, (float)(h*(1-t*t/d2)));
}

PVector gamma_tangent(double t) 
{
    return new PVector(1, -(float)(alpha2*t));
}

PVector gamma_normal(double t) 
{
    return new PVector((float)(alpha2*t), 1);
}

// Numerical method to solve for t given sValue, using Newton's method
double solve_for_t(double sValue) {
    if (Math.abs(sValue)<1e-5)
       return 0;
    final double eps = 1e-5; // Tolerance level for precision
    double t = sValue/ds(0);
    // Newton's method loop
    for (int i = 0; i < 10; i++) { // Limit iterations to prevent infinite loops
        double t1 = t - (s(t)-sValue)/ds(t);
        if (Math.abs(t1-t) <= eps*Math.abs(t1))
          return t1;
        t = t1;
    }
    throw new ArithmeticException("Newton's method did not converge for "+sValue);
}


float constrain_freq(float f)
{
  if (Float.isNaN(f) || f == Float.NEGATIVE_INFINITY)
      return 0.0f;
  else
      return f;
}

void settings()
{
  size(frameWidth, frameHeight, P3D);
  smooth(8);
}

void setup() {
  fac = width/600.0;
  frameRate(1000);

  //cam = new PeasyCam(this,  0*fac,  0*fac,  0*fac,  550*fac);
  // cam = new PeasyCam(this,  69.31465 *fac,  16.644705 *fac,  69.999 *fac,  411.6035061837172 *fac);
  // cam.setRotations( 0.8818339 ,  -0.69729364 ,  0.7555049 );
  //cam = new PeasyCam(this,  69.31465 *fac,  16.644705 *fac,  69.999 *fac,  411.603515625 *fac);
  //cam.setRotations( 0.8818339 ,  -0.69729364 ,  0.7555049 );
  //cam = new PeasyCam(this,  69.31465 *fac,  16.644705 *fac,  69.999 *fac,  411.603515625 *fac);
  //cam.setRotations( -0.0 ,  0.0 ,  -0.0 );
  // cam = new PeasyCam(this,  0.0 *fac,  0.0 *fac,  0.0 *fac,  650.0 *fac);
  // cam.setRotations( 1.5072619 ,  -0.6465731 ,  1.4747462 );
  //cam = new PeasyCam(this,  222.64006 *fac,  415.07318 *fac,  428.07187 *fac,  1273.063676001737 *fac);
  //cam.setRotations( 0.48491418 ,  -0.37045673 ,  0.2441992 );
  cam = new PeasyCam(this,  244.18518 *fac,  601.8838 *fac,  -264.27887 *fac,  1660.7410047886365 *fac);
  cam.setRotations( 1.0536625 ,  0.96110266 ,  -1.0081584 );
 
  courier = createFont("Courier New",24*fac,true);

  // Initialize the reader with the path to your binary file
  reader = new SpectrogramReader(dataPath("beethoven.bin"),false);

  // Initialize the tracks
  tracks = new track[track_freqs.length];
  for(int i=0; i<track_freqs.length; i++)
  {
    tracks[i] = new track(track_freqs[i], max_samples);
    for(int j=0; j<tracks[i].max_samples; j++)
      tracks[i].store_sample(-100.0);
    tracks[i].ampf = 1.2;
    tracks[i].fq = pow(filterq,1/fac);

  }
  tracks[0].ampf=2.1;
    
}

void draw() 
{
  t = 1.0*(frameCount-1)/(reader.totalPackets+padding);
  background(0);
  
  pushMatrix();
  rotateY(HALF_PI);
  //translate(0, height/2);
  
  final float rho = 7;
  //rotateY( constrain(map(1.0*frameCount/reader.totalPackets, 0, 0.1, 0, PI/3),0,PI/3) );
  //rotateY(-PI/3*sin((frameCount-1)*rho));
  //rotateX(-PI/12*sin((frameCount-1)*rho));
  // translate(0, 0, -track_dist*fac*tracks.length/4*sin((frameCount-1)*rho));
  float t1 = ease(constrain(rho*t,0,1), 2.2);
  float t2 = ease(constrain(rho*t-1, 0, 1), 2.2);
  float t3 = ease(constrain(rho*t-2, 0, 1), 2.2);
  float t4 = ease(constrain(rho*t-3, 0, 1), 2.2);
  float t5 = ease(constrain(rho*t-4, 0, 1), 2.2);
  float t6 = ease(constrain(rho*t-5, 0, 1), 2.2);
  float t7 = ease(constrain(rho*t-6, 0, 1), 2.2);
  
  rotateY( map(t1, 0, 1, 0, 1.05*PI) );
  rotateX( map(t1, 0, 1, 0, PI/6) );
  rotateZ( map(t1, 0, 1, 0, PI/6) );
  translate(0,-width*t2,-width*t2);
  rotateY( map(t3, 0, 1, 0, -PI) );
  rotateX( map(t3, 0, 1, 0, -PI/6) );
  translate(0, width*t4, 0);
  rotate(PI/4*t4);
  rotateX(PI/3*t5);
  translate(0,-width*t5,0);
  translate(-width*t6,width*t6,0);
  translate(4*width*t7,0,-4*0.6*width*t7);
  
  
  if (reading && reader.readNextPacket())
  {
    for(track t: tracks)
    {
      float sample = constrain_freq(reader.get(t.freq));
      float dBValue = 10 * (float) Math.log10(sample + 1e-10);
      if (Float.isNaN(dBValue) || dBValue == Float.NEGATIVE_INFINITY) {
         dBValue = -100.0f;
      }
      dBValue = constrain(dBValue, -100.0f, 160.0f);
      t.store_sample(dBValue);
    }
  }
  else
     reading = false;
     
  if (!reading)
  {
    padding_count++;
    for(track t: tracks)
    {
      t.store_sample(-100.0f);
    }
  }
  
 
  {
    float y = 0;
    float x = -width/2;
    float xw = width/2; // waterfall x
    float yw = (float)h;
    float z = -track_dist*fac*tracks.length/2;
    for(track t: tracks)
    {
      float sample = constrain_freq(reader.get(t.freq));
      float dBValue = 10 * (float) Math.log10(sample + 1e-10);
      if (Float.isNaN(dBValue) || dBValue == Float.NEGATIVE_INFINITY) {
         dBValue = -100.0f;
      }
      dBValue = constrain(dBValue, -100.0f, 160.0f);
      t.store_sample(dBValue);
      
      float[] vals = t.get_samples();
      int n = t.get_num_samples();
      
      // Draw the waterfall
      stroke(255);
      fill(0);
      strokeWeight(2*fac);
      beginShape();
      // vertex(x, y+100*fac, z);
      float last_val = -100.0;
      for(float i=0; i<(float)d; i++)
      {
        int k = int(map(i, 0, (float)d-1, min(n, max_samples)-1, 0));
        last_val = vals[k];
        stroke(jetColormap(map(last_val,-24,60,0,1)));
        PVector pg = gamma(i);
        PVector png = gamma_normal(i).normalize().mult(last_val*fac);
        pg = PVector.add(pg, png);
        vertex(xw-pg.x, yw-pg.y, z);
      }
      stroke(0);
      vertex(xw, yw, z);
      endShape();
      
      // Draw the pegs
      stroke(jetColormap(map(vals[n-1],-24,60,0,1)));
      noFill();
      beginShape();
      vertex(x+width,y+100*fac, z);
      vertex(x+width+hook_size*fac, y+100*fac, z);
      vertex(x+width+hook_size*fac, y-75*fac, z);
      endShape(); 
      
      // Draw the thin lines
      strokeWeight(0.5);
      stroke(80);
      noFill();
      line(xw, y+100*fac, z, xw+4*width, y+100*fac, z); // thin lines
      line(xw, yw, z, xw-4*width, yw, z);
      
      // Draw the peg's texts
      pushMatrix();
      translate(width/2,y,z);
      textSize(14*fac);
      rotateY(-HALF_PI);
      rotateZ(-HALF_PI);
      noStroke();
      fill(255);
      String ffreq = String.format("%.1f Hz", t.freq);
      text(ffreq,80*fac,5*fac,-hook_size*fac);
  
      popMatrix();
      z += track_dist*fac;
    }
  }
  popMatrix();
  
  // texts
  /* Need to compute the position of a frame sitting in the field of view
  pushMatrix(); pushStyle();
  float[] rot = cam.getRotations();
  float[] pos = cam.getPosition();
  //rotateX(rot[0]);
  //rotateY(-rot[1]);
  //rotateZ(rot[2]);
  
  translate(pos[0],pos[1],pos[2]-260*fac);
  fill(255);stroke(255);
  textAlign(CENTER, CENTER);
  textFont(courier);
  textSize(7*fac);
  text("@infinitymathart • 2024",0.0*width,0.0*height);
  popStyle(); popMatrix();  */
  
  if (recording)
  {
    saveFrame("/tmp/f/frame_#####.png");
    if (stamping) 
      println(frameCount,"/",reader.totalPackets);
    if (padding_count >= padding)
    {
      reader.close();
      exit();
    }
  }

}

void mousePressed()
{
  float[] rotations = cam.getRotations();
  float[] look = cam.getLookAt();
  double dist = cam.getDistance();
  stamping = false;
  println("..:");
  println("  cam = new PeasyCam(this, ", look[0]/fac, "*fac, ", look[1]/fac, "*fac, ", look[2]/fac, "*fac, ", dist/fac,"*fac);" );
  println("  cam.setRotations(", rotations[0], ", ", rotations[1], ", ", rotations[2],");\n\n");
}

/* License:
 *
 * Copyright (c) 2024 Waldeck Schutzer
 *
 * All rights reserved.
 *
 * This code after the template and the related animations are the property of the
 * copyright holder. Any reproduction, distribution, or use of this material,
 * in whole or in part, without the express written permission of the copyright
 * holder is strictly prohibited.
 */

SpectrogramReader.pde

import java.io.*;
import java.nio.*;

class SpectrogramReader {
    float[] frequencies; // Array of frequencies in Hz
    float[] minValues; // Min values for each frequency bin
    float[] maxValues; // Max values for each frequency bin
    int numFrequencies; // Number of frequency bins
    int sampleRate;
    int totalPackets;
    DataInputStream dataInputStream;
    boolean entireFileLoaded = false;
    float[][] allPackets; // Optional: store all packets if loading the entire file
    float[] currentPacket; // Current packet of frequency data

    // Constructor to load the binary file and initialize the frequency array
    SpectrogramReader(String filePath, boolean loadEntireFile) {
        try {
            // Open the binary file
            dataInputStream = new DataInputStream(new FileInputStream(filePath));
            
            // Read the header (128 bytes)
            sampleRate = Integer.reverseBytes(dataInputStream.readInt());
            numFrequencies = Integer.reverseBytes(dataInputStream.readInt());
            totalPackets = Integer.reverseBytes(dataInputStream.readInt());

            // Skip the padding bytes to complete the 128-byte header
            dataInputStream.skipBytes(128 - 12);

            // Read the frequencies array as little-endian floats
            frequencies = new float[numFrequencies];
            for (int i = 0; i < numFrequencies; i++) {
                frequencies[i] = readLittleEndianFloat(dataInputStream);
            }

            // Read the min and max values for each frequency bin
            minValues = new float[numFrequencies];
            maxValues = new float[numFrequencies];
            for (int i = 0; i < numFrequencies; i++)
                minValues[i] = readLittleEndianFloat(dataInputStream);
            for (int i = 0; i < numFrequencies; i++)
                maxValues[i] = readLittleEndianFloat(dataInputStream);
        
            // If requested, load the entire file into memory
            if (loadEntireFile) {
                //loadAllPackets();
                entireFileLoaded = true;
            }

            // Initialize the current packet array
            currentPacket = new float[numFrequencies];
        } catch (IOException e) {
            System.out.println("Error loading file: " + e.getMessage());
        }
    }

    // Method to read the next packet of spectrum data
    boolean readNextPacket() {
        if (entireFileLoaded) {
            return false; // If the entire file is loaded, use the stored data instead
        }

        try {
            // Read the spectrum values as little-endian floats
            for (int i = 0; i < numFrequencies; i++) {
                currentPacket[i] = readLittleEndianFloat(dataInputStream);
            }
            return true; // Packet read successfully
        } catch (EOFException e) {
            System.out.println("End of file reached.");
            return false; // Properly handle the end of file
        } catch (IOException e) {
            System.out.println("Error reading packet: " + e.getMessage());
            return false;
        }
    }

    // Method to get the value of the i-th frequency bin
    float get(int i) {
        if (i >= 0 && i < numFrequencies) {
            return currentPacket[i];
        } else {
            System.out.println("Index out of bounds: " + i);
            return 0.0f;
        }
    }

    // Method to get the minimum value of the i-th frequency bin
    float getMin(int i) {
        if (i >= 0 && i < numFrequencies) {
            return minValues[i];
        } else {
            System.out.println("Index out of bounds: " + i);
            return 0.0f;
        }
    }

    // Method to get the maximum value of the i-th frequency bin
    float getMax(int i) {
        if (i >= 0 && i < numFrequencies) {
            return maxValues[i];
        } else {
            System.out.println("Index out of bounds: " + i);
            return 0.0f;
        }
    }

    // Method to get the value for a specific frequency with interpolation
    float get(float freq) {
        int lowerIndex = binarySearch(frequencies, freq);
        if (lowerIndex == -1 || lowerIndex == numFrequencies - 1) {
            System.out.println("Frequency out of range: " + freq);
            return 0.0f;
        }

        int upperIndex = lowerIndex + 1;
        float lowerFreq = frequencies[lowerIndex];
        float upperFreq = frequencies[upperIndex];
        float lowerValue = currentPacket[lowerIndex];
        float upperValue = currentPacket[upperIndex];

        // Linear interpolation
        float proportion = (freq - lowerFreq) / (upperFreq - lowerFreq);
        return lowerValue + proportion * (upperValue - lowerValue);
    }

    // Method to get the minimum value for a specific frequency with interpolation
    float getMin(float freq) {
        int lowerIndex = binarySearch(minValues, freq);
        if (lowerIndex == -1 || lowerIndex == numFrequencies - 1) {
            System.out.println("Frequency out of range: " + freq);
            return 0.0f;
        }

        int upperIndex = lowerIndex + 1;
        float lowerFreq = frequencies[lowerIndex];
        float upperFreq = frequencies[upperIndex];
        float lowerValue = minValues[lowerIndex];
        float upperValue = minValues[upperIndex];

        float proportion = (freq - lowerFreq) / (upperFreq - lowerFreq);
        return lowerValue + proportion * (upperValue - lowerValue);
    }

    // Method to get the maximum value for a specific frequency with interpolation
    float getMax(float freq) {
        int lowerIndex = binarySearch(maxValues, freq);
        if (lowerIndex == -1 || lowerIndex == numFrequencies - 1) {
            System.out.println("Frequency out of range: " + freq);
            return 0.0f;
        }

        int upperIndex = lowerIndex + 1;
        float lowerFreq = frequencies[lowerIndex];
        float upperFreq = frequencies[upperIndex];
        float lowerValue = maxValues[lowerIndex];
        float upperValue = maxValues[upperIndex];

        float proportion = (freq - lowerFreq) / (upperFreq - lowerFreq);
        return lowerValue + proportion * (upperValue - lowerValue);
    }

    // Method to average the minimum values within a frequency range
    float getMin(float startFreq, float endFreq) {
        return getRangeValue(startFreq, endFreq, minValues);
    }

    // Method to average the maximum values within a frequency range
    float getMax(float startFreq, float endFreq) {
        return getRangeValue(startFreq, endFreq, maxValues);
    }

    // Method to average the values within a frequency range
    float get(float startFreq, float endFreq) {
        return getRangeValue(startFreq, endFreq, currentPacket);
    }

    // Helper method to average values within a range
    private float getRangeValue(float startFreq, float endFreq, float[] values) {
        if (startFreq > endFreq) {
            float temp = startFreq;
            startFreq = endFreq;
            endFreq = temp;
        }

        int startIndex = binarySearch(frequencies, startFreq);
        int endIndex = binarySearch(frequencies, endFreq);

        if (startIndex == -1 || endIndex == -1) {
            System.out.println("Frequency range out of bounds: " + startFreq + " - " + endFreq);
            return 0.0f;
        }

        if (startIndex > endIndex) {
            int temp = startIndex;
            startIndex = endIndex;
            endIndex = temp;
        }

        float sum = 0.0f;
        int count = 0;

        for (int i = startIndex; i <= endIndex; i++) {
            sum += values[i];
            count++;
        }

        return count == 0 ? 0.0f : sum / count;
    }

    // Helper method to perform binary search to find the closest frequency bin
    private int binarySearch(float[] array, float target) {
        int left = 0;
        int right = array.length - 1;

        while (left <= right) {
            int mid = left + (right - left) / 2;

            if (array[mid] == target) {
                return mid; // Exact match
            } else if (array[mid] < target) {
                left = mid + 1;
            } else {
                right = mid - 1;
            }
        }

        if (right < 0 || left >= array.length) {
            return -1;
        }

        return right;
    }

    // Method to read a little-endian float from the DataInputStream using ByteBuffer
    private float readLittleEndianFloat(DataInputStream dis) throws IOException {
        byte[] bytes = new byte[4];
        dis.readFully(bytes);
        return ByteBuffer.wrap(bytes).order(ByteOrder.LITTLE_ENDIAN).getFloat();
    }

    // Close the file stream when done
    void close() {
        try {
            dataInputStream.close();
        } catch (IOException e) {
            System.out.println("Error closing file: " + e.getMessage());
        }
    }
}