After watching the excellent YouTube lecture: "Compressive Sensing" — https://youtu.be/rt5mMEmZHfs?si=gEvXWluUf3dBwT2Q this gist is an attempt on R translation of the lecture’s MATLAB code (works)

gistfile1.txt

# After watching the excellent YouTube lecture:
# "Compressive Sensing" — https://youtu.be/rt5mMEmZHfs?si=gEvXWluUf3dBwT2Q
# Which I found in refernces from even more awesom video:
# "X-ray backscatter with compressed sensing algorithm"
# https://www.youtube.com/watch?v=EuVgGrun1V0 (time to start building The Thing!-)
#
# Below is an R translation of the lecture’s MATLAB code.

# Original MATLAB line #1: clear all; close all; clc
# Equivalents of: clear all; close all; clc
rm(list = ls(all.names = TRUE))    # Clear workspace
graphics.off()                     # Close all plots
cat("\014")                        # Clear console

# Ensure 'pracma' is available (for DCT)
# NOTE: In pracma 2.4.4 the symbol `dct` is NOT exported, so pracma::dct() fails.
# We keep this note for context; the script below uses 'fftw' or 'gsignal' instead.
# If you still want pracma for other utilities, you can install it:
# if (!requireNamespace("pracma", quietly = TRUE)) install.packages("pracma")
# library(pracma)

# ---------------- Installation notes (Ubuntu) ----------------
# Preferred options for DCT:
#   Option A (fast, C-backed): R package 'fftw' which provides DCT/IDCT.
#     Ubuntu system dependency REQUIRED before installing the R package:
#       sudo apt-get update && sudo apt-get install --yes libfftw3-dev
#     Then in R: install.packages("fftw")
#   Option B (pure R, no system deps): R package 'gsignal' (has dct/idct).
#     In R: install.packages("gsignal")
# The code below automatically uses fftw if available, otherwise gsignal,
# and finally a portable FFT-based fallback (scaling may differ by a constant).
# -------------------------------------------------------------

# ----------------------------- Parameters and signal -------------------------
# Original MATLAB line #3: n = 5000;
n <- 5000
# Original MATLAB line #4: t = linspace(0, 1/8, n);
t <- seq(0, 1/8, length.out = n)   # MATLAB linspace(0, 1/8, n)
# Original MATLAB line #5: f = sin(1394*pi*t) + sin(3266*pi*t);
f <- sin(1394 * pi * t) + sin(3266 * pi * t)

# ------------------------------ Robust DCT-II helper (final) -----------------
dct_1d <- function(x) {
  # Validate input up front
  stopifnot(is.numeric(x), length(x) > 0)
  
  # Preferred: fftw::DCT (Type-II as in MATLAB dct)
  if (requireNamespace("fftw", quietly = TRUE)) {
    return(fftw::DCT(x, type = 2))
  }
  
  # Alternative: gsignal::dct (defaults to Type-II)
  if (requireNamespace("gsignal", quietly = TRUE)) {
    return(gsignal::dct(x))
  }
  
  # Fallback: Compute DCT-II via mirrored FFT (scaling may differ by a constant)
  n_local <- length(x)
  y <- c(x, rev(x))                 # Mirror to length 2n
  Y <- fft(y)                       # Complex FFT
  k <- 0:(n_local - 1)
  return(Re(2 * Y[k + 1] * exp(-1i * pi * k / (2 * n_local))))
}
# ---------------------------- End robust DCT-II helper -----------------------

# Original MATLAB line #8: ft = dct(f);
ft <- dct_1d(f)

# ============================= Single-window layout ==========================
# (We switch to a single device with two stacked panels: top = time signal (+samples, +recon),
#  bottom = DCT spectrum (+recon spectrum). This avoids device juggling in RStudio.)

