OrangeTreeSamples needs harmonic extractor and custom resynthesis

[elanhickler]

SampleTailExtender.zip

Explore the code base, you'll easily find the harmonic analysis and extraction stuff. We need you to repurpose this for extracting individual harmonics.

Right now the goal is to extract harmonics and measure the amplitude decay to then resynthesize it, except that the amplitude decay needs to be smoothed out because there's unwanted harmonic beating or distortion that is causing unwanted modulation.

[elanhickler]

From your commits it looks like you're doing a lot of coding, no idea why you need to do all this extra coding, I thought the SampleTailExtender had most of the code that we needed for phaselocking already done... but you're the expert!

[RobinSchmidt]

no idea why you need to do all this extra coding, I thought the SampleTailExtender had most of the code that we needed for phaselocking already done

hmm...it does some basic amplitude tracking of the partials, so maybe the resulting amplitude tracks can be used to do a fixed frequency resynthesis. but how would you be supposed to obtain the residual, non-harmonic part? for doing that via time domain subtraction, we need to be able to resynthesize the signal completely, with correct time-varying frequencies, amplitudes, and phases. for doing it via frequency domain domain subtraction, we would need spectrogram resynthesis

[RobinSchmidt]

pleeez Robin do you have an ETA? This is the most important thing you will do for me till the end of time!

yea - this is a very research heavy thing once again. but i think, when the basic system is in place, it will allow to build great sample manipulation algorithms on top of it. it's always hard to predict what stumbling blocks will appear down the road. at the moment, it looks like i'm closing in. maybe by end of the week the general system should be up and running - i hope.

[elanhickler]

well I thought it had a basic harmonic removal tool, and it had the ability to analyze phase/frequency/amplitude or at least add/subtract the harmonics. I've seen rudimentary FFT tools do these kinds of things pretty easily.

[RobinSchmidt]

i don't see any of that in the code. can you be more specific where this functionality can be found? maybe you have another version of the code? what sort of tools are you talking about?

ooh---wait....

[elanhickler]

I can't find it but I remember there was a simple puredata patch that could track/isolate harmonics . I remember building Max/MSP patch that could do it... so long ago I don't have any examples. Ok, looking at the code, I can't find a harmonic extraction function, but consider this: In order to paste on harmonics to the end of a guitar samples, you'd have to know the phase of the harmonics around the end of the sample, the amplitude, and the frequency. The tool does a crossfade to go from original to synthesized audio, and it is quite seamless.

[RobinSchmidt]

ok - there are these 3 functions:

    double getMagnitudeForSinusoid (const std::vector<double>& magnitude, int index);
    double getPhaseForSinusoid (const std::vector<double>& phase, int index);
    double getTunedFrequency (const std::vector<double>& magnitude, int index, double sampleRate);

which are called in SampleTailExtender::extendSample like this

    for (auto& h : harmonics)
    {
        double mag = getMagnitudeForSinusoid (magnitudeSpectrum, h);
        double phase = getPhaseForSinusoid (phaseSpectrum, h);
        double frequency = getTunedFrequency (magnitudeSpectrum, h, sampleRate);
     }

but this magnitudeSpectrum, phaseSpectrum is just the spectrum of a single frame - not a spectrogram. it seems to just take a single STFT snapshot at the splicing point from which it then synthesizes the tail. but there's no resynthesis of the section before the tail - that section is just used as is from the original signal

[RobinSchmidt]

In order to paste on harmonics to the end of a guitar samples, you'd have to know the phase of the harmonics around the end of the sample, the amplitude, and the frequency. The tool does a crossfade to go from original to synthesized audio, and it is quite seamless.

yep - that's exactly what it does. and then also apply an amplitude modulation that it seems to get from measuring the amplitude trajectories of the section before the tail. well, it analyzes the amplitude trajectories to extract some high-level parameters (decay-time, beating-freq, beating amount). maybe that amplitude trajectory could be used for a rudimentary sinusoidal resynthesis. but don't expect it to match the original signal to such a degree to allow time-domain subtraction for obtaining a residual. that will almost certainly not work

