frequency equation for curve fitting

Hello All,

Is there any equation frequency based for exponential decaying cosine signals to fit on measure data for optimization or curve fitting method?

Thank you

Answer 1

[Star Strider]

To do a nonlinear fit in the frequency domain, you will need to calculate a function that is the Fourier transform of the time-domain function you want to fit.

Try this:

syms k1 k2 t Tmax w w0 
y(t,w0,k1,k2) = exp(k1*t) * cos(k2*w0*t)
Y = int(y*exp(1j*w*t), t, 0, Tmax)
Y = simplify(Y, 'Steps',250)
Yf = matlabFunction(Y)                                      % Individual Parameters
Yf = matlabFunction(Y, 'Vars',{[w0,k,Tmax],w})              % Parameter Vector ‘in1’

producing:

Yf = @(in1,w) -(in1(:,2)+w.*1i)./((in1(:,2)+w.*1i).^2+in1(:,1).^2)+(exp(in1(:,3).*in1(:,2)+in1(:,3).*w.*1i).*(in1(:,1).*sin(in1(:,3).*in1(:,1))+cos(in1(:,3).*in1(:,1)).*(in1(:,2)+w.*1i)))./((in1(:,2)+w.*1i).^2+in1(:,1).^2);

where ‘k’ (‘in(:,2)’) is the exponential decay constant, ‘w0’ (‘in(:,1)’) is the frequency of the cosine function, and ‘Tmax’ (‘in(:,3)’)is the maximum (end) time of the vector, and ‘w’ is the independent variable and radian frequency. Note that this produces a complex result, so it might be easiest to take the abs() of it and use it with the abs() of the Fourier transform. Also, as written here, ‘in’ is defined as a row vector, and will throw an error if you use a column vector as your initial parmeter estimates.

I have never done any nonlinear parameter estimation in the frequency domain, so I have no experience with it. I also did not test this function with the Fourier transform of a decaying cosine function, so I have no idea if it will succeed.

Have fun!

[Ill ch]

Dear Star Strider,

Thank you so much for your spontaneous reply. Actually there are two problem here:

Is it possible to use this function without syms as i dont have symbolic math toolbox and second thing is it possible to use direct frequency domain equation rather transformation from time domain?

Thank you again for your reply. In time domain my equation looks like this $X(t)=e(-d*t)*cos(wt) where A amplitude and another similar to your equation.

Your suggestions will be helpful for me

[Star Strider]

My pleasure.

You do not need to use the Symbolic Math Toolbox. I posted the ‘Yf’ anonymous function. That is all you need for the nonlinear parameter estimation. It should work with every nonlinear parameter estimation function available in MATLAB.

‘is it possible to use direct frequency domain equation rather transformation from time domain?’

The ‘Yf’ function is the frequency-domain function, created by taking the Fourier transform of ‘y’:

syms A delta phi t Tmax w w0 
y(t,w0,A,delta,phi) = A*exp(-delta*t) * cos(w0*t-phi)
Y = int(y*exp(1j*w*t), t, 0, Tmax)  

It is necessary to begin with the time-domain function in order to create the frequency-domain function.

‘In time domain my equation looks like this ... where A amplitude and another similar to your equation.’

It would have been very helpful to have known that at the outset! The new matlabFunction call is:

Yf = matlabFunction(Y, 'Vars',{[w0,A,delta,phi,Tmax],w})    % Parameter Vector ‘in1’

so the ‘in1’ vector is in order: [w0,A,delta,phi,Tmax], so ‘in1(:,1)’ is ‘W0’ and ‘in1(:,5)’ is ‘Tmax’.

The new ‘Yf’ function is:

