2014-04-23 75 views
0

尽管我可以绘制我的FFT(快速傅里叶变换)图(X,Y),但是我无法绘制我的拟合线f(x)以及我的FFT 。从曲线拟合工具收集方程,我取十个最佳拟合并求平均值以得到方程f(x)。在同一图形/ matlab上绘制两个不同的公式

我必须对f执行(X)与数(频率)来绘制和日志(POW)

代码:因此

row=69; 
col=69; 

colormap gray 
whitebg('black') 

iterations=10^3; 


Next=zeros(row,col); 
laplacian=zeros(row,col); 
critical=zeros(row,col); 

B= zeros(row,col); 
lums=zeros(1000); 


flw=0.5; 

u=0.1; 

crit=5; 
%bns=200; 
bns=1000; 

for k=1:iterations 
    B=B+(rand(row,col)-0.5); 
    Next=B; 

rns=5.; 
for i=1:row 
for j=1:col 

rfromns=(col+rns-j); 
     critical(i,j)=0; 
     if i<=2 left=row; else left=(i-1);end 
     if i==row right=1; else right=(i+1);end 
     if (j<=2) top=1; else top=(j-1);end 
     if (j==col) bottom=j; else bottom=(j+1);end 
     l=B(i,top)+B(left,j)+B(right,j)+B(i,bottom)+0.5*(B(left,top)+B(right,top)+B(left,bottom)+B(right,bottom)); 
     l=l-6.*B(i,j); 
     laplacian(i,j)=l; 
     bfromns=bns/rfromns^3; 
     if (abs(l)/((abs(B(i,j))+bfromns)+0.01))>crit; critical(i,j)=1; end 
     %if abs(l)>crit; critical(i,j)=1; end 
end 
end 



      for j = 1:col 
      if (j==col) lum=0.;end 
      for i = 1:row 

      if (j>1) Next(i,j)=(1-flw)*B(i,j)+flw*B(i,j-1); end; 
      if (j==1) Next(i,j)=(1-flw)*B(i,j); end; 

      if (critical(i,j)==1)&& (j>1) Next(i,j)=B(i,j)*(1-flw-flw*u)+flw*B(i,j-1)+(flw*u)*B(i,j)/5.; end; 
      if i<2 left=row; else left=(i-1);end 
      if i==row right=1; else right=(i+1);end 
      if (j<=2) top=1; else top=(j-1);end 
      if (j==col) bottom=j; else bottom=(j+1);end 

      if (critical(left,j)==1) Next(i,j)=Next(i,j)+flw*u*B(left,j)/5.;end 
      if (critical(right,j)==1) Next(i,j)=Next(i,j)+flw*u*B(right,j)/5.;end 
      if (critical(i,top)==1) Next(i,j)=Next(i,j)+flw*u*B(i,top)/5.;end 
      if (critical(i,bottom)==1) Next(i,j)=Next(i,j)+flw*u*B(i,bottom)/5.;end 

      if (j==col) lum=lum+u*B(i,j)*B(i,j)+abs(laplacian(i,j)); end 
      end 
      end 

lums(k)=lum; 

B=Next;     
%Matplot(B) 
%if (k>00) 
surf(B); 
%plot(lums) 
%view([0 90]) 
%pause(0.001) 
%end 
end 


c=fft(lums(129:iterations)); 
pow=abs(c).^2; 
pow=pow(2:(iterations-128)/2); 
freq=(2:(iterations-128)/2); 

X=log(freq); 
Y=log(pow); 

%x=length(X); 

x=0.6:.1:6.; 

%Linear model Poly2 
a1 = -0.155; 
a2 = 0.2154; 
a3 = 15.1; 
af(x) = a1*x.^2 + a2*x + a3; 

%Linear model Poly3 
b1 = 0.01805; 
b2 = -0.3687; 
b3 = 0.9874; 
b4 = 14.29; 
bf(x) = b1*x.^3 + b2*x.^2 + b3*x + b4; 


%General model Power2 
c1 = -0.09124; 
c2 = 2.179; 
c3 = 15.34; 
cf(x) = c1*x.^c2+c3; 


%General model Rat02 
d1 = 727.3; 
d2 = -3.447; 
d3 = 51.6; 
df(x) = (d1)/(x.^2 + d2*x + d3); 


%General model Gauss1 
e1 = 15.01; 
e2 = 1.346; 
e3 = 8.152; 
ef(x) = e1*exp(-((x-e2)/e3).^2); 


%General model Gauss2 
w1 = 1.737; 
w2 = 3.46; 
w3 = 2.333; 
w4 = 30.03; 
w5 = -23.14; 
w6 = 28.23; 
wf(x) = w1*exp(-((x-w2)/w3).^2) + w4*exp(-((x-w5)/w6).^2); 


%General model Sin1 
g1 = 15.11; 
g2 = 0.1526; 
g3 = 1.428; 
gf(x) = g1*sin(g2*x+g3); 



%Linear model Poly4 
h1 = 0.0179; 
h2 = -0.252; 
h3 = 1.047; 
h4 = -1.97; 
h5 = 16.23; 
hf(x) = h1*x.^4 + h2*x.^3 + h3*x.^2 + h4*x + h5; 


%General model Fourier1 
m1 = 11.05; 
m2 = 3.31; 
m3 = 2.104; 
m4 = 0.3644; 
mf(x) = m1 + m2*cos(x*m4) + m3*sin(x*m4); 


%Linear model 
p1 = 0.815; 
p2 = 0.1061; 
p3 = 8.904; 
pf(x) = p1*(sin(x-pi)) + p2*((x-10).^2) + p3; 


f(x)=(af(x)+bf(x)+cf(x)+df(x)+ef(x)+wf(x)+gf(x)+hf(x)+mf(x)+pf(x))/10; 



plot(X,Y) 
plot(f(x)) 
+2

那里有一个'hold'吗? – Schorsch

+0

你是什么意思?绘制f(x)是我遇到的问题,但它似乎有一些参数问题x – Crisp

+0

@Crisp:您需要设置'hold on',否则第二个plot命令将删除第一个线。 http://www.mathworks.de/de/help/matlab/ref/hold.html – Daniel

回答

1

基于Matlab /倍频使用hold关键字。

如果您想在一张图上绘制几件事情,您可以使用hold on开始您的程序,然后执行一个或多个绘图命令,并使用hold off完成。

实施例:

hold on; 
x = -10:0.1:10; 
plot (x, sin (x)); 
plot (x, cos (x)); 
hold off; 

文档:https://www.gnu.org/software/octave/doc/interpreter/Manipulation-of-Plot-Windows.html


作为文档描述,plot通常会调用newplot命令,去除先前的积结果,hold on;这种行为被防止。

+0

好的谢谢你,将这样做 – Crisp

+0

好轻微不同,但我得到了这个 plot( (f(x)) plot(f(x)) – Crisp

相关问题