[elanhickler]

Hmm, well you found the code. I can't offer any more insight.

[RobinSchmidt]

i'm getting close to a perfect resynthesis of a sinusoid - but i'm still having some artifacts. the phase of the resynthesized signal tempoarily desynchronizes from the original phase, leading to bursts in the residual (black: original, blue: resynthesized, green: residual - the fade-in in the resynthesis is due to the abrupt start of the original sinusoid - which may not be an issue in real-world signals). i have some ideas how to fix this phase desync issue and will implement them now. i post this anyway because my ideas may fix this issue for this and many other cases, but maybe not for every possible case - so should such an issue someday re-appear someday with a real-world signal, we know immediately what the problem is and where to look at for fixes:

btw: in this case, the original sound actually contained a 2nd sinusoid at around 10% higher frequency (freqs were around 1 kHz and 1.1 kHz @ 44.1 fs) to make the analysis more challenging - but here i have plotted only one (original and resynthesized) sinusoid

[RobinSchmidt]

i've just listened to the resynthesized sound - this phase anomaly is audible, even though it is really hard to see when just looking at the waveform - is sounds like a sort of soft-reset or something in the sinsusoid - which can perhaps be used for an interesting effect: insert such phase-anomalies at regular intervals and you may get some sort of "robotization" effect (i guess) - and maybe use different intervals/periodicities for the phase-anomaly intervals for different partials. ha! this will be fun to play with!

i mean: i'm not suggesting that this bug is a feature - i mean, i'll fix it but we'll keep in mind that a feature may be constructed from the same idea later - a phase-bumper effect or something

[elanhickler]

We'll probably have to require a high samplerate for low frequency instruments, but it won't be uncommon to have to deal with low frequency instruments where the harmonic spacing is ~50hz

so your 1100hz vs 1000hz shouldn't be too challenging.

[elanhickler]

I thought this stuff was easy because all the FFT programs I've used, if you isolate a band of frequencies in the FFT program and then do a time domain subtraction of those isolated frequencies outside of the FFT program, you get perfect subtraction. This is how denoising works as well. So you could select the harmonics in an FFT program for example, get sinusoids, and then subtract them in the time domain. It's always been easy. Is it just that the programmers behind the FFT programs have solved all things you are working on now?

[RobinSchmidt]

...hmm...well...with more such bumps, it sounds more like phase-modulation, which makes sense and is perhaps not that interesting...we'll see

[RobinSchmidt]

We'll probably have to require a high samplerate for low frequency instruments, but it won't be uncommon to have to deal with low frequency instruments where the harmonic spacing is ~50hz

so your 1100hz vs 1000hz shouldn't be too challenging.