Yf = @(in1,w) -(in1(:,2).*(in1(:,3).*cos(in1(:,4))-w.*cos(in1(:,4)).*1i+in1(:,1).*sin(in1(:,4))))./(in1(:,3).*w.*2.0i-in1(:,3).^2+w.^2-in1(:,1).^2)+(in1(:,2).*exp(-in1(:,5).*in1(:,3)).*(cos(in1(:,5).*w)+sin(in1(:,5).*w).*1i).*(in1(:,3).*cos(in1(:,4)-in1(:,5).*in1(:,1))-w.*cos(in1(:,4)-in1(:,5).*in1(:,1)).*1i+in1(:,1).*sin(in1(:,4)-in1(:,5).*in1(:,1))))./(in1(:,3).*w.*2.0i-in1(:,3).^2+w.^2-in1(:,1).^2);

Nonlinear parameter estimation depends on appropriate initial parameter estimates, so this usually requires some experimentation.

If you would prefer to do a time-domain parameter estimation (that would likely be much more reliable), consider the approach in: How to filter noise from time-frequency data and find natural frequency of a cantilever? That code also calculates good initial parameter estimates for the nonlinear regression.

[Star Strider]

You most likely need to choose different initial parameter estimates.

I would have to see your code and your data. I cannot help you without them.

[Ill ch]

Dear Star Strider,

Thank you very much again

here attached is my code what i have tried.

Thank you Very much again for your help

[Star Strider]

My pleasure.

I was able to get a reasonable but obviously not a perfect fit because your data and the function are not entirely the same.

My code:

Measured_signal = readmatrix('MyArtificial.xlsx');
Fs = 12000;                                                     % sampling frequency
L = numel(Measured_signal);

t = linspace(0, 1, numel(Measured_signal))/Fs;

figure
plot(t, Measured_signal)
grid

Fn = Fs/2;
FT_s = fft(Measured_signal)/L;
Fv = linspace(0, 1, fix(L/2)+1)*Fn;
Iv = 1:numel(Fv);

figure
plot(Fv, abs(FT_s(Iv))*2)
grid

w = 2*pi*Fv';
y = abs(FT_s(Iv))*2;

fit = @(in1,w)-(in1(:,2).*(in1(:,3).*cos(in1(:,4))-w.*cos(in1(:,4)).*1i+in1(:,1).*sin(in1(:,4))))./(in1(:,3).*w.*2.0i-in1(:,3).^2+w.^2-in1(:,1).^2)+(in1(:,2).*exp(-in1(:,5).*in1(:,3)).*(cos(in1(:,5).*w)+sin(in1(:,5).*w).*1i).*(in1(:,3).*cos(in1(:,4)-in1(:,5).*in1(:,1))-w.*cos(in1(:,4)-in1(:,5).*in1(:,1)).*1i+in1(:,1).*sin(in1(:,4)-in1(:,5).*in1(:,1))))./(in1(:,3).*w.*2.0i-in1(:,3).^2+w.^2-in1(:,1).^2); %b(1)*exp(-b(2).*w).*cos(b(3).*w-b(4));%.*(cos(2*pi*x./b(2) + 2*pi/b(3))) .* exp(b(4).*x) + b(5);     % Function to fit
fcn = @(in1) sum((abs(fit(in1,w)) - y).^2);                                              % Sum-squared-error cost function
[s, fval] = fminsearch(fcn, [750,  1,  200,  0, 0.5 ]);

figure
plot(Fv, abs(FT_s(Iv))*2)
hold on
plot(Fv, abs(fit(s,w)))
hold off
grid

fprintf(1, 'W_0\t\t= %12.6f\nA\t\t= %12.6f\ndelta\t= %12.6f\nphi\t\t= %12.6f\nTmax\t= %12.6f\n', s)

producing these parameters:

W_0		=  1444.369504
A		=     0.012463
delta	        =   -98.871919
phi		=     0.056982
Tmax	        =     0.186358

and this plot:

That is an acceptable fit, considering that your data are not actually an exponentially-decaying cosine.

[Ill ch]

Good Morning Star Strider

Thank you very much.It works and too great codes. I have one last question if i want to get from this result to againtime based signal then can i do simple ifft like

New_Time_Based=ifft(fit(s,w),L,'symmetric') ? or New_Time_Based=ifft(abs(fit(s,w)),L,'symmetric') ?

