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

问题描述

| 我有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,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],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,1.5),grid
xlabel(\'frequency (Hz)\');ylabel(\'amplitude (MW)\')
subplot(212),phi,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,1.5);grid on
hold on;
plot(t,Pmodel,\'r\',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,1.5);grid on
hold on;
plot(t+1,\'model\');xlabel(\'time (hours )\');ylabel(\'Power (MW)\')

为什么会发生这种变化的问题，我该如何解决？,问题出在 2号线 \“ t = [1：ts：24]; \” 它应该是\“ t = 0：ts：23; \”