Использование анализа Фурье для подгонки функции к данным

Question

У меня 24 значения для Y и соответствующие 24 значения для значений Y измеряются экспериментально,

, в то время как t имеет значения: t=[1,2,3........24]

Я хочу найти связь между Y иВ качестве уравнения с использованием анализа Фурье

я попытался сделать следующее:

Я написал следующий код 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

результаты гистограммыis:

Теперь я хочу найти уравнение для этих гармонических составляющих, которые представляют данные.После этого я хочу нарисовать исходные данные Y с данными, найденными из функции подгонки, и две кривые должны быть близки друг к другу.

Должен ли я использовать cos или sin или -sin или -cos?

Другими словами, каково правило для представления этих гармоник как функции: Y = f (t)?

Answer 1

Пример работы с вашими данными и 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}]

YouПридется немного поэкспериментировать ...

Answer 2

Проблема с строка 2 "t = [1: ts: 24];" должно быть "t = 0: ts: 23;"

Answer 3

Спасибо за вашу помощь.

Я нашел решение, к которому стремился, но почему-то все сдвинуто на 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)')

Почему произошла эта проблема, и как я могу ее решить?

Answer 4

Зависит от того, что MATLAB вернул вам.Это либо синус и косинус, либо сложная экспонента.

Большинство известных мне БПФ-алгоритмов обычно требуют, чтобы число точек данных было целой степенью двух.Самый близкий для вашего набора данных - 32, поэтому вы должны заполнить его нулями.