as my purpose is to generate signale based on frequency domain to time domain. Now i have frequency domain signal now this signal i have to convert into timedomain so that should be same like original signal

Thank you and Best Regards

Ill

[Star Strider]

As always, my pleasure.

This is probably closest to being correct:

New_Time_Based=ifft(fit(s,w),L,'symmetric')

however it would likely be best to evaluate it on the interval of (-w,+w) to give it some symmetry before you do the ifft.

[Ill ch]

Dear Star Strider,

Thank you ,attached are my figures however i do not know how to evaluate it with interval of (-w,+w) ,if you have time please guide me for it.

Your help will be never forgotten. Thank you so much

Regards

Ill

[Star Strider]

I was thinking of something like this:

New_Time_Based=ifft(fit(s,[-w,w]),L,'symmetric');

figure
plot(t, New_Time_Based)
grid

That is likely as close as you can get with the function fit to your data. The model does not accurately reflect the data, and there are certain assumptions with respect to the Fourier transform that your data do not exactly conform to.

You may have to experiment to get the result you want.

[Ill ch]

Thank you for all

[Star Strider]

As always, my pleasure.

[Ill ch]

Dear Star Strider,

One more suggestion i need from you. Could you please guide me which frequency domain formula should be best to fit and to get after it time domain SIgnal that should be similar like original time domain signal. as a example at attached photo(original time domain)

Thank you

[Star Strider]

As always, my pleasure.

The plot you posted does not resemble the plot I get of the data you attached (‘MyArtificial.xlsx’) as the original time-domain signal:

In any event, the ifft of the fitted model is unlikely to faitthfully reproduce the original time-domain model, largely because it is of a fitted model.

[Ill ch]

Dear Star Strider

yes i checked it. I do not know why in your case data is plotting in that form when we use xlxc. However by me its in exponentially decaying form. I am using csv. as well xlxc data file for same values. with csv. data results are coming as i have attached figure. Today i optimize values and results are coming well good.

Is it good idea to optimize value more than one time? like

initial values what i am providing in first trial and getting new optimize values...now if i give these optimize values as a initial condition then result coming good.

if yes then could you please suggest me how i can do here loop to run two to three time by changing initial to new_1_optimize values and again new_1_optimize value to new_2_optimize value

your suggestions and guidance will be helpful for me

Thank you Sir

[Star Strider]

As always, my pleasure.

There is a significant problem with your ‘MyArtificial.xlsx’ file.

See for example lines 109 and 110 in the Excel file:

0.9177
10,267

The same problem appears other places as well. This may be a problem with the decimal separator going from ‘.’ to ‘,’, however I cannot say for sure. Please provide a consistent data vector, attach it to a Comment, and I will take another look at your data.

[Star Strider]

I was hoping for a file with that information, so I created one froim it. I am attaching it here.

You can use any appropriate optimisation function to estimate the parameters. I usually start with fminsearch with Answers posts (unless it is obvious that the person has the appropriate Toolbox functions) because everyone has it.

You can certainly use a loop experiment with different initial parameter estimates. There are built-in functions that do this. If you have the Global Optimization Toolbox, experiment with MultiStart (and similar functions, linked to in and at the end of that documentation page).

I still get a good fit with the new data, (that I still had to edit). The ifft of the evaluated model does not reproduce your original signal because the model itself does not describe your signal.

We also do not appear to be using the same data, since the time-domain plot of these data is:

This does not closely resemble the image of the data that you posted:

We are obviously not using the same data.

[Ill ch]

I am sorry there were still problem with *.* and *,* now here i am giving mat file . that is what exact what i am using without errors. Please could you have look on it. I do not have multistart function as i have not global optimization toolbox

Thank you for your kind support

[Star Strider]

