{"spec_id":"spectrogram-mel","library":"ggplot2","language":"r","code":"#' anyplot.ai\n#' spectrogram-mel: Mel-Spectrogram for Audio Analysis\n#' Library: ggplot2 3.5.1 | R 4.4.1\n#' Quality: 88/100 | Created: 2026-06-03\n\nlibrary(ggplot2)\nlibrary(scales)\nlibrary(ragg)\n\nset.seed(42)\n\n# --- Theme tokens ---\nTHEME       <- Sys.getenv(\"ANYPLOT_THEME\", \"light\")\nPAGE_BG     <- if (THEME == \"light\") \"#FAF8F1\" else \"#1A1A17\"\nELEVATED_BG <- if (THEME == \"light\") \"#FFFDF6\" else \"#242420\"\nINK         <- if (THEME == \"light\") \"#1A1A17\" else \"#F0EFE8\"\nINK_SOFT    <- if (THEME == \"light\") \"#4A4A44\" else \"#B8B7B0\"\nINK_MUTED   <- if (THEME == \"light\") \"#6B6A63\" else \"#A8A79F\"\n\n# --- Audio parameters ---\nsample_rate <- 22050L\nduration    <- 4.0\nn_fft       <- 2048L\nhop_length  <- 512L\nn_mels      <- 128L\n\n# --- Synthesize audio: C-major arpeggio with 8 harmonics + percussive transients ---\nn_samples  <- as.integer(sample_rate * duration)\nt_vec      <- seq(0, duration, length.out = n_samples + 1L)[seq_len(n_samples)]\nnote_freqs <- c(261.63, 329.63, 392.00, 523.25, 659.25, 523.25, 392.00, 329.63)\nnote_dur   <- duration / length(note_freqs)\n\naudio    <- numeric(n_samples)\nharm_amp <- 0.5 * 0.6^(0:7)  # 8 harmonics with exponential amplitude decay\n\nfor (i in seq_along(note_freqs)) {\n  t0  <- (i - 1L) * note_dur\n  idx <- which(t_vec >= t0 & t_vec < t0 + note_dur)\n  f   <- note_freqs[i]\n  dt  <- t_vec[idx] - t0\n  env <- pmax(0, (1 - exp(-300 * dt)) * exp(-4 * dt))\n  wave <- Reduce(\"+\", lapply(seq_along(harm_amp), function(h) {\n    harm_amp[h] * sin(2 * pi * h * f * t_vec[idx])\n  }))\n  # Percussive onset burst: broadband noise decaying over ~25 ms\n  burst_len <- min(as.integer(0.025 * sample_rate), length(idx))\n  burst_env <- c(exp(-150 * dt[seq_len(burst_len)]), numeric(length(idx) - burst_len))\n  audio[idx] <- env * wave + 0.18 * burst_env * rnorm(length(idx))\n}\naudio <- audio + 0.012 * rnorm(n_samples)\n\n# --- STFT: Hann-windowed power spectrogram ---\nn_fft_half <- n_fft %/% 2L + 1L\nhann_win   <- 0.5 * (1 - cos(2 * pi * seq(0L, n_fft - 1L) / (n_fft - 1L)))\nn_frames   <- floor((n_samples - n_fft) / hop_length) + 1L\n\nstft_power <- matrix(0.0, nrow = n_fft_half, ncol = n_frames)\nfor (i in seq_len(n_frames)) {\n  s <- (i - 1L) * hop_length + 1L\n  stft_power[, i] <- Mod(fft(audio[s:(s + n_fft - 1L)] * hann_win)[seq_len(n_fft_half)])^2\n}\n\n# --- Mel filterbank ---\nhz_to_mel <- function(f) 2595 * log10(1 + f / 700)\nmel_to_hz <- function(m) 700 * (10^(m / 2595) - 1)\n\nf_min     <- 80.0\nf_max     <- as.numeric(sample_rate) / 2.0\nmel_pts   <- seq(hz_to_mel(f_min), hz_to_mel(f_max), length.out = n_mels + 2L)\nhz_pts    <- mel_to_hz(mel_pts)\nfft_freqs <- seq(0, f_max, length.out = n_fft_half)\n\nmel_fb <- matrix(0.0, nrow = n_mels, ncol = n_fft_half)\nfor (m in seq_len(n_mels)) {\n  rising  <- (fft_freqs - hz_pts[m])      / (hz_pts[m + 1L] - hz_pts[m])\n  falling <- (hz_pts[m + 2L] - fft_freqs) / (hz_pts[m + 2L] - hz_pts[m + 1L])\n  mel_fb[m, ] <- pmax(0.0, pmin(rising, falling))\n}\n\n# --- Mel spectrogram in dB, normalized to 0 dB peak ---\nmel_spec    <- mel_fb %*% stft_power\nmel_spec_db <- 10 * log10(mel_spec + 1e-10)\nmel_spec_db <- mel_spec_db - max(mel_spec_db)\n\n# --- Long-format data frame ---\nmel_centers <- mel_to_hz(mel_pts[2:(n_mels + 1L)])\ntime_axis   <- ((seq_len(n_frames) - 1L) * hop_length + n_fft / 2L) / sample_rate\n\ngrid_idx <- expand.grid(mel_band = seq_len(n_mels), time_idx = seq_len(n_frames))\ndf <- data.frame(\n  time_s   = time_axis[grid_idx$time_idx],\n  mel_band = grid_idx$mel_band,\n  db       = mel_spec_db[cbind(grid_idx$mel_band, grid_idx$time_idx)]\n)\n\n# Y-axis: key frequency labels at representative mel-band positions\nkey_freqs  <- c(100, 250, 500, 1000, 2000, 4000, 8000)\nkey_bands  <- sapply(key_freqs, function(f) which.min(abs(mel_centers - f)))\nkey_labels <- ifelse(key_freqs >= 1000, paste0(key_freqs / 1000, \"k Hz\"), paste0(key_freqs, \" Hz\"))\n\n# Annotation reference positions\nnote_onsets <- (seq_along(note_freqs) - 1L) * note_dur\nband_1k     <- which.min(abs(mel_centers - 1000))\nband_2k     <- which.min(abs(mel_centers - 2000))\nt_max       <- max(time_axis)\n\ntitle_str <- \"spectrogram-mel · r · ggplot2 · anyplot.ai\"\n\n# --- Plot ---\np <- ggplot(df, aes(x = time_s, y = mel_band, fill = db)) +\n  geom_tile() +\n  # Note onset markers — reveal rhythmic structure of the arpeggio\n  geom_vline(\n    xintercept = note_onsets[-1],\n    color      = INK_MUTED,\n    linewidth  = 0.3,\n    linetype   = \"dotted\"\n  ) +\n  # Perceptual boundary: 1 kHz separates fundamental region from overtones\n  geom_hline(\n    yintercept = band_1k,\n    color      = INK_SOFT,\n    linewidth  = 0.45,\n    linetype   = \"dashed\"\n  ) +\n  annotate(\"text\",\n    x = t_max * 0.97, y = band_1k + 2.5,\n    label = \"1 kHz\", color = INK_SOFT,\n    size = 2.3, hjust = 1\n  ) +\n  # Frequency region labels for interpretive guidance\n  annotate(\"text\",\n    x = 0.10, y = 5,\n    label = \"Fundamentals\", color = INK_MUTED,\n    size = 2.3, hjust = 0, fontface = \"italic\"\n  ) +\n  annotate(\"text\",\n    x = 0.10, y = band_2k + 4,\n    label = \"Overtones\", color = INK_MUTED,\n    size = 2.3, hjust = 0, fontface = \"italic\"\n  ) +\n  scale_fill_gradient(\n    name   = \"dB\",\n    low    = \"#009E73\",\n    high   = \"#4467A3\",\n    limits = c(-80, 0),\n    oob    = scales::squish,\n    breaks = c(0, -20, -40, -60, -80)\n  ) +\n  scale_x_continuous(\n    name   = \"Time (s)\",\n    expand = c(0, 0)\n  ) +\n  scale_y_continuous(\n    name   = \"Frequency (Hz)\",\n    breaks = key_bands,\n    labels = key_labels,\n    expand = c(0, 0)\n  ) +\n  guides(fill = guide_colorbar(barheight = 7, barwidth = 0.7, ticks = TRUE)) +\n  labs(\n    title    = title_str,\n    subtitle = \"C-major arpeggio · 8 harmonics · mel scale compresses perceptual distances\"\n  ) +\n  theme_minimal(base_size = 8) +\n  theme(\n    plot.background   = element_rect(fill = PAGE_BG,     color = PAGE_BG),\n    panel.background  = element_rect(fill = PAGE_BG,     color = NA),\n    panel.grid.major  = element_blank(),\n    panel.grid.minor  = element_blank(),\n    panel.border      = element_blank(),\n    axis.title        = element_text(color = INK,        size = 10),\n    axis.text         = element_text(color = INK_SOFT,   size = 8),\n    axis.line         = element_line(color = INK_SOFT,   linewidth = 0.4),\n    plot.title        = element_text(color = INK,        size = 12, face = \"bold\"),\n    plot.subtitle     = element_text(color = INK_SOFT,   size = 8),\n    legend.background = element_rect(fill = ELEVATED_BG, color = INK_SOFT, linewidth = 0.3),\n    legend.text       = element_text(color = INK_SOFT,   size = 8),\n    legend.title      = element_text(color = INK,        size = 10),\n    plot.margin       = margin(16, 16, 16, 16)\n  )\n\n# --- Save ---\nggsave(\n  filename = sprintf(\"plot-%s.png\", THEME),\n  plot     = p,\n  device   = ragg::agg_png,\n  width    = 8,\n  height   = 4.5,\n  units    = \"in\",\n  dpi      = 400\n)\n"}