Using Fourier Analysis to fit function to data

问题:

I have 24 values for Y and corresponding 24 values for the Y values are measured experimentally,

while t has values : t=[1,2,3........24]

I want to find the relationship between Y and t as an equation using Fourier analysis,

what I have tried and done is:

I wrote the following MATLAB code:

Y=[10.6534
    9.6646
    8.7137
    8.2863
    8.2863
    8.7137
    9.0000
    9.5726
   11.0000
   12.7137
   13.4274
   13.2863
   13.0000
   12.7137
   12.5726
   13.5726
   15.7137
   17.4274
   18.0000
   18.0000
   17.4274
   15.7137
   14.0297
   12.4345];

ts=1; % step

t=1:ts:24; % the period is 24 

f=[-length(t)/2:length(t)/2-1]/(length(t)*ts); % computing frequency interval

M=abs(fftshift(fft(Y)));

figure;plot(f,M,'LineWidth',1.5);grid % plot of harmonic components

figure;

plot(t,Y,'LineWidth',1.5);grid % plot of original data Y

figure;bar(f,M);grid % plot of harmonic components as bar shape

the results of the bar figure is:

Now, I want to find the equation for these harmonic components which represent the data. After that I want to draw the original data Y with the data found from the fitting function and the two curves should be close to each other.

Should I use cos or sin or -sin or -cos?

In another way, what is the rule to represent these harmonics as a function: Y = f (t) ?

回答1:

An example done with your data and Mathematica using Discrete sine transform. Hope you can extrapolate to Matlab:

n = 24;
xg = N[Range[n]]/n
fg = l                             (*your list *)

fp = ListPlot[Transpose[{xg, fg}], PlotRange -> All] (*points plot*)

coef = FourierDST[fg, 1]/Sqrt[n/2]; (*Fourier transform*)

Show[fp, Plot[Sum[coef[[r]]*Sin[Pi r x], {r, n - 1}], {x, -1, 1}, 
  PlotRange -> All]]

The coefficients are:

{16.6411,    -4.00062,    5.31557, -1.38863,    2.89762,    0.898562,
  1.54402,   -0.116046,   1.54847,  0.136079,   1.16729,    0.156489,   
  0.787476,  -0.0879736,  0.747845, 0.00903859, 0.515012,   0.021791,   
  0.35001,    0.0159676,  0.215619, 0.0122281,  0.0943376, -0.00150218}

More detailed view:

Edit

However, as an even function seems to be better, I made also a discrete fourier cosine transform of type 3, which works much better:

In this case the coefficients are:

{14.7384,  -8.93197,   4.56404,  -2.85262,   2.42847,   -0.249488, 
  0.565181,-0.848594,  0.958699, -0.468337,  0.660136,  -0.317903, 
  0.390689,-0.457621,  0.427875, -0.260669,  0.278931,  -0.166846, 
  0.18547, -0.102438,  0.111731, -0.0425396, 0.0484102, -0.00559378}

And the plotting of coeffs and function are obtained by:

coef  = FourierDCT[fg, 3]/Sqrt[n];(*Fourier transform*)
f[x_]:= Sum[coef[[r]]*Cos[Pi (r - 1/2) x], {r, n - 1}]

You'll have to experiment a little ...

回答2:

Depends on what MATLAB gave you back. It's either sine and cosine or a complex exponential.

Most FFT algorithms that I know of usually demand that the number of data points be an integer power of two. The closest one for your data set is 32, so you should pad it out with zeros.

回答3:

Thanks for your help.

I found the solution I was aiming to get but for some reason everything is shifted by 1

Here is the code:

ts = 1; % time step
t = [1:ts:24]; 
fs = 1/ts; % frequency step
f=[-length(t)/2:length(t)/2-1]/(length(t)*ts); % frequency formula

%data
P=[10.7083
    9.7003
    8.9780
    8.4531
    8.1653
    8.2633
    8.8795
    9.9850
   11.3289
   12.5172
   13.2012
   13.2720
   12.9435
   12.6647
   12.8940
   13.8516
   15.3819
   17.0033
   18.1227
   18.3039
   17.4531
   15.8322
   13.9056
   12.1154];

plot(t,P,'LineWidth',1.5);grid
xlabel('time (hours)');ylabel('Power (MW)')
title('Power Profile for 2nd Feb, 1998')

% fourier transform analysis
P1 = fft(P)/length(t);
P2=fftshift(P1);
amp=abs(P2); % amplitude
phi = angle(P2); % phase angle

figure
subplot(211),stem(f,amp,'LineWidth',1.5),grid
xlabel('frequency (Hz)');ylabel('amplitude (MW)')
subplot(212),stem(f,phi,'LineWidth',1.5),grid
xlabel('frequency (Hz)');ylabel('phase angle (rad)')


% NOW, I WILL CONSTRUCT THE MODEL FROM THE FIGURE
% THE STRUCTURE IS:
% Pmodel=Ai*COS(i*w*t+phii)
% where,  w=2*pi/24  and  i is the harmonic order
% Here, up to the third harmonic is enough
% and using Parseval's Theorem, the model is:

% PP=12.6635+2*(1.9806*cos(w*tt+1.807)+0.86388*cos(2*w*tt+2.0769)+0.39683*cos(3*w*tt-    1.8132));

w=2*pi/24;

Pmodel=12.6635+2*(1.9806*cos(w*t+1.807)+0.86388*cos(2*w*t+2.0769)+0.39686*cos(3*w*t-1.8132));

figure
plot(t,P,'LineWidth',1.5);grid on
hold on;
plot(t,Pmodel,'r','LineWidth',1.5)
legend('original','model');xlabel('time (hours )');ylabel('Power (MW)')

% But here is a problem, the modeled signal is shifted
% by 1 comparing to the original one
% I redraw the two figures together by plotting Pmodeled vs t+1
% Actually, I don't know why it is shifted, but they are 
% exactly identical with shifting by 1

figure
plot(t,P,'LineWidth',1.5);grid on
hold on;
plot(t+1,Pmodel,'r','LineWidth',1.5)
legend('original','model');xlabel('time (hours )');ylabel('Power (MW)')

Why has this shifting problem happened, and how can I solve it?

回答4:

The problem is with line 2 "t = [1:ts:24];" it should be "t= 0:ts:23;"