// // ffttrack.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 // ------------------------------------------------------------------------ // // ffttrack.sci // // First step to create a frequency track for additive synthesis // // 1. Reads the source data file containing the entire note as an ascii // data file. // // 2. Reads a data file that is a single column of sample numbers, between // which the FFT's will be calculated (sampletrack.dat). // // 3. Performs an FFT on each section of the data between the sample // numbers. // // 4. Creates a matrix of frequencies in which each row is the FFT of // data between the sample points. The first column of the matrix is // the times corresponding to each FFT. Then the columns are interleaved // frequency, real amplitude, and imaginary amplitude. // // 5. Writes the matrix of frequencies to a file called freqtrack.dat. // // Usage: // // 1. Modify the line, // A = read('../resample/asn5a.dat', -1, 2); // to read your input data file (the sound file as an ASCII file). // 2. Put your sample number boundaries in a file called sampletrack.dat. // 3. Change the sample rate, Fs, if needed. // 4. Change the maximum number of frequencies, Fqz, to keep, if needed. // 5. Change the noise threshold, threshold, to optimize the elimination // of sidebands around the resonance lines. // 6. Change maxfreq if needed. // 7. At the Scilab prompt, type // exec('ffttrack.sci'); // // // // the sample rate // Fs = 44100; // // number of frequencies to track over time // Fqz = 50; // number of frequencies to keep threshold = 0.04; // smallest relative amplitude to keep maxfreq = 7000; // highest absolute frequency to keep smallblock = 1; debugger = 0; // set to zero to turn off, 1 to turn on // // the input sound data already converted into ascii numbers // A = read('../resample/asn5a.dat', -1, 2); tx = A(:,1); y = A(:,2); m = max(size(tx)); tn = 1:m; // // the data file, a single column, of sample numbers marking the sections // B = read('sampletrack.dat', -1, 1); sampmax = max(B); dataF = []; dataD = []; // DEBUGGER dataD2 = []; // DEBUGGER row = 0; // // number of time frames // tmframes = max(size(B)) - 1; fnummax = 0; frowmax = 0; gnummax = 0; growmax = 0; // // loop through each time section // for k = 1 : tmframes, // // get the current section from the data // n1 = B(k); m1 = B(k+1); q1 = n1:m1; M1 = max(size(q1)); // size of the cut y1 = y(q1); // y data from the cut tx1 = tx(q1); // time axis tn1 = 1:M1; // sample num axis fy1 = fft(y1, -1); // FFT of the cut data ff1 = Fs / M1; // frequency increments in FFT fm1 = (tn1-1)*ff1; // frequency axis without DC fn1 = [0, fm1]; // frequency axis complete fn1 = fn1(1:M1); // frequencies // throw away DC fy1(1) = 0 + %i*0; // find the loudest frequency mxx = max(abs(fy1)); // // normalize the frequencies (divide by the loudest) // real_array = [real(fy1)/mxx]'; imag_array = [imag(fy1)/mxx]'; norm_array = [fy1 / mxx]'; // // // zero-out the noise // // // indices in the array of those amplitudes greater than the threshold // ind = find(abs(norm_array) > threshold); // // get how many indices were found // indsize = max(size(ind)); // // if we found more than the requested length, truncate the indices // if (indsize > Fqz) then ind2 = ind(1:Fqz); else ind2 = ind; end // // the number of indices could actually be less than Fqz (but not greater) // ind2size = max(size(ind2)); // // use only the first 'Fqz' number of valid frequencies // rarray = zeros(1,Fqz); iarray = zeros(1,Fqz); farray = zeros(1,Fqz); rarray(1:ind2size) = real_array(ind2); iarray(1:ind2size) = imag_array(ind2); farray(1:ind2size) = [fm1(ind2)]; if (debugger == 1) then dataD = [dataD ; farray]; end // // save only those frequencies in the lower half of the fft // ind3 = []; for qk = 1:ind2size, if (farray(qk) <= maxfreq) then ind3 = [ind3, qk]; else break; end end ksz = max(size(ind3)); farray = [farray(ind3), zeros(1,Fqz-ksz)]; rarray = [rarray(ind3), zeros(1,Fqz-ksz)]; iarray = [iarray(ind3), zeros(1,Fqz-ksz)]; if (debugger == 1) then dataD2 = [dataD2 ; farray]; end // // look for missing harmonics (zeros in between the nonzero freqs) // checking for nonzeros is ok because the arrays are packed toward // the left // ind4 = find(farray > 0); // these numbers not always contiguous i4m = max(size(ind4)); // size of ind4 if (i4m > fnummax) then fnummax = i4m; frowmax = k; save_farray = farray; end ind4max = max(ind4); // index of highest frequency if (ind4max > gnummax) then gnummax = ind4max; growmax = k; end n4 = 1:ind4max; // count up to the highest frequency tmparray = zeros(1, Fqz); tmparray(n4) = farray(n4); farray = tmparray; tmparray = zeros(1, Fqz); tmparray(n4) = rarray(n4); rarray = tmparray; tmparray = zeros(1, Fqz); tmparray(n4) = iarray(n4); iarray = tmparray; Fqx = Fqz; // // new data row in the output matrix // row = row+1; // // create a new row of zeros // dataF = [dataF; zeros(1, (Fqx*3)+1)]; // // indices for interleaving the data // Fqzz = 1:3:(Fqx*3); // // starting time for this frequency change // dataF(row, 1) = n1 / Fs; // // interleave the frequency values // dataF(row, Fqzz+1) = farray; // // interleave the real amplitues corresponding to these frequencies // dataF(row, Fqzz+2) = rarray; // // interleave the imag amplitues corresponding to these frequencies // dataF(row, Fqzz+3) = iarray; end // // save the data as an ASCII file // since this is normalized data, we don't need "engineering notation" // anything smaller than 0.000100 is noise anyway // // the format of the output data is // 1. Each row is a different time in the life of the note // 2. The first number in each row is the time in seconds since the // beginning of the note // 3. The 2nd through Nth column is the interleaved freuqency components // in groups of 3 numbers: // --- the frequency in cycles per second // --- the amplitude (-1 -- +1) of the real component // --- the amplitude (-1 -- +1) of the imag component // // example: // time freq1 real1 imag1 freq2 real2 imag2 freq3 real3 imag3 ...... // time freq1 real1 imag1 freq2 real2 imag2 freq3 real3 imag3 ...... // . . . . . . . . . . // . . . . . . . . . . // // if (smallblock == 1) then mprintf('fnummax= %d at row %d : gnummax= %d at row %d', fnummax, frowmax, gnummax, growmax); NewCols = (gnummax * 3) + 1; dataF = dataF(:,1:NewCols); end fprintfMat('freqtrack.dat', dataF, '%.6f'); if (debugger == 1) then fprintfMat('debugger.dat', dataD, '%.6f'); fprintfMat('debugger2.dat', dataD2, '%.6f'); end mprintf('file freqtrack.dat has %d columns', NewCols);