The result you get with a parameter estimation of the new data is likely as good as you can hope for. Remember, your initial time-domain equation for your exponentially-decaying sinusiod does not actually describe the data you are modeling, so it will not actually fit your data in either the time-domain or frequency domain. I tried a time-domain regression with the new data, and even using the ga (genetic algorithm) function that searches the entire parameter spece, it failed to fit your function to your data. Unfortunately, it did not even come close.

[Ill ch]

Dear Star Strider,

Thank you very much for your polite way to explain and all help. I tried my best to get good result and i am getting as in attached figure.

I feel its better i cut initial part of original signal from starting so that my result will be close to original.from my side idea was to take original signal from first peak beacuse initial points are nearly 0.00123xx. i want to cut that part and start with first peak. could you please help me? how to do with loop i am trying it since yesterday but i dont know how i can do with loop.

i have noted in sig.jpg figure from that point i want to start my signal and before that point i would like to cut the signal, could you please help me to write better loop

Thank you very very much

[Star Strider]

As always, my pleasure.

To locate that point, then plot it:

StartPt = find(Measured_signal >= 0.43, 1, 'first')

figure
plot(t, Measured_signal, t(StartPt),Measured_signal(StartPt),'+')
grid

The ‘StartPt’ value is the index of that value in the vector. It is not the value itself.

Experiment to get the result you want.

[Ill ch]

Thank you very much for your reply

actually i wanted to move that point to initial or viceversa i add in synthetic signal zeros to check similarity between both.

i like to do but it automatically

could you have please look on figure attached

Thank you

[Star Strider]

If you want to use ‘StartPt’ as the starting point for your vectors, define your vectors as:

t2 = t(StartPt:end);
Measured_signal2 = Measured_signal(StartPt:end);

Then use ‘t2’ and ‘Measured_signal2’ in your analyses. You can of course define ‘StartPt’ however you want.

[Ill ch]

Thank you very much

[Star Strider]

As always, my pleasure.

[Star Strider]

I do not understand what you want to do with that.

[Star Strider]

I have tried everything I can think of, including a genetic algorithm, and I cannot fit your function to your data. I can get a good fit to the frequency, but not to the decaying exponential enveloope, even when I add a time-offset term.

Also, you need to use in both the exponential and cosine arguments. I did, and it still did not make any difference. I still could not get it to fit.

[Ill ch]

Dear Star Strider,

Thank you very much for your support and kind help.

I will have discussion with my senior supervisisors.In case i have further doubts then i will upload here.

Thank you so much and have a great weekend :)

[Star Strider]

As always, my pleasure.

You, too!

If I think of a way to fit your data, I will post back here.

[Ill ch]

sure Thank you.

[Star Strider]

