Curve fitting a frequency sweep (user defined function)

Share icon

I am attempting to develop a curve fitting function for a frequency sweep (a sin wave with a time varying phase function, logarithmic in this case).

the equation i use is the following:

f0= ln(K0)
f1= (ln(k1)-f0)/StimDur
f2=2*pi*e^f0/f1
f3= mod(f2,2*pi)
f(x) = Amp*sin(f2*e^(f1*x)- f3)

this equation was actually taken from one of the sample igor files.

Fit coefficients:
K0 is the start frequency (Hz)
K1 is the end frequency (Hz)
StimDur is the stimulus duration (s)

I believe F3 is there to ensure that the sin wave always begins at zero and goes positive.

i can successfully generate a wave with this function and use the same function to get a good fit to the generated wave.

my problem arises when i try to fit the function to some real data (acquired at the output of an amplifier, see attached file).

the fit i obtain does not even come close despite the fact that it only takes a couple of iterations to reach the cut off chi squared delta.

my goal is to remove/notch out/subtract the fundamental from the amplifier output and obtain a total harmonic distortion measurment (THD). Since THD measurements don't include phase errors the curve fit needs to be able handle the linear phase shifts that amplifiers produce (on the order of a few degrees). The frequencies that the are typically used are 20hz to 20khz.

thank you in advance for your help. trialF_test_2.ibw (885.35 KB)
I would suggest trying to construct a smooth time-varying phase function using the zero-crossings of data (with mean or baseline subtracted). These should be easy to find using a user-defined function. Then use the zero-crossings to implement non-linear time scaling, such that the time-scaling forces the data to approximate a single sinusoid on a y vs scaled time (not uniformly spaced points) axis. Then interpolate the y-data to a uniformly-spaced time scale, and do standard spectral analysis. THD should appear as sidebands or harmonics, in addition to the dominant spectral peak.

Stephen R. Chinn

January 29, 2009 at 11:05 am - Permalink

Your model for the phase is going to have to match the actual data really closely. Even a small mismatch will result in the phase being alternately correct and a 180 degrees out of phase. When it's 180 out of phase, the amplitudes of the data and the curve fit will be opposed. My guess is that the fit is just finding a local minimum where the amplitude is about right.

Steve is a pretty smart guy. You might try his solution. I think fitting this thing will be very difficult.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

January 29, 2009 at 12:47 pm - Permalink

thanks for your feedback. it took me a while to implement what stephen suggested and now i have a couple of follow up questions. First i'll give a summary of what i did.

1. using find_levels i found all the zero crossings. from this i determined the time between crossings, took the inverse and plotted it against the original time scale.

2. i then fit the curve with a sigmoid- it provided the best fit, but was still out by several hz in some locations on the curve

3. then, using the coefficients of the curve fit i applied the sigmoid equation to the original evenly spaced time data (in this case the sample interval was 1/44100).

4. i then plotted the original frequency sweep against the new time data. this gave a wave with a single fundamental frequency (i had some trouble figuring out just what frequency it was)

5. i used the igor interpolate function to evenly space the time data. i chose a delta that maintained the highest sampling rate of the pre-interpolated data. i used linear interpolation.

6. i ran a FFT on the data and plotted the first column of the magnitude output wave.

i did this on real data as well as on a ideal frequency sweep (generated with the equation in the first post).

In both cases i did not obtain any sidebands. rather i had the fundamental with a relatively high noise floor, only 50 or 60dB. i can't really explain why the noise floor would be so high on the "ideal frequency sweep". is it because my sampling rate is too low for the higher frequencies (sampling at 44.1 khz and the sweep goes to 20 khz) or is it something else.

also, if the curve fit is not perfect would this widen the fundamental or show up in some other way?

lastly i thought of another way to fit the frequency sweep - what are your thoughts on fitting each cycle of the sweep with the sin function, then concatenating each fit to generate a fit for the entire frequency sweep? i realize this is a very ugly and slow way of doing it but do you think the end result would be usable?

thanks again for your help - I've already learned a lot.

February 9, 2009 at 07:57 pm - Permalink

I can't quite follow what you have done. I don't think any curve fitting should be necessary. Here is a bare bones example of a procedure that illustrates my proposed method for a simple linear frequency chirp:
(construct your own graphs, where indicated)
function doit()
    make/O/N=800 wphase, wphase2, wtest
    setscale x, 0, 1,"s", wphase, wphase2, wtest
    wphase  =    2*pi*(20*x + 40*x^2)          //   frequency is time ('x') derivative of phase; linear chirp
    wphase2 =    2*pi*(40*x + 80*x^2)          //  weak second harmonic added
    wtest   =   sin(wphase) + 1e-2 * cos(wphase2)   // construct simulated test data;  DISPLAY THIS
       
    FindLevels/Q  wtest, 0   // find zero crossings;
    WAVE w_FindLevels   //  might also use extrema
   
    variable Npts = numpnts(w_FindLevels) //    create a linear time wave whose
    make/O/N=(Npts) wxlin   //  initial frequency is close to the data's
    setscale/P x, w_Findlevels[0], w_FindLevels[1]-w_FindLevels[0],"s", wxlin
    wxlin   =   x   // DISPLAY W_FindLevels vs wxlin
   
    Make/O/N=1024 wdechirp, xinterp  //  create a dechirped wave from the data; use a number of points
    setscale x, wxlin[0], wxlin[Npts-1], "s" wdechirp, xinterp    // convenient for FFT
    xinterp = interp(x, wxlin, w_FindLevels)    // these steps are the key; wdechirp has a linear time scale,
    wdechirp    =   wtest(xinterp)  // but its y-values come from the 'data' evaluated from
// the non- linearly scaled times;  DISPLAY THIS; it shows a quasi single-frequency wave
    FFT/DEST=FFTout/OUT=4/WINF =hanning wdechirp    //  find its power spectra; DISPLAY FFTout
end


Try this yourself. Please forgive any screwed up formatting; I'm new to HTML tags.

Stephen R. Chinn

February 10, 2009 at 10:13 am - Permalink

I have just tried my time-axis scaling code on the previous binary wave you had posted with logarithmic chirp. I eliminated a small subset of problematic points at the beginning and end of the wave, retaining 400,000 points. I increased the scaled, interpolated wave and FFT sizes to 16,384 and did no windowing. I found an essentially single frequency spectrum, with greater than 63dB signal-to-baseline. For general data, without special circumstances on frequencies and number of cycles, I think you will be hard pressed to improve on this FFT SNR. I noticed your data type was INT16, and that may have had some effect.

If this SNR is not sufficient for you, another method will be necessary.

Stephen R. Chinn

February 11, 2009 at 11:29 am - Permalink