# Prevent "figure margins too large" in small Plots pane; open external window if needed
ensure_plotting_area <- function(min_w_in = 6.5, min_h_in = 6.0) {
  sz <- try(grDevices::dev.size("in"), silent = TRUE)
  need_external <- inherits(sz, "try-error") || any(is.na(sz)) ||
    sz[1] < min_w_in || sz[2] < min_h_in
  if (need_external) {
    if (Sys.getenv("RSTUDIO") == "1" && capabilities("X11")) {
      grDevices::dev.new(width = max(min_w_in, 7), height = max(min_h_in, 6), noRStudioGD = TRUE)
    } else if (capabilities("X11")) {
      grDevices::x11(width = max(min_w_in, 7), height = max(min_h_in, 6))
    }
  }
}
ensure_plotting_area()

op <- par(no.readonly = TRUE)
par(mfrow = c(2, 1),
    mar   = c(3.0, 3.2, 1.6, 0.8),
    mgp   = c(1.8, 0.5, 0),
    tcl   = -0.25,
    xaxs  = "i", yaxs = "i",
    cex.axis = 0.9, cex.lab = 0.95)

# Helper to render both panels in one call (idempotent re-render)
render_panels <- function(t, f_original, ft_original,
                          samples = NULL, f_recon = NULL, ft_recon = NULL) {
  # ---- Top panel: time-domain signal ----
  # Original MATLAB line #6: plot(t, f, 'LineWidth', [2])
  plot(t, f_original, type = "l", lwd = 2,
       xlab = "t (s)", ylab = "f(t)",
       main = "Sum of Two Sine Waves (Time Domain)")
  if (!is.null(samples)) {
    # Original MATLAB line #20: hold on, plot(tr, fr, 'ro', 'LineWidth', [2])
    points(samples$tr, samples$fr, pch = 1, col = "red", lwd = 2)
  }
  if (!is.null(f_recon)) {
    # (Used later: overlay reconstruction)
    lines(t, f_recon, lwd = 2, col = "blue")
  }
  legend_entries <- c("Original")
  legend_cols    <- c("black")
  if (!is.null(samples)) {
    legend_entries <- c(legend_entries, "Samples")
    legend_cols    <- c(legend_cols, "red")
  }
  if (!is.null(f_recon)) {
    legend_entries <- c(legend_entries, "Reconstruction")
    legend_cols    <- c(legend_cols, "blue")
  }
  legend("topright", legend = legend_entries, col = legend_cols, lwd = 2, bty = "n")
  
  # ---- Bottom panel: DCT spectrum ----
  # Original MATLAB line #9: figure(2)
  # Original MATLAB line #10: plot(ft, 'LineWidth', [2])
  plot(ft_original, type = "l", lwd = 2,
       xlab = "k", ylab = "Amplitude",
       main = "DCT Spectrum (Type-II)")
  if (!is.null(ft_recon)) {
    lines(ft_recon, lwd = 2, col = "blue")
    legend("topright", legend = c("Original DCT", "Reconstructed DCT"),
           col = c("black", "blue"), lwd = 2, bty = "n")
  }
}

# ------- Adding transparancy (begin) ---------
# =========================== Multi-layer renderer with transparency ===========================
# Color palette with transparency (works in RStudio/X11 devices)
col_original <- grDevices::adjustcolor("gray30",     alpha.f = 0.90)
col_samples  <- grDevices::adjustcolor("red",        alpha.f = 0.65)
col_wrong    <- grDevices::adjustcolor("blue",       alpha.f = 0.80)
col_cvx      <- grDevices::adjustcolor("darkgreen",  alpha.f = 0.98)

