//
//  fft3.sci
//
//  ------------------------------------------------------------------------
//  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
//  ------------------------------------------------------------------------



//
//  fft3.sci
//
Fs = 44100;
A = read('../resample/asn5a.dat', -1, 2);
tx = A(:,1);
y = A(:,2);
m = max(size(tx));
tn = 1:m;

//
// section 1  -- a cur with constant amplitude
//
n1 = 31039;
n2 = 37708;
q1 = n1:n2;
m1 = max(size(q1));  // size of the cut
y1 = y(q1);          // y data from the cut
tx1 = tx(q1);
tn1 = 1:m1;
fy1 = fft(y1, -1);  // FFT of the cut data
T1 = m1 / Fs;
ff1 = 1 / T1;       // frequency increments in FFT

//
// zero-out the noise
//
ns = find(abs(fy1) < 5);
zz1 = zeros(fy1);
fyy1 = fy1;
fyy1(ns) = zz1(ns);
fyy1(1) = 0 + %i * 0;   // zero out DC
fyy1 = fyy1 * 10;   // increase SNR
//
subzeros = zeros(1,40) + (%i * zeros(1,40)); // remove looping noise
fyy1(1:40) = subzeros';
//
//  time stretch by scaling the indices
//
M2 = m1 * 10;               // size of stretched data
fyyy1 = zeros(1,M2)  +  %i * zeros(1,M2);  // placeholder for larger spectrum
nq = find(abs(fyy1) > 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 * 10)+1;     // multiply the factor and shift DC back
T2 = M2 / Fs;
ff2 = 1 / T2;       // frequency increments in stretched FFT
//
// copy the (nonzero) frequencies
//
for i = 1:snq2,
    //
    //  first half
    //
    fyyy1(nqq(i))  =  fyy1(nq(i));
    //
    //  2nd half (mirror image of first half)
    //
    re1 = real(fyy1(nq(i)));
    im1 = imag(fyy1(nq(i)));
    fyyy1(M2 - nqq(i) + 2)  =  re1 - %i * im1;
end


// inverse fft
yyy1 = fft(fyyy1, 1);
yyy1 = real(yyy1);

wavwrite(yyy1, Fs, 'asn5a_stretch.wav');