Non-linear curve fitting for forced vibration acceleration data

I'm trying to fit a curve to my data using MATLAB's lsqcurvefit function. The data is from an accelerometer (in g's), and I'm trying to fit curves at each of the resonance modes.

The function that the non-linear solver is using is the standard one for amplitude as a function of frequency:

I also measured the force, which had a max value of around 3.5 N, and mass was around 3-4 kg of the whole system. I applied a FFT to the data and added a Hanning window, so now I essentially have amplitude v.s. frequency (Hz). When I try giving lsqcurvefit initial conditions which are close to what I would expect to final values, it gets it completely wrong (see below). I'm also only displaying the datapoints close to the peak for brevity.

On the other hand, when I played around with the initial estimate, I was able to get an ok fit but with completely arbitrary values:

Below is my code:

% Using equation x = F./(m*sqrt(((2*omega_0*damping).^2)+((omega_0.^2 - omega.^2).^2)));
% x(1) = omega_0;
% x(2) = damping;
% x(3) = F;
% x(4) = m;
pa_avg = [0.0034;0.0033;0.0036;0.0043;0.0071;0.0113;0.0130;0.0117;0.0114;0.0123;0.0130;0.0168...
        0.0228;0.0218;0.0147;0.0122;0.0173;0.0257;0.0321;0.0321;0.0320;0.0406;0.0533;0.0629...
        0.0714;0.0787;0.0763;0.0619;0.0429;0.0259;0.0165;0.0145;0.0144;0.0134;0.0114;0.0088...
        0.0070;0.0057;0.0042;0.0028;0.0021;0.0021;0.0023;0.0028;0.0030;0.0024;0.0018;0.0015...
        0.0015;0.0015;0.0016;0.0014;0.0012;0.0010;0.0010;0.0009;0.0009;0.0009;0.0009;0.0009;0.0008];
f = 890:10:1490; 
w = 2*pi*f(90:150); % Convert to radians/sec
func = @(x,xdata) x(3)./(x(4).*sqrt(((2*x(1)*x(2)).^2)+((x(1).^2 - xdata.^2).^2)));
x0 = [2*pi*1140, 0.3, 10000, 50]; % Initial guess which results in an ok fit
x = lsqcurvefit(func,x0,w,pa_avg(90:150));
semilogy(w/(2*pi),func(x,w))
hold on
semilogy(w/(2*pi), pa_avg(90:150), "o")
legend("Fitted function", "Actual Data")

I tried integrating the acceleration data twice to get displacement but this didn't work at all.

I also tried treating F_0/m as a single variable, acceleration, but this didn't make any difference. Any help is much appreciated.

Answer 1

[Star Strider]

There are some errors.

First, ‘w’ should be calculated by multiplying it by , or using the deg2rad function. It also needs to be converted to a column vector to match ‘pa_avg’ that also has problems, and should be:

pa_avg = [0.0034;0.0033;0.0036;0.0043;0.0071;0.0113;0.0130;0.0117;0.0114;0.0123;0.0130;0.0168;...
        0.0228;0.0218;0.0147;0.0122;0.0173;0.0257;0.0321;0.0321;0.0320;0.0406;0.0533;0.0629;...
        0.0714;0.0787;0.0763;0.0619;0.0429;0.0259;0.0165;0.0145;0.0144;0.0134;0.0114;0.0088;...
        0.0070;0.0057;0.0042;0.0028;0.0021;0.0021;0.0023;0.0028;0.0030;0.0024;0.0018;0.0015;...
        0.0015;0.0015;0.0016;0.0014;0.0012;0.0010;0.0010;0.0009;0.0009;0.0009;0.0009;0.0009;0.0008];

I used the ga (genetic algorithm) function and quickly got a reasonably decent fit (at least by my criteria) with:

x =
          19.8404873275757
              0.3168359375
          35.9852020359039
          34.2824745559692

The fitness function I used:

          ftns = @(x) norm(pa_avg - func(x,w));

producing:

Note that ‘x(1)’ is in rad/sec, so it would convert to 1136.76942273353 °/sec.

.

[Matija Milenovic]

Hi Star Strider, thanks for your answer.

I was wondering why you're using degrees in the solution? My frequency measurements are in cycles per second (Hz), so shouldn't the conversions be done by multiplying/dividing by 2*pi?

Otherwise, your result of 1136 degrees per second matches my 1136 Hz, so I'm guessing both are two ways of arriving at the same solution. I also realise that your x(3) and x(4) are approximately "g" times larger than the real values, which is good.

I will try using the genetic algorithm since it seems lsqcurvefit can't handle the non-linearity. Cheers!

[Star Strider]

My pleasure.

I was not aware of what you are doing.

Using:

w = 2*pi*f(:);% Convert to radians/sec

the estimated parameters are:

x =
                  7143.125
                  111.6795
                  2678.245
                      0.02

producing:

These are constrained only to be greater than 0.

Different runs produced:

x =
                7144.40625
                105.136125
                 8447.9525
        0.0662499999999999

and:

x =
                 7144.1875
                   105.057
                  8929.715
                      0.07

with essentially the same result.

.

[Matija Milenovic]

I've played around with this, and it seems that your initial answer was probably the best! If it was an accident that you converted to degrees, then it's a really great accident to happen!

Further, I also agree with the use of ga. Getting similar results with lsqcurvefit is virtually impossible, at least for this particular function. I could only get converged values when I plugged something very close to your

x =
          19.8404873275757
              0.3168359375
          35.9852020359039
          34.2824745559692

The values in your second answer are close to my initial results, which unfortunately aren't close to the real-life ones.

However, many thanks again!

[Star Strider]

As always, my pleasure!