% x = input signal
% fs = sample rate
% ti = input duration (window length) in seconds
% to = output duration in seconds
% d = average grain density
% w = window type (1=triangular, 2=Parzen, 3=Nuttal, 4=sinusoidal, 5=Welch)
function y = freeze(x, fs, ti, to, d, w)

  % if x is row, transponse it
  if (columns(x) == 1)
    x = x';
  endif
  
  % compute durations in samples
  N = round(ti * fs);
  No = round(to * fs);

  % generate velvet noise
  nx = rand(1, No);
  low = d / (2 * N);
  high = 1 - low;
  n(nx<=high) = 0;
  n(nx>high) = 1;
  n(nx<low) = -1;

  % compute input level
  tau = 0.06514417228548776;	% 300 ms for 99% excursion
  mA1 = exp(-1 / (fs * tau));
  B0 = 1 - mA1;
  lx = filter(B0, [1, -mA1], 20 * log10(max(abs(x), 1e-6)), -120);
  lx = lx(end);

  % granulate input
  x = x(end-N+1:end) .* window(0:N-1, N, w);
  
  % convolution
  y = fftconv(x, n);
  y = y(1:No);

  % compensate output level
  ly = filter(B0, [1, -mA1], 20 * log10(max(abs(y), 1e-6)), 20 * log10(d / ti));
  y = y .* (10 .^ (min(lx - ly, 6) / 20));

endfunction

% n in [0,N)
function y = window(n, N, w)
  if (w == 1)
    % triangular
    y = 1-abs(2*n/N-1);
  elseif (w == 2)
    % Parzen
    n = n-N/2;
    y(abs(n)<N/4) = 1-6*(2*n(abs(n)<N/4)/N).^2.*(1-2*abs(n(abs(n)<N/4))/N);
    y(abs(n)>=N/4) = 2*(1-2*abs(n(abs(n)>=N/4))/N).^3;
  elseif (w == 3)
    % Nuttal
    a0=0.355768;
    a1=0.487369;
    a2=0.144232;
    a3=0.012604;
    y = a0 - a1 * cos(2*pi*n/(N-1)) + a2 * cos(4*pi*n/(N-1)) - a3 * cos(6*pi*n/(N-1));
  elseif (w == 4)
    % sin
    y = sin(pi*n/(N-1));
  else
    % Welch
    y = 1-((n-(N-1)/2)/((N-1)/2)).^2;
  endif
endfunction
