function y = fft_stretch(data, ifactor, nsthresh, pr)
//
//  y = fft_stretch(data, ifactor)
//
//  ------------------------------------------------------------------------
//  copyright 2007  Alfred Steffens Jr.
//  -----------------------------------------------------------------------
//  -----------------------------------------------------------------------
//  Copying Permission:
//
//    This is free software; you can redistribute it and/or modify
//    it under the terms of the GNU General Public License as published by
//    the Free Software Foundation.
//
//    This software is distributed in the hope that it will be useful,
//    but WITHOUT ANY WARRANTY; without even the implied warranty of
//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//    GNU General Public License for more details.
//
//    You should have received a copy of the GNU General Public License
//    along with this software (see file called "COPYING"); if not, write to
//    the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
//    MA 02111-1307  USA
//  ------------------------------------------------------------------------
//
//
//  A utility function used in code for time-stretching a sound array.
//  1.  "data" is an row vector of real data that the calling code has cut from
//      the source sound data.
//  2.  The time is stretched by increasing the number of sample points by
//      an integral factor, "ifactor".
//  3.  "nsthresh" is the maximum (fractional) absolute amplitude of
//      frequencies to keep.
//  4.  "pr" is a parity adjustment depending on where there are an odd or even
//      number of samples.  Set this to 0, -1, or 1 to get rid of modulation
//      resulting from misplacement of the upper frequency bins.  When in
//      doubt, try a value of 1.
//
//  Returns y = the stretched sound data.
//
//  -----------------------------------------------
//  copyright Alfred Steffens Jr., 2007
//  GNU public license
//  -----------------------------------------------
//
//

m1 = max(size(data));
//
fydata = fft(data, -1);
rmax = max(real(fydata));
imax = max(imag(fydata));
amax = max([rmax, imax]);
ns = find((abs(fydata)/amax) < nsthresh);
zz1 = zeros(fydata);
dcsave = fydata(1);
fydata(ns) = zz1(ns) +  %i * zz1(ns);   // zero the noise
fydata(1) = 0  +  %i * 0;               // clear the DC value
//
//  scale the time
//
M2 = m1 * ifactor;      // new array size
fyydata = zeros(1,M2)  +  %i * zeros(1,M2);     // placeholder for new freqs
nq = find(abs(fydata) > 0);     // indices of the nonzero data
//
//  take the first 1/2 of indices
//
snq = max(size(nq));
snq2 = int(snq/2);
//
// shift over the dc component before multiplying
//
nq2 = nq(1:snq2) - 1;
nqq = (nq2 * ifactor)+1;     // multiply the factor and shift DC back
//
// copy the (nonzero) frequencies
//
fyydata(1) = dcsave;
for i = 1:snq2,
    //
    //  first half
    //
    fyydata(nqq(i))  =  fydata(nq(i));
    //
    //  2nd half (mirror image of first half)
    //
    re1 = real(fydata(nq(i)));
    im1 = imag(fydata(nq(i)));
    fyydata(M2 - nqq(i) + 1 + pr)  =  re1 - %i * im1;
end
//
//
// inverse fft
//
yyy1 = fft(fyydata, 1);
//
y = real(yyy1);