# Render both panels with multiple overlays; overlay_mode = "overlay" or "staggered"
render_panels_multi <- function(t, f_original, ft_original,
                                samples = NULL,
                                f_wrong = NULL,  ft_wrong = NULL,
                                f_cvx  = NULL,   ft_cvx  = NULL,
                                overlay_mode = c("overlay", "staggered")) {
  overlay_mode <- match.arg(overlay_mode)
  n_local <- length(t)
  
  # ---- Top panel: time domain ----
  plot(t, f_original, type = "l", lwd = 2, col = col_original,
       xlab = "t (s)", ylab = "f(t)", main = "Sum Of Two Sine Waves (Time Domain)")
  
  # Samples (semi-transparent filled circles)
  if (!is.null(samples)) {
    points(samples$tr, samples$fr, pch = 21, lwd = 1.1,
           col = col_samples, bg = grDevices::adjustcolor("red", alpha.f = 0.35), cex = 0.8)
  }
  
  if (!is.null(f_wrong) || !is.null(f_cvx)) {
    if (overlay_mode == "overlay") {
      if (!is.null(f_wrong)) lines(t, f_wrong, lwd = 2.2, col = col_wrong)
      if (!is.null(f_cvx))   lines(t, f_cvx,   lwd = 2.6, col = col_cvx)
    } else { # "staggered"
      i1 <- 1:floor(n_local / 3)
      i23_start <- floor(n_local / 3) + 1
      if (!is.null(f_wrong)) lines(t[i1], f_wrong[i1], lwd = 2.4, col = col_wrong)
      if (!is.null(f_cvx))   lines(t[i23_start:n_local], f_cvx[i23_start:n_local], lwd = 3.0, col = col_cvx)
    }
  }
  
  # Legend (top)
  legend_lbl <- c("Original"); legend_col <- c(col_original); legend_lwd <- c(2); legend_pch <- c(NA)
  if (!is.null(samples)) { legend_lbl <- c(legend_lbl, "Samples"); legend_col <- c(legend_col, col_samples); legend_lwd <- c(legend_lwd, NA); legend_pch <- c(legend_pch, 21) }
  if (!is.null(f_wrong))  { legend_lbl <- c(legend_lbl, if (overlay_mode=="overlay") "Wrong Recon" else "Wrong Recon (Left 1/3)"); legend_col <- c(legend_col, col_wrong); legend_lwd <- c(legend_lwd, 2.2); legend_pch <- c(legend_pch, NA) }
  if (!is.null(f_cvx))    { legend_lbl <- c(legend_lbl, if (overlay_mode=="overlay") "CVX Recon" else "CVX Recon (Right 2/3)");   legend_col <- c(legend_col, col_cvx);   legend_lwd <- c(legend_lwd, 2.6); legend_pch <- c(legend_pch, NA) }
  legend("topright", legend = legend_lbl, col = legend_col, lwd = legend_lwd, pch = legend_pch, pt.cex = 0.9, bty = "n")
  
  # ---- Bottom panel: DCT spectrum ----
  par(mfg = c(2, 1))
  plot(ft_original, type = "l", lwd = 2, col = col_original,
       xlab = "k", ylab = "Amplitude", main = "DCT Spectrum (Type-II)")
  if (!is.null(ft_wrong)) lines(ft_wrong, lwd = 2.2, col = col_wrong)
  if (!is.null(ft_cvx))   lines(ft_cvx,   lwd = 2.6, col = col_cvx)
  
  legend2_lbl <- c("Original DCT"); legend2_col <- c(col_original)
  if (!is.null(ft_wrong)) { legend2_lbl <- c(legend2_lbl, "Pseudo-inverse Spectrum (x)"); legend2_col <- c(legend2_col, col_wrong) }
  if (!is.null(ft_cvx))   { legend2_lbl <- c(legend2_lbl, "CVX Spectrum");                 legend2_col <- c(legend2_col, col_cvx) }
  legend("topright", legend = legend2_lbl, col = legend2_col, lwd = 2, bty = "n")
}
# ========================= End multi-layer renderer with transparency ==========================

# ------- Adding transparency (end) -----------

# Original MATLAB lines #6–#10 effect (first render)
render_panels(t, f, ft)

