FFT interpolation example

Summary

Introduction

Example application

• The example program FFT_interpolation_example.m
• The function mn_FFT_interpolation.m
• The input audio file in.wav
• Here is also the resulting audio file: out.wav

FFT_interpolation_example.m

close all; clear all;
% Load audio file
[y, fs, bits]=wavread('in.wav');
% nChan is 2 for a stereo file: Independent left / right channel
[n, nChan]=size(y);

yOut=[];
for index=1:nChan
  % extract the channel as row vector
  yc=transpose(y(:, index));
    
  % Interpolate: Lower pitch by one semitone
  yc=mn_FFT_interpolation(yc, 'factor', power(2, 1/12));
  
  % overwrite original data
  % real() is needed to remove imaginary part caused by finite numerical precision
  yOut(:, index)=transpose(real(yc));
end

% In extreme cases, ringing might raise peaks. 
% Therefore, normalize to 1
peak=max(max(abs(yOut)));
yOut=yOut*1/peak;

% write output file
wavwrite('out.wav', yOut, fs, bits);

mn_FFT_interpolation.m

function data=mn_FFT_interpolation(data, varargin)

% ****************************************************************
% param check I
% ****************************************************************
if mod(length(varargin), 2) > 0
  error('wrong number of arguments (need option, value pairs)');
end
[nRow, nCol]=size(data);
if nRow ~= 1 
  error('data must be row vector');
end

% ****************************************************************
% default values
% ****************************************************************
nPad=-1; % invalid

% ****************************************************************
% parameter handling
% ****************************************************************
for argindex=1:2:length(varargin)
  arg=varargin{argindex};
  param=varargin{argindex+1};
  switch(arg)
   case 'factor'
    nPad=floor((param-1)*nCol+0.5);
   case 'nSamplesAdded'
    nPad=nCol+param;
   case 'nSamplesTotal'
    nPad=param-nCol;
   otherwise
    error(['invalid parameter at position ', mat2str(argindex)]);
  end % switch
end % for argindex

% ****************************************************************
% param check II
% ****************************************************************
if nPad < 0
  error('invalid / incomplete parameters');
end

if nPad == 0
  return % nothing to do
end
  
realValued=isreal(data);

% ****************************************************************
% to frequency domain
% ****************************************************************
data=fft(data);

% ****************************************************************
% Zero-pad frequency domain data
% ****************************************************************
if mod(nCol, 2)==1
  % ****************************************************************
  % odd length
  % ****************************************************************
  lside=data(1:(nCol+1)/2); 
  rside=data((nCol+1)/2+1:nCol);
  data=[...
      lside, ... % positive frequencies
      zeros(1, nPad), ... % zero padding
      rside ... % negative frequencies
       ];
else
  % ****************************************************************
  % even length
  % The center bin needs to be split in two 
  % since it contains sum of exp(-i pi k) and exp(i pi k) terms
  % ****************************************************************  
  lside=data(1:nCol/2); 
  center=data(nCol/2+1);
  rside=data(nCol/2+2:nCol);
  data=[...
      lside, ... % positive frequencies
      center/2, ... % previous Nyquist frequency (positive half),
      zeros(1, nPad-1), ... % zero padding
      center/2, ... % previous Nyquist frequency (negative half),
      rside ... %negative frequencies
      ];
end

% ****************************************************************
% to time domain
% scale, so that average amplitude remains as it was
% ****************************************************************
f=(nPad+nCol)/nCol;
data=ifft(f*data);

% ****************************************************************
% Chop off imaginary part, if the input data was real-valued
% ****************************************************************
if realValued > 0
  data=real(data);
end

Comments? Questions?