a smaller frequency spacing actually just implies longer time windows, so the absolute frequency spacing is irrelevant - everything will just scale up and down accordingly (freq-resolution goes up -> time resolution goes down, by the same factor). i just wanted to make sure that my analysis can estimate the sine parameters correctly even in the presence of another sinewave (and if that works well - which it does, then probably adding even more sines into the picture won't change the situation very much anymore (i hope))

[RobinSchmidt]

I thought this stuff was easy because all the FFT programs I've used, if you isolate a band of frequencies in the FFT program and then do a time domain subtraction of those isolated frequencies outside of the FFT program

are you talking about spectrogram editors here? ....like, selecting a rectangle in the spectrogram and then synthesizing only the selected area?

[elanhickler]

yes, or iZotopeRX has tools for being more precise than a rectangle. But a rectangle can be sufficient sometimes as well, like a bunch of small rectangles that follow the harmonic.

[RobinSchmidt]

as said - spectrogram resynthesis actually works already (well, not quite exactly the way i want it to, but still). with an appropriate gui, it should be possible to select rectangles and resynthesize only those parts of a spectrogram. but a sinusoidal model is actually a higher level model of abstraction of the sound. the skeleton of the data-structure looks like this:

/** Data structure to hold and edit the instantaneous parameters of a sinusoidal partial at one 
instant of time. */
class rsInstantaneousSineParams
{
  T time = 0;      // in seconds
  T freq = 0;      // in Hz
  T gain = 0;      // as raw factor - rename to amp
  T phase = 0;     // radians in [-pi, pi]
};

/** Data structure to hold and edit the information about one single sinusoidal partial. This is 
effectively an array of rsInstantaneousSineParams objects ordered by time-stamp.  */
class rsSinusoidalPartial  
{
  std::vector<rsInstantaneousSineParams<T>> instParams;
};

/** Data structure to hold the information of a sinusoidal model for a complete sound. This is 
effectively an array of rsSinusoidalPartial objects. */
class rsSinusoidalModel
{
  std::vector<rsSinusoidalPartial<T>> partials;
};

so, as you can see, we literally represent the sound as bunch of frequency, phase and amplitude trajectories in time. the rsInstantaneousSineParams structure is a record of a single datapoint, a partial is a sequence of such datapoints, ordered in time and a complete sound is a set of such partials. now, the task is to write these two classes:

/** Creates a sinusoidal model of an input sound by means of identifying and tracking stable 
sinusoidal partials in its spectrogram. */
class SinusoidalAnalyzer
{
  /** Analyzes the given input sound and returns the sinusoidal model object that models the given 
  sound. */
  rsSinusoidalModel<T> analyze(T* sampleData, int numSamples, T sampleRate);
};

/** A class for synthesizing a sound from a sinusoidal model via interpolating the instantaneous 
phases and amplitudes of the partials up to the sample-rate and driving an oscillator bank with 
these interpolated values.  */
class SinusoidalSynthesizer
{
  /** Synthesizes a sound from a sinusoidal model and returns it as std::vector. The sound will be
  synthesized at a sample-rate that you can set up via setSampleRate. */
  std::vector<T> synthesize(const RAPT::rsSinusoidalModel<T>& model) const;
};

meanwhile, i have solved the phase desync burst problem. i was on the wrong track what caused it two times - implemented a lot of fancy code that actually doesn't solve the problem and even made some things worse - until i found a silly bug in phase-interpolation code. ...that fixed it easily. ..however - i hope this new code may still become useful at some stage - we'll see. i could get a model from two simultaneous sinusoids and resynthesize them very accurately, that the residual is between -80 and -90 dB. i hope to improve this further still - but now another bump appeared: it has problems handling sweeping frequencies. somehow - again - the phase gets messed up. handling the phase correctly is really tricky...

[RobinSchmidt]

did you watch these 3 videos? at the beginning of the 1st, he mentions the STFT - that's just another word for spectrogram: https://www.coursera.org/lecture/audio-signal-processing/sinusoidal-model-1-gjiP7 ...it's a lot of math, but you don't need to deal with these details and may nevertheless get a high level overview

ahh - actually only the first two are relevant - i'm not doing this iFFT synthesis at the moment - i'm just using and oscillator bank for the time being

[elanhickler]

In the case of the SampleTailExtender the "rectangle selection" is done by the harmonic analyzer. So you are doing something higher level, creating a sinewave or single frequency rather than a block of partials via resynthesis.

So is decent resynthesis easier that sinusood models? If so, why don't I use that for harmonic extraction?

[RobinSchmidt]

So is decent resynthesis easier that sinusood models? If so, why don't I use that for harmonic extraction?

well, (re)synthesizing a signal from a spectrogram is easier...or let's say, it's solved - but when you want to start manipulating the spectrogram - like enforcing a particular phase at a particular instant of time for a particular frequency, it's not that straightforward. you cannot just manipulate the phases of the spectrogram at one particular frequency of interest and resynthesize and expect a good result. first, there's no one-to-one correspondence between spectral bins and frequencies - each bin has a certain bandwidth that depends on the window function. ...soo i would say - identity resynthesis from a spectrogram is kinda easy but as soon as you want to apply modifications, there will likely be artifacts

[elanhickler]

when you want to start manipulating the spectrogram - like enforcing a particular phase at a particular instant of time for a particular frequency, it's not that straightforward.

I don't think we need that feature. If resynthesis round trip is solved, then we should be able to extract the bandwidths of frequencies per frame from the original signal. The bandwidths will be small enough to extract harmonics. Resynthesize the original signal without those bins. That will remove the harmonics. Then we synthesize a sinewave per harmonic at a desired frequency and phase while transferring over the amplitude.

[elanhickler]

yay or nay?

[RobinSchmidt]

i'm not sure, if i understand completely what you are suggesting. do you mean, we should forget about time-domain subtraction to obtain the residual and instead do the subtraction in the frequency domain, i.e. just zero out all bins that are associated with the partial's frequency? ...and the advantage would be that we don't need to worry about getting the phase correct? what do you mean by bandwidth of the harmonic? i usually think of the harmonic of having an infinitely narrow (possibly time-varying) line spectrum and the apparent bandwidth in the DFT is only an artifact of the windowing. so what you see is not the spectrum of the signal itself but a convolution of all these infinitely narrow lines with that orange function here (in case of a Hamming window):

so for each partial, we'd get an appropriately scaled and shifted version of that orange thing, all overlapping and summing together. we usually adjust the window length such that the mainlobes of the function just barely don't overlap for a given spacing of the incoming harmonics - but they will typically almost touch - now imagine zeroing out all the mainlobes and think of what remains left then? next to nothing

[RobinSchmidt]

this is how a typical short-time spectrum looks like: do you want to zero out all the lobes and just leave that little stuff in between untouched? or do you want to increase the window size (thereby sacrificing time-resolution), such that the lobes will be much narrower than that?

[elanhickler]

I'd like to try time domain method of subtracting harmonics at some point but for now we should first try frequency domain subtraction as that has proven to work for certain signals with all FFT tools I've tried. It may not work for the lowest frequency signals. That's where time domain subtraction could be crucial.

We want to leave the in-between harmonic content (inharmonic) as untouched as possible. We may need to try different settings such as bin size and window function to see what gives best results.

[elanhickler]

Blackman window function seems ideal.

HAMMING:

HANN:

BLACKMAN:

BLACKMAN with time overlap and smaller bin size:

or there's even narrower functions

[elanhickler]

Here's harmonics removed:

Original:

[RobinSchmidt]

judging from the spectrogram plots, here, the window size must be at least 3 or 4 times as large as would be strictly necessarry for sinusoidal parameter estimation (the empty space between the harmonics is about 3-4 times the width of the harmonics themselves). hmmm...spectral subtraction/zeroing would be another approach. we would need a means to flag each bin as belonging to a harmonic or not (or maybe some continuous measure instead of a boolean flag - something like a mask to multiply the spectrogram with)

i was actually targeting time domain subtraction and was somehow stumped by wanting to refine the frequency estimates by using the phase information (the standard approach is to fit a parabola to the dB-values of 3 bins and solve for the maximum - which works quite well when you look at the numbers, but may give a biased estimate, the error of which may accumulate over time, leading to out-of-phase resynthesis at some point). my idea is that the phase at each frame (and frequency) should be given by (up to some factor) the average frequency in the interval between the two frames (i.e. just their arithmetic mean) multiplied by the time-delta between the frames. because the phase estimates are quite accurate, i could compute refined frequencies from actually measured phases. when doing the math, i arrived at a pentadiagonal system of linear equations - and i don't yet have facilities in my library to solve such systems (i have only code for tridiagonal systems - and dense general systems which is inefficient in this case). so i'm looking into options how to add such a facility..

[RobinSchmidt]

there are lots of linear algebra libraries around, and i'm comparing features, licensing issues etc. ...as well as the option of rolling my own code for this problem (actually maybe a more general case of band-diagonal systems). but i tend to think, in this case, it may not be a good idea to re-invent the wheel and instead add some tried and true linear algebra solutions to the library - which will very probably be useful for many other things down the road. ...there's a translation of linpack for c - but it's lgpl, so it doesn't fit - i'm considering to translate (a subset of) the original fortran sources of linpack or lapack to c++ ...maybe with the help of a fortran-to-c translator (there are at least two of these available)

[elanhickler]

What solution can we start testing sooner rather than later? I have 3 clients all needing my services, which depend on this, and need to get going as soon as possible.

[RobinSchmidt]

i think, i will first try to implement my own pentadiagonal solver based on that paper: https://arxiv.org/pdf/1409.4802.pdf the pseudocode looks straightforward enough to implement. i hope, it will work out (there may be conditions under which the algo may fail - it doesn't seem to use any kind of pivoting, so it may in general encounter divisions by zero - but hopefully that won't occur in this case - we'll see...)

if that fails, we may explore other options. ...at some point, we may want to look into more advanced linear algebra algos/libraries anyway, but maybe that's not yet necessarry for this particular problem

[elanhickler]

...that might as well be japanese... hahaha.

[RobinSchmidt]

the pseudocode looks straightforward enough to implement.

pfff...it was a mess and there were even obvious errors in the paper (the solution vector on rhs of the matrix equation (printed as column) was different from the one printed as row a few lines later) - so who knows, what other - not-so-obvious - errors are still hiding there. i've given up on that and sat down and re-derived the algorithm - just did the Gaussian elimination thing with taking into account the special structure of the matrix and arrived at code that is much neater and more concise - here's a prototype:

std::vector<double> solvePentaDiagnonalSystem(
  std::vector<double>& M, std::vector<double>& L,
  std::vector<double>& D,
  std::vector<double>& U, std::vector<double>& V,
  std::vector<double>& B)
{
  int N = (int) D.size();

  // Gaussian elimination without pivot-search - we just always use D[i] as pivot element:
  int i;
  double k;
  for(i = 0; i < N-2; i++) {
    k       = L[i]/D[i];
    D[i+1] -= k*U[i];
    B[i+1] -= k*B[i];    
    U[i+1] -= k*V[i];
    k       = M[i]/D[i];
    L[i+1] -= k*U[i];
    D[i+2] -= k*V[i];
    B[i+2] -= k*B[i];
  }
  k       = L[i]/D[i];     // a final partial step outside the loop
  D[i+1] -= k*U[i];
  B[i+1] -= k*B[i];

  // Gaussian elimination is done - now do the backsubstitution to find the solution vector:
  std::vector<double> x(N);
  x[N-1] =  B[N-1]                  / D[N-1];
  x[N-2] = (B[N-2] - U[N-2]*x[N-1]) / D[N-2];
  for(i = N-3; i >= 0; i--)
    x[i] = (B[i] - U[i]*x[i+1] - V[i]*x[i+2]) / D[i];

  return x;
}

i tested it and returns a correct result for my test problems. however - i feel not very good about the numerical precision of the algo due to the lack of pivot-search. with a 9x9 matrix/system, i got already errors in the last 3-decimal-digits - and in the frequency estimation problem we are dealing with matrices of sizes of several hundreds. so i'm now indeed looking into porting LAPACK to C++ https://en.wikipedia.org/wiki/LAPACK ...so far, it's going surprisingly smoothly. i already converted a couple of fortran files to c and could build them. i hope, it continues that way. if it does, we'll soon have really advanced matrix algorithms available in RAPT

[elanhickler]

Ok, what is a new ETA for extracting harmonics?

[RobinSchmidt]

puuh - that would have to be a wild guess. i translated the main high-level routine that is relevant for this particular solver and the two subroutines, it calls - now on to their subroutines and then their subsubroutines and so on. it's a bit hard to predict how many branches and twigs this tree will turn out to have. could be just a dozen or a hundred or more (and i guess, it's probably something in between). i'll focus only on the routines that are relevant for this solver that we need for this problem (the generalized band-diagonal solver...well, actually my matrix is symmetric, and i think, there's another, specialized solver for this case - but we may also use the more general one anyway - maybe later the optimized version for the symmetric case)