# ============================ Sampling block (top) ===========================
# Original MATLAB line #13: m = 200;      (Note: the lecture then indexes temp(1:500))
# Original MATLAB line #14: temp = randperm(5000);
# Original MATLAB line #15: ind = temp(1:500);
# Original MATLAB line #16: tr = t(ind);
# Original MATLAB line #17: fr = sin(1394*pi*tr) + sin(3266*pi*tr);
# Original MATLAB line #20: hold on, plot(tr, fr, 'ro', 'LineWidth', [2])
set.seed(42)                       # Make sampling reproducible (choose any integer)
m <- 500                           # Number of samples to take (choose 500 to match the plot)
temp <- sample.int(n)              # MATLAB randperm(n)
ind <- temp[1:m]                   # First m indices
tr <- t[ind]                       # Sampled time instants
# NOTE: For overlay reliability, sample 'f' by index rather than recomputing sin()
fr <- f[ind]                       # Sampled signal values

# Re-render both panels so that the red circles appear on the TOP panel
samples <- list(tr = tr, fr = fr)
render_panels(t, f, ft, samples = samples)

# ======================= "Wrong" pseudo-inverse approach (as in lecture) =======================
# Original MATLAB line #21: D = dct(eye(n, n));
# Original MATLAB line #22: A = D(ind, :);
# Build only the required rows of the DCT-II transform matrix (no full n×n allocation).
if (!exists("build_dct_rows")) {
  build_dct_rows <- function(n, rows) {
    stopifnot(length(rows) > 0, all(rows >= 1), all(rows <= n))
    k  <- as.numeric(rows - 1)         # 0..n-1 for requested rows
    nn <- 0:(n - 1)                    # time indices 0..n-1
    # Unnormalized DCT-II kernel; scale mismatch is fine for this demo
    t(2 * cos((pi / n) * outer(nn + 0.5, k, "*")))
  }
}
A_wrong <- build_dct_rows(n, ind)      # m × n matrix analogous to D(ind, :)

# Original MATLAB line #25: x = pinv(A) * fr';
# Moore–Penrose pseudoinverse (MASS::ginv). Install once if needed:
if (!requireNamespace("MASS", quietly = TRUE)) install.packages("MASS")
x_wrong <- as.numeric(MASS::ginv(A_wrong) %*% matrix(fr, ncol = 1))

# Original MATLAB line #26: sig1 = dct(x);
sig1_wrong <- dct_1d(x_wrong)          # Uses the robust helper defined earlier

# Original MATLAB line #31: figure(3)
# Original MATLAB line #32: plot(t, f, t, sig1)
# Top panel: draw original time signal and overlay the (wrong) reconstruction in blue.
render_panels(t, f, ft, samples = samples, f_recon = sig1_wrong)

# Original MATLAB line #34: figure(4)
# Original MATLAB line #35: plot(ft, 'b'), hold on, plot(x, 'r')
# Bottom panel: overlay x (the “spectrum” from pseudo-inverse) in red.
par(mfg = c(2, 1))
lines(x_wrong, lwd = 2, col = "red")
legend("topright",
       legend = c("Original DCT (ft)", "Pseudo-inverse spectrum (x)"),
       col    = c("black", "red"), lwd = 2, bty = "n")

# ======================= Final (correct) sparse recovery via CVXR =============================
# NOTE: Lecturer edited starting from MATLAB line #22 (see below). We mirror that here in R.
#
# Ubuntu install hint (one-time): CVXR usually installs from CRAN without extra system libs.
# If compilation errors pop up about compilers, ensure you have build tools:
#   sudo apt-get update && sudo apt-get install --yes build-essential pkg-config libopenblas-dev liblapack-dev
# Then in R:
#   install.packages("CVXR")

if (!requireNamespace("CVXR", quietly = TRUE)) install.packages("CVXR")

