2014-04-23 112 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

相關問題