儘管我可以繪製我的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))
那裏有一個'hold'嗎? – Schorsch
你是什麼意思?繪製f(x)是我遇到的問題,但它似乎有一些參數問題x – Crisp
@Crisp:您需要設置'hold on',否則第二個plot命令將刪除第一個線。 http://www.mathworks.de/de/help/matlab/ref/hold.html – Daniel