I tried to fit your data in the time-domain using the code I referred you to earlier (in How to filter noise from time-frequency data and find natural frequency of a cantilever?) however, even including a time-shift term (so that ‘t’ was replaced by ‘(t-b(6))’ in both the exponential and cosine terms, I was unable to get a decent fit. I used fminsearch, lsqcurvefit, and ga, all without success. The only parameter I was able to estimate with any accuracy was the frequency of the cosine term.

[Ill ch]

hmm I used also same code what you had referred but unable to fit perfectly.

[Ill ch]

Dear Star Strider,

I am sorry i have asked many questions but thank you too for your kind support.

I had discussion with my senior. result what i got its sufficient however i need to add after getting inverse fft additionally new function as a example to give in strating initial values as a zero with angle . similarly both signal should have same angle as attached photo you can see. Similiarly attached you can check my code too.

I will be very thank ful to you if you can help me.

[Star Strider]

I recently considered:

instead of:

which is a definite improvement. However I still cannot model the delay correctly.

[Star Strider]

There is not much to solve. The error is that ‘fit’ is undefined (so likely Inf or NaN) with the initial parameter estimates. (You can determine that easily enough by running ‘fun’ using those initial parameter estimates as the values for ‘in1’.) The solution is to change them so that the value is finite.

[Star Strider]

Use this vector:

[1245,80,150,0,ym,0.033]

as the ‘in1’ argument in ‘fun’ to understand what the problem is. Then you can change it as necessary so ‘fun’ produces a finite result. Then use that vector as your initial parameter estimates.

[Ill ch]

Dear Star,

One small question: in the following code Tmax also should be Totalnumber of sample(N)/ Fs ? I mean to say here also same rule like Ts=1/Fs. or i can use any Tmax?

Thank you

syms k1 k2 t Tmax w w0 
y(t,w0,k1,k2) = exp(k1*t) * cos(k2*w0*t)
Y = int(y*exp(1j*w*t), t, 0, Tmax)
Y = simplify(Y, 'Steps',250)
Yf = matlabFunction(Y)                                      % Individual Parameters
Yf = matlabFunction(Y, 'Vars',{[w0,k,Tmax],w})              % Parameter Vector ‘in1’

[Star Strider]

If ‘Ts’ is the sampling interval and ‘Fs’ is the sampling frequency, ‘Ts=1/Fs’ is always true.

If you define ‘Tmax’ as being the maximum time, then it is likely true that ‘Totalnumber of sample(N)/Fs’ is the correct expression for it. This would be the same as ‘Totalnumber of sample(N)*Ts’. (This of course assumes that the sampling interval is the same for all the samples.)

As always, my pleasure.

[Ill ch]

thank you..... Awesome...

[Star Strider]

Again, my pleasure!

[Star Strider]

See: Import or Export a Sequence of Files and: How can I process a sequence of files?. Referring you to these is the best I can do.

[Ill ch]

Thank you

[Star Strider]

As always, my pleasure.

[Ill ch]

Hello Star ,

I have small quary if you can guide me. In the following i have equation and i want to make for it absolute. Is it possible to make absolute of the equation without values? or to evaluate imagenary and real part in in this form a+bi?

% A value always >0; w is vector; I am taking one as a variable eg.A and another all are constant Tmax, delta,phi w0 viceversa 
syms A delta phi w w0 Tmax Yf
eqn = Yf ==(A*(delta*cos(phi) + w0*sin(phi) + w*cos(phi)*1i))/(delta^2 + delta*w*2i - w^2 + w0^2) + (A*exp(-Tmax*delta)*(sin(Tmax*w) + cos(Tmax*w)*1i)*(delta*cos(phi - Tmax*w0) + w*cos(phi - Tmax*w0)*1i + w0*sin(phi - Tmax*w0))*1i)/(delta^2 + delta*w*2i - w^2 + w0^2)
 

[Star Strider]

The easiest way to find out is to ask MATLAB!

eqn2 = Yf == abs((A*(delta*cos(phi) + w0*sin(phi) + w*cos(phi)*1i))/(delta^2 + delta*w*2i - w^2 + w0^2) + (A*exp(-Tmax*delta)*(sin(Tmax*w) + cos(Tmax*w)*1i)*(delta*cos(phi - Tmax*w0) + w*cos(phi - Tmax*w0)*1i + w0*sin(phi - Tmax*w0))*1i)/(delta^2 + delta*w*2i - w^2 + w0^2));
eqnRe = Yf == real((A*(delta*cos(phi) + w0*sin(phi) + w*cos(phi)*1i))/(delta^2 + delta*w*2i - w^2 + w0^2) + (A*exp(-Tmax*delta)*(sin(Tmax*w) + cos(Tmax*w)*1i)*(delta*cos(phi - Tmax*w0) + w*cos(phi - Tmax*w0)*1i + w0*sin(phi - Tmax*w0))*1i)/(delta^2 + delta*w*2i - w^2 + w0^2));
eqnIm = Yf == imag((A*(delta*cos(phi) + w0*sin(phi) + w*cos(phi)*1i))/(delta^2 + delta*w*2i - w^2 + w0^2) + (A*exp(-Tmax*delta)*(sin(Tmax*w) + cos(Tmax*w)*1i)*(delta*cos(phi - Tmax*w0) + w*cos(phi - Tmax*w0)*1i + w0*sin(phi - Tmax*w0))*1i)/(delta^2 + delta*w*2i - w^2 + w0^2));

These do not automatically throw errors, however it is not possible to determine the results without knowing how you want to use them in calculations. They will likely require numerical values at some point, however they may not do what you expect them to do.

[Ill ch]

Great! Thank you very much. you are the best.:)

