Ideal Spectral Interpolation

The ideal interpolation kernel for a uniformly sampled spectrum computed by a length $N$ DFT was derived in §4.2. The result was the aliased sinc function, derived in Eq.$\,$(1.3):

$$\hbox{asinc}_N(\omega) \isdef \frac{1}{N}\frac{\sin(\omega N/2)}{\sin(\omega/2)}$$

This result was quickly derived by performing an inverse DFT followed by a DTFT.

Figure F.1 lists a matlab function for performing ideal spectral interpolation. Note that in this implementation, the DFT is treated as causal rather than zero phase, so the precise interpolation kernel used is

$$\frac{1}{N}\frac{1-e^{j\omega N}}{1-e^{j\omega}} = e^{-\omega N/2}\hbox{asinc}_N(\omega).$$

This causal variation of the ideal spectral interpolation kernel has a linear phase term corresponding to a delay of $N/2$ samples.

Figure F.1: Matlab function specinterp for performing ideal spectral interpolaion.

 

function [si] = specinterp(spec,indices);
%SPECINTERP     
%          [si] = specinterp(spec,indices);
%
%          Returns spectrum sampled at indices (from 1).
%          Assumes a causal rectangular window was used.
%          The spec array should be the WHOLE spectrum,
%          not just 0 to half the sampling rate.
%          An "index" is a "bin number" + 1.
% EXAMPLE:
%          N = 128;
%          x = sin(2*pi*0.1*[0:N-1]);
%          X = fft(x);
%          Xi = specinterp(X,[1:0.5:N]); % upsample by 2
%          xi = ifft(x); % should now see 2x zero-padding

N = length(spec);
ibins = indices - 1;
wki = ibins*2*pi/N;

M = length(indices);
si = zeros(1,M);
wk = [0:N-1]*2*pi/N;
for i=1:M
  index = indices(i);
  ri = round(index);
  if abs(index - ri) < 2*eps
    si(i) = spec(ri);
  else
    w = wki(i);
    % ideal interp kernel:
    wts = (1-exp(j*(wk-w)*N)) ./ (N*(1-exp(j*(wk-w))));
    si(i) = sum(wts .* spec);
  end
end

Note that resampling a spectrum using the ideal interpolation kernel requires ${\cal O}(N^2)$ operations, This computational burden can be reduced in a number of ways:

• A window $w(n)$ other than the rectangular window can be applied to the data frame, so that $W(e^{j\omega})$ can be used as an interpolation kernel in place of the aliased sinc function. If the window sidelobes are sufficiently small, they can be neglected. An example would be $w(n)$ set to a Hann window, with 50% overlap between frames, as commonly used in the STFT.
• In pure frequency-warping applications, windows can be used for both analysis and reconstruction, as is done in apparently all perceptual audio coding methods such as MPEG [15]. For example, the square root of the Hann window can be applied to the analysis frame as well as the final overlap-add frame. This scheme has the advantage of tapering any discontinuities at the frame boundary more gracefully to zero, thus reducing the audibility of artifacts.
• Specially optimized narrow-support kernels can be designed by a number of methods [62,201]. If they are constrained to be nonnegative, their square root can be used twice like the Hann window in the previous point. Asymmetric factorizations are also possible when the kernel changes sign.