# Provide idct_1d if not defined yet (paired with our dct_1d).
if (!exists("idct_1d")) {
  idct_1d <- function(X) {
    stopifnot(is.numeric(X), length(X) > 0)
    if (requireNamespace("fftw", quietly = TRUE)) {
      return(fftw::IDCT(X, type = 2))
    }
    if (requireNamespace("gsignal", quietly = TRUE)) {
      return(gsignal::idct(X))
    }
    # Fallback IDCT-III matching the common (unnormalized) DCT-II convention
    n <- length(X)
    nn <- 0:(n - 1)
    k  <- 0:(n - 1)
    M  <- cos((pi / n) * outer(nn + 0.5, k, "*"))
    M[, 1] <- M[, 1] * 0.5
    as.numeric((M %*% X) / n)
  }
}

# Build rows of the IDCT synthesis matrix for the sampled time indices:
# This implements A := S * IDCT_matrix, where S picks rows 'ind' (time samples).
build_idct_rows <- function(n, time_rows) {
  stopifnot(length(time_rows) > 0, all(time_rows >= 1), all(time_rows <= n))
  nn <- as.numeric(time_rows - 1)     # time indices 0..n-1 to KEEP (rows)
  k  <- 0:(n - 1)                     # DCT coefficient indices 0..n-1 (columns)
  M  <- cos((pi / n) * outer(nn + 0.5, k, "*"))
  M[, 1] <- M[, 1] * 0.5              # 1/2 on the DC term
  M / n                               # Unnormalized IDCT scaling
}

# Final section: Implements the lecturer’s intended approach.
# L1 minimization (Basis Pursuit) recovers a signal that is sparse in the DCT basis from few time-domain samples.
# The recovered time-domain signal and its DCT spectrum are correctly overlaid on the top and bottom panels, respectively.

# ======================= Final part (correct recovery) =======================
# Note: Lecturer got back and edited Matlab starting from line #22.
# Solves: minimize ||x||_1  subject to  A * x = fr,  then f_recon = idct(x).

# Original MATLAB line #22: A = D(ind, :);             (In the correct model A = S * IDCT_matrix)
A_cvx <- build_idct_rows(n, ind)                       # m × n

# Original MATLAB lines #25–#30:
#   cvx_begin
#     variable x(n);
#     minimize(norm(x,1));
#     subject to
#       A*x == fr';
#   cvx_end

# install.packages(c("CVXR","ECOSolveR","SCS","OSQP"), dependencies = TRUE, repos = "https://cloud.r-project.org")
# Ensure CVXR and solvers are available (ECOS/SCS/OSQP). Install if missing.
# Ubuntu one-time note: if compilation fails, run in terminal:
#   sudo apt-get update && sudo apt-get install --yes build-essential pkg-config libopenblas-dev liblapack-dev
if (!requireNamespace("CVXR", quietly = TRUE)) {
  install.packages(c("CVXR","ECOSolveR","SCS","OSQP"), dependencies = TRUE, repos = "https://cloud.r-project.org")
}

# ---------------- Ensure CVXR and its dependencies are available ----------------
# Ubuntu one-time system deps for Rmpfr (required by CVXR):
#   sudo apt-get update && sudo apt-get install --yes libgmp-dev libmpfr-dev
# Then in R (once): install.packages(c("gmp","Rmpfr","CVXR","ECOSolveR","SCS","OSQP"), dependencies = TRUE)