i started this new repo: https://github.com/RobinSchmidt/LAPackCPP

give me at least a few days for the translations. it's a big pile of grunt work

[RobinSchmidt]

my task just became sooo much easier - there is already a c translation available: https://www.netlib.org/clapack/index.html fortunately, i did not yet translate too much fortran code myself. and having done a bit of it, i have aquired the skill to do so - which may become useful one day because this netlib site is a real treasure trove of math-code: https://www.netlib.org/liblist.html ...but much of it is written in fortran - this language seems to be still a thing in the scientific computing community. how could i not know that website for all these years? :-O i still need to adapt the lapack c-code to c++ to make it work nicely with my template system in rapt - but at least the fortran-to-c translation step is already done...which is very good because for some of the lapack .f files, the f2c translator produced inscrutable error messages.

[RobinSchmidt]

this netlib site is a real treasure trove of math-code...but much of it is written in fortran

for example - this: https://en.wikipedia.org/wiki/SLATEC !!! ...and all in the public domain (not everything on the netlib site though - but a lot)

[RobinSchmidt]

my task just became sooo much easier - there is already a c translation available

hmm.. actually, to be honest, if the f2c translation works well, this additional step is a small one - and the c-version doesn't seem to have the most up-to-date versions of the routines - so i'll probably stick to converting myself and use the pre-translated (older) c-versions only when the converter fails. the text reformatting and templatization is necessarry anyway