[Star Strider]

As always, my pleasure!

[Ill ch]

Dear star,

I need your small favour to solve my confusion. In this discussion portal we had found the curve fitting equation in frquency domain: (you can see your starting comments)

syms A delta phi t Tmax w w0 
y(t,w0,A,delta,phi) = A*exp(-delta*t) * cos(w0*t-phi)
Y = int(y*exp(1j*w*t), t, 0, Tmax)  

I am thinking since a day what is the unit of Amplitude A in frequency domain equation?

Basically my measuremets data are in voltage (in time domain). This data i am converting in frequency domain to fit the above shows frequency domain equation.

Delta has no unit. Time can be taken in second (i suppose), phase in rad but what about Amplitude?

Thaank you very much for your all helps which u have did....

Thank you in advance

[Star Strider]

As always, my pleasure.

The amplitude units are the same for frequency and time domain signals, unless they are deliberately changed. (For example, to use the Fourier transform for power spectral density, the amplitude is squared, so the units are power, not voltage. That is not done here, so I am using this only as an illustration.)

Time is in whatever units it is defined, and frequency is the inverse of those units, so frequency is in units of cycles/(time unit). The arguments to the exponential function (and for that matter, all transcendental functions) must by definition be dimensionless, so ‘delta’ is in units of inverse time (for example 1/second if ‘t’ is in seconds).

[Ill ch]

Hello Star,

First of all Thank you very much for your feedback and detailed explanation. I was thinking delta is dimensionless as it is just exponential coefficient.

Regarding amplitude as shown in the following figure curve fitting Amplitude A=84 (for black curve) and in the figure shows max amplitude 0.2xx on y axis!! I do not understand here what Amplitude mean !?

for remind our past discussion: short explanation on figure: 1. Measured data (voltage) converted to fft 2. then frequency domain formula were applied 3. figure shows the result with best fit.

[Star Strider]

I am not certain what the problem is. The fit appears to be reasonable.

I do not understand where the ‘A=84’ comes from or what it means.

[Ill ch]

There is no any problem. Just i am confuse with Amplitude A =84 which comes from curve fitting optimized result from the following equation.

syms A delta phi t Tmax w w0 
y(t,w0,A,delta,phi) = A*exp(-delta*t) * cos(w0*t-phi)
Y = int(y*exp(1j*w*t), t, 0, Tmax)  

Yf = matlabFunction(Y, 'Vars',{[w0,A,delta,phi,Tmax],w})    % Parameter Vector ‘in1’

I am confuse with amplitude unit. Why on y axis amplitude is very small. Or the amplitude which stets on y-axis is not the amplitude which i am getting from the equation from above?

from your sentense:

The amplitude units are the same for frequency and time domain signals, unless they are deliberately changed. (For example, to use the Fourier transform for power spectral density, the amplitude is squared, so the units are power, not voltage. That is not done here, so I am using this only as an illustration.)

That means in my case Amplitude unit is Voltage? A=84 seems very high voltage in equation.

Thank you again

[Star Strider]

I have long since lost track of what we were doing.

I have no idea why the amplitude should be 84 when it otherwise appears to be as small as it is in other analyses, such as the Fourier transform.

[Ill ch]

Hello Star Strider,

Almost, after one year again on matlab query :). I wish you are fine and very well in this COVID 19 time.

I have uploaded one more question on my current isssues. I hope you can help me there.

as i had forgotten this id so created new and uploaded there.

Here is link:https://de.mathworks.com/matlabcentral/answers/774192-desired-input-signal-generation-based-on-input-and-measured-output-and-expected-output-with-data-com?s_tid=srchtitle

Thank you very much in advance

Best regards,

Ic