ensure_cvxr_available <- function() {
  # Try to load CVXR; if missing, attempt install (R-side only). System libs must be preinstalled.
  if (!requireNamespace("CVXR", quietly = TRUE)) {
    # Ensure Rmpfr dependency chain is present
    if (!requireNamespace("gmp", quietly = TRUE)) {
      message("Installing 'gmp' (requires libgmp-dev on Ubuntu)...")
      install.packages("gmp", repos = "https://cloud.r-project.org")
    }
    if (!requireNamespace("Rmpfr", quietly = TRUE)) {
      message("Installing 'Rmpfr' (requires libmpfr-dev on Ubuntu)...")
      install.packages("Rmpfr", repos = "https://cloud.r-project.org")
    }
    message("Installing 'CVXR' and common solvers (ECOS/SCS/OSQP)...")
    install.packages(c("CVXR","ECOSolveR","SCS","OSQP"), dependencies = TRUE, repos = "https://cloud.r-project.org")
  }
  
  # Pick a stable solver
  solvers <- tryCatch(CVXR::installed_solvers(), error = function(e) character())
  if (length(solvers) == 0) return(NA_character_)
  if ("ECOS" %in% solvers) return("ECOS")
  if ("SCS"  %in% solvers) return("SCS")
  if ("OSQP" %in% solvers) return("OSQP")
  solvers[[1]]
}

use_solver <- ensure_cvxr_available()

# Check use_solver is set
if (!use_solver %in% tryCatch(CVXR::installed_solvers(), error = function(e) character())) {
  alt <- intersect(c("SCS","OSQP"), tryCatch(CVXR::installed_solvers(), error = function(e) character()))
  use_solver <- if (length(alt)) alt[[1]] else NA_character_
}

x_var <- CVXR::Variable(n)
prob  <- CVXR::Problem(CVXR::Minimize(CVXR::norm1(x_var)),
                       list(A_cvx %*% x_var == matrix(fr, ncol = 1)))

# Use explicit solver if available (avoids first-run autodetect hiccups)
res <- if (!is.na(use_solver)) {
  CVXR::solve(prob, solver = use_solver)
} else {
  CVXR::solve(prob)  # fallback if solver detection failed
}

# Sparse DCT-coefficient estimate
x_cvx <- as.numeric(res$getValue(x_var))

# Time-domain reconstruction via IDCT
f_cvx  <- idct_1d(x_cvx)
ft_cvx <- dct_1d(f_cvx)

# Small normalization guard: libraries use different DCT scalings.
# Rescale f_cvx to match the observed samples at 'ind' (least-squares one-parameter fit).
alpha  <- sum(fr * f_cvx[ind]) / max(1e-12, sum(f_cvx[ind]^2))
f_cvx  <- alpha * f_cvx
ft_cvx <- alpha * ft_cvx

# ---------------- Overlay the correct reconstruction in dark green ----------------
# Ensure the two-panel layout is still active; then add lines to both panels:
par(mfg = c(1, 1))  # Top panel (time domain)
lines(t, f_cvx, lwd = 2, col = "darkgreen")

par(mfg = c(2, 1))  # Bottom panel (DCT spectrum)
lines(ft_cvx, lwd = 2, col = "darkgreen")

# Refresh legends to include all layers (Original, Samples, Wrong, CVX)
par(mfg = c(1, 1))
legend("topright",
       legend = c("Original", "Samples", "Wrong recon", "CVX recon"),
       col    = c("black",    "red",      "blue",       "darkgreen"),
       lwd    = c(2,          NA,         2,            2),
       pch    = c(NA,         1,          NA,           NA),
       bty    = "n")

par(mfg = c(2, 1))
legend("topright",
       legend = c("Original DCT (ft)", "Pseudo-inverse spectrum (x)", "CVX spectrum"),
       col    = c("black",             "red",                         "darkgreen"),
       lwd    = c(2,                   2,                              2),
       bty    = "n")

# ========================= Final transparent refresh (top + bottom) ============================
render_panels_multi(t, f, ft,
                    samples = samples,
                    f_wrong = sig1_wrong,  ft_wrong = x_wrong,
                    f_cvx  = f_cvx,        ft_cvx  = ft_cvx,
                    overlay_mode = "staggered")   # try "overlay" if you prefer full-length overlays

# ===============================================================================================

# (Optional) Restore par settings at the very end if you will plot further elsewhere
# try(par(op), silent = TRUE)