[RobinSchmidt]

34 routines translated, so far. i already can use a simple version of the banddiagonal-solver (gbsv) since yesterday - and it returns a correct (but numerically rather imprecise) result. so, at least, i'm doing the right thing. but i need more precision - there are more sophisticated versions of the solver, called gbsvx and gbsvxx - i'm still in the process translating all their lower level dependencies - but i think, i should hit the bottom level soon.

[elanhickler]

Running out of time. Can you just explain how I can zero out some bins? We already have the harmonic detection code based on the simple one found in the tail resynthesis. Even if this isn't an adequate solution, I should still be testing this because I need to get this process started, the tooling, the learning, the infrastructure, the problem solving, etc.

[RobinSchmidt]

if you want to experiment with zeroing out bins in the spectrogram, look at the function spectrogramSine() in the test project. this code creates a sinewave and sends this signal through a spectrogram-analysis-resynthesis roundtrip. in between analysis and resynthesis, you can manipulate the spectrogram (a 2D array, i.e. a matrix - the first index is the frame index and the second index is the bin index). the matrix returned from the spectrogram-analyzer is complex valued, but you can easily convert to magnitude and phase and back (the forward conversion is even done in this test for plotting - the conversion back from magnitude/phase back to real/imag is: real = magnitude cos(phase), imag = magnitude sin(phase) ...or you could just use std::polar) - but if you want to just zero out elements, you can do that directly on the complex valued matrix as well.

