**You are now following this question**

- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.

# frequency equation for curve fitting

60 views (last 30 days)

Show older comments

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

##### 0 Comments

### Accepted Answer

Star Strider
on 9 Sep 2019

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!

##### 53 Comments

Ill ch
on 9 Sep 2019

Edited: Ill ch
on 28 Oct 2019

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
on 9 Sep 2019

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
on 9 Sep 2019

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.

Star Strider
on 9 Sep 2019

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
on 10 Sep 2019

Edited: Ill ch
on 10 Sep 2019

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
on 10 Sep 2019

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.

Star Strider
on 10 Sep 2019

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.

Star Strider
on 10 Sep 2019

As always, my pleasure.

Ill ch
on 11 Sep 2019

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
on 11 Sep 2019

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
on 11 Sep 2019

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
on 11 Sep 2019

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
on 12 Sep 2019

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
on 12 Sep 2019

Edited: Ill ch
on 12 Sep 2019

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
on 12 Sep 2019

Ill ch
on 13 Sep 2019

Edited: Ill ch
on 13 Sep 2019

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
on 13 Sep 2019

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
on 13 Sep 2019

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
on 13 Sep 2019

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.

Star Strider
on 13 Sep 2019

As always, my pleasure.

Star Strider
on 14 Sep 2019

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

Star Strider
on 14 Sep 2019

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
on 14 Sep 2019

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
on 14 Sep 2019

As always, my pleasure.

You, too!

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

Star Strider
on 16 Sep 2019

Ill ch
on 20 Sep 2019

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
on 20 Sep 2019

I recently considered:

instead of:

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

Star Strider
on 23 Sep 2019

Star Strider
on 23 Sep 2019

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
on 8 Oct 2019

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
on 8 Oct 2019

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.

Star Strider
on 8 Oct 2019

Again, my pleasure!

Star Strider
on 12 Oct 2019

Star Strider
on 13 Oct 2019

As always, my pleasure.

Ill ch
on 15 Oct 2019

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
on 15 Oct 2019

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.

Star Strider
on 15 Oct 2019

As always, my pleasure!

Ill ch
on 19 Dec 2019

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
on 19 Dec 2019

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
on 19 Dec 2019

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
on 19 Dec 2019

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
on 19 Dec 2019

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
on 19 Dec 2019

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
on 16 Mar 2021

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.

Thank you very much in advance

Best regards,

Ic

### More Answers (0)

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!**An Error Occurred**

Unable to complete the action because of changes made to the page. Reload the page to see its updated state.

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list:

## How to Get Best Site Performance

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

### Americas

- América Latina (Español)
- Canada (English)
- United States (English)

### Europe

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)

### Asia Pacific

- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)