% A Matlab script to demonstrate parabolic interpolation of spectral peaks,
% which are subsequently used for modal synthesis. See the following link
% for details on the parabolic interpolation forumulae:
%
% - https://ccrma.stanford.edu/~jos/sasp/Quadratic_Interpolation_Spectral_Peaks.html
%
% by Gary Scavone, McGill University, 2019-2022.

[y, fs] = audioread( 'triangle.wav' );
N = 8192;             % length of FFT
iStart = 10000;       % starting sample of the FFT
Y = fft(y(iStart:N+iStart-1));

f = (0:N-1)*fs/N;
plot(f, 20*log10(abs(Y)))
xlim([0 15000])
grid
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)')

Y = abs(Y(1:ceil(N/2))); % throw away redundant half of fft
Xi = [];
threshold = 10^(30/20);  % linear gain corresponding to 30 dB
for n = 2:length(Y)-2    % search for peak indices above threshold
  if Y(n+1) <= Y(n) && Y(n-1) < Y(n) && Y(n) > threshold
    Xi = [Xi, n];
  end
end

% Perform parabolic interpolation of peaks
Yi = zeros( size( Xi ) );
for n = 1:length(Xi)
  alpha = Y(Xi(n)-1);
  beta = Y(Xi(n));
  gamma = Y(Xi(n)+1);
  Xii = 0.5*(alpha-gamma) / (alpha - 2*beta + gamma); % interpolated bin value (-1/2 to +1/2)
  Yi(n) = beta - 0.25 * (alpha-gamma) * Xii; % interpolated peak value
  if Xii > 0
    Xi(n) = f(Xi(n)) + Xii*(f(Xi(n)+1)-f(Xi(n))); % bin value in Hz
  else
    Xi(n) = f(Xi(n)) + Xii*(f(Xi(n))-f(Xi(n)-1)); % bin value in Hz
  end
end

% Choose first K peaks to resynthesize
K = 4;
Xi = Xi(1:K);
Yi = Yi(1:K);
for n = 1:K, xline(Xi(n), 'r--'); end
B = 0.2*[2.4 2 6.5 2.5]; % manually estimated 3dB bandwidths with fudge factor
r = exp(-pi*B/fs);
b0 = Yi;
a1 = -2*r.*cos(2*pi*Xi/fs);
a2 = r.^2;

% Create impulse response of filter bank
N = length(y);
x = [1; zeros(N-1, 1)];
y2 = filter( b0(1), [1 a1(1) a2(1)], x );
for n = 2:K
  y2 = y2 + filter( b0(n), [1 a1(n) a2(n)], x );
end

useResidual = false;
if useResidual
  % Derive residual of estimate from original
  Y = fft(y);
  Y2 = fft(y2);
  R = Y ./Y2;
  r = real(ifft(R));

  figure(2)
  N = 8192;             % length of FFT
  Rnew = fft(r(iStart:N+iStart-1));
  plot(f, 20*log10(abs(Rnew)))
  xlim([0 15000])
  grid
  xlabel('Frequency (Hz)')
  ylabel('Magnitude (dB)')
  title('Residual spectrum');

  % Now use residual to drive filters
  r = 30*r; % arbitrary level increase of residual
  y2 = filter( b0(1), [1 a1(1) a2(1)], r );
  for n = 2:K
    y2 = y2 + filter( b0(n), [1 a1(n) a2(n)], r );
  end
end

y2 = 0.5 * y2 / max(abs(y2));
sound([y; zeros(10000, 1); y2], fs)