[RobinSchmidt]

it's probably obvious, but just in case: if you use zero-padding, the number of bins that you want to zero out per harmonic should be scaled up by the zero-padding factor. be sure to also read the comments in the spectrogramSine function - i just updated them with some more notes what to take care of.

if you want to experiment directly with the code, modify my tests, plot results, conduct your own etc., you may want to use this: https://github.com/RobinSchmidt/GNUPlotCPP/blob/master/Documentation/UserManual.pdf (just follow (some of the) the steps in the "Preliminaries" section - the code for GNUPlotter code is already integrated in the test-codebase - so the last step is already done). with this, you can have an experimentation workflow similar to matlab but directly in c++. this is how i do my experiments mostly

[elanhickler]

I need you to update your UI library to latest JUCE so we can get back and sync so you can look at some UI crashes (due to displaying the parameter value, updating the UI with the value as the parameter is changed, only happens on my plugins, and your library is using async update in a bad way that Jim pointed out, forget the details, we can deal with it once you're ready, it's an easy fix). But before that you need to get the harmonic extraction resynthesis functions working. Don't stop working on this.

[elanhickler]

I'm actually going on a vacation for the next week so I won't have time to even try to do harmonic extraction. Also in a week I will be risking losing work / failing to meet a deadline. If you can get this harmonic extraction stuff done and easy for me to use with a few function calls, then things might work out. Keep going! thank you.

