使用傅里叶分析使函数适合数据

| 我有24个Y值,相应的24个Y值是通过实验测量的, 而t的值是:
t=[1,2,3........24]
我想使用傅立叶分析找到Y和t之间的关系作为方程式, 我尝试过的是: 我编写了以下MATLAB代码:
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
条形图的结果是: 现在,我想找到这些代表数据的谐波分量的方程式。之后,我想用从拟合函数中找到的数据绘制原始数据Y,两条曲线应彼此靠近。 我应该使用cos或sin还是-sin或-cos? 换句话说,将这些谐波表示为函数的规则是什么:is2ѭ?
已邀请:
使用离散正弦变换对数据和Mathematica进行处理的示例。希望您可以推断到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]]
系数为:
{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}
更详细的视图: 编辑 但是,由于偶数函数似乎更好,所以我还进行了类型3的离散傅立叶余弦变换,其效果更好: 在这种情况下,系数为:
{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}
并通过以下公式获得系数和函数的图:
coef  = FourierDCT[fg, 3]/Sqrt[n];(*Fourier transform*)
f[x_]:= Sum[coef[[r]]*Cos[Pi (r - 1/2) x], {r, n - 1}]
您将不得不尝试一下...
取决于MATLAB给您的回报。它可以是正弦和余弦,也可以是复杂的指数。 我知道的大多数FFT算法通常都要求数据点的数量为2的整数次幂。您的数据集最接近的一个是32,因此您应该用零填充它。
谢谢你的帮助。 我找到了我想要的解决方案,但由于某种原因,一切都偏移了1 这是代码:
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)\')
为什么会发生这种变化的问题,我该如何解决?
问题出在 2号线 \“ t = [1:ts:24]; \” 它应该是\“ t = 0:ts:23; \”

要回复问题请先登录注册