2011-03-26 28 views
1

我对Y代表Y值24个值和对应的24个值实验测得,利用傅立叶分析,以适应功能,数据

而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 

酒吧图的结果是:

enter image description here

现在,我想找到用于表示该数据,这些谐波分量的方程。之后,我想用拟合函数中找到的数据绘制原始数据Y,两条曲线应该彼此接近。

我应该使用cos还是sin或-sin或-cos?

换句话说,将这些谐波表示为函数的规则是什么:Y = f (t)

回答

0

取决于MATLAB给你回来的东西。它可以是正弦和余弦,也可以是复数指数。我知道大多数FFT算法通常要求数据点的数量是2的整数次幂。数据集最近的一个是32,所以你应该用零填充。

+1

用0填充会显着改变FFT的结果,这可能会隐藏有用的信号。如果你可以根据你的理由抽取2的幂数,但是如果你有数据,使用你的数据。然而,这表示,http://en.wikipedia.org/wiki/Bluestein%27s_FFT_algorithm允许任何大小的FFT(包括素数)在时间O(n log(n))内发生。所以你实际上不必使用2的幂。 – btilly 2011-03-27 03:01:10

4

使用离散正弦变换完成您的数据和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]] 

enter image description here

的系数为:

{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} 

更详细的视图:

enter image description here

编辑

然而,作为偶函数似乎更好,我还做了3类型,它的效果要好得多的离散傅立叶余弦变换:

enter image description here

在这种情况下,系数为:

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

而且coeffs和功能的标绘通过获得

你必须尝试一点点...

0

感谢您的帮助。

,我发现我的目标是获得解决方案,但由于某种原因,一切都是由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)') 

为什么会发生这种转移的问题发生了,我该怎么解决它?

0

问题出在 第2行 “t = [1:ts:24];” 它应该是“t = 0:ts:23;”