[RobinSchmidt]

i have just finished the translation work for the band-diagonal solvers. boy - was this a big pile of grunt work! i'll file this under preparations/preliminaries so i didn't write work-hours for that. now my test-system gets solved correctly up to the last digit (which the simple solvers get only right up to 8 (decimal) digits or so). perfect. couldn't be better. it's impressive how much algorithmic sophistication goes into making sure full numerical precision in lapack. i'll have to check a couple of things - if i did them really right (i had to make some guesses - it seems to work, but i'll have to double-check), then write a wrapper class for convenient use and then copy the code over to the rapt codebase - ....and then actually invoke it for my frequency-estimation refinement equation - the one with the pentadiagonal system

[RobinSchmidt]

x: original/target x2: simple solver result (gbsv) x4: sophisticated solver result (gbsvxx) there's also a gbsvx with an intermediate level of sophistication - in case you wonder, why there's no x3 shown

[elanhickler]

nice, not sure what this solves, but I see x1 matches x4

[RobinSchmidt]

i have set up an example system of equations A*x = b with known matrix A and vector x and computed the right-hand side b via a matrix-vector multiply, like this:

  // 11 12 13 00 00 00 00 00 00     1     74
  // 21 22 23 24 00 00 00 00 00     2    230
  // 31 32 33 34 35 00 00 00 00     3    505
  // 41 42 43 44 45 46 00 00 00     4    931
  // 00 52 53 54 55 56 57 00 00  *  5 = 1489
  // 00 00 63 64 65 66 67 68 00     6   2179
  // 00 00 00 74 75 76 77 78 79     7   3001
  // 00 00 00 00 85 86 87 88 89     8   3055
  // 00 00 00 00 00 96 97 98 99     9   2930

then, i gave the solver the matrix A and resulting vector b in order to recover the vector x - this is the "solving the system of linear equations" step - the hard step - the thing which all of this is about. the matrix A has zero entries everywhere except on the main diagonal and (in this case) 2 superdiagonals and 3 subdiagonals (this is the special band-diagonal structure which allows for an O(N) algorithm for the solver - instead of the O(N^3) what a regular general-system solver would take - and we also only need O(N) sotrage space for the matrix instead of O(N^2) for general matrices). the black art in all these messy calculations seems to be to ensure that the result has the numerical precision that the machine's number representation actually allows - and not much less (which is what you get when you just naively implement a solver with standard gaussian elimination, say)

[elanhickler]

How many days you think until I have something like this to use?:

/*
audioDataPerHarmonicOut: return audio data for each harmonic, needed for
resynthesis

amplitudeDataPerHarmonicDataOut: return amplitude data for each harmonic, needed 
for resynthesis

noiseOnlyOut: return audio with harmonics removed, needed for resynthesis
*/ 
vector<vector<double>> extractHarmonics(double fundamentalFreq,
 vector<vector<double>> & audioDataPerHarmonicOut
 vector<vector<double>> & amplitudeDataPerHarmonicDataOut
 vector<double> & noiseOnlyOut) 

/*
Other possible parameters, ideally I wouldn't have to set these, but if the algorithm needs help:
    minAmplitude
    maxFrequencyToExtract
*/

Or should we go straight for single function that you create?

/* Returns resnyhesized audio.
Goes through each channel of audio and harmonic in each channel and gets an amplitude
value over time. Rebuilds each harmonic discarding frequency data. Starting phase for
each harmonic is retrieved from phasePerAuioChannelPerHarmonic.
*/
vector<vector<double>> phaseLockAudio(vector<vector<double>> multichannelAudio, vector<vector<double>> phasePerAuioChannelPerHarmonic);