MATLABによるマイケルソン干渉計のシミュレーション

お世話になっております。

低コヒーレンス光を用いたマイケルソン干渉計による光路差の測定のシミュレーションをMATLABでやっています。簡単に言うと干渉信号の位相差から光路差を求めるのですが、その位相差がうまく算出されません。

まず、理論の説明です。
簡単のために、位相差によって光路差が求められることの確認をします。
低コヒーレンス光は波長780nm~820nmまでの41本の波の合成波を使います。
但し、この波のスペクトル分布はガウス分布に従います。(ソースコードfigure(4)参照)

(スペクトル分布、FFT、位相分布は横軸がλの逆数として、1/820~1/780が該当する部分を切り取っています、単位がμmなので1/0.82〜1/0.78ですけど…)

今、初期光路差が参照光と物体光で0とした時、参照ミラーが-400μm=400μmまで動いたと考えて、例えば波長780nmの波は0μmの地点で必ず光強度最大となります。もちろん、この波長のみでしたら他の地点でも最大光強度はあります。
ここで、シミュレーションの都合上、41本すべての波に関してそれぞれ干渉信号を求め、その41本の波を足し合わせると0μmの地点で光強度最大となります。(ソークコードfigure(1)参照)

これをFFTして最大振幅を求め、その最大振幅の位相の傾きから位相差を求めることができますが、今回は初期光路差が0なので当然FFTの最大の点での位相αは0です。

位相αを用いてα=2πd/λより光路差dは0と求められます。

つまり、初期光路差が0の場合はFFTの最大点での位相は0でなければならないのですが、下のプログラムのように書いて実行するとギザギザした感じになってしまいます。(ソースコードfigure(3)参照)
FFTまではうまくいっていると思うのですが、どうもよくわかりません。

わかる方がいらっしゃいましたらご教授お願いいたします。

以下ソースコードです。
コピーして動かして頂くと詰まっているところがわかっていただけるかと思います。
変数の名前が完全に適当で申し訳ございません。

＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿
追記:位相を求めるところって合っていますか？atan2とatanを使った時で値が変わるんですけど…
＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿

MATLAB
1
2for i = 780:1:820;
3    test(1,(i-779)) = 1000./i
4end
5a = mean(test);
6b = sqrt(var(test));
7
8test1 = normpdf(test,a,b/2)  %スペクトル分布
9for i = 780:1:820;
10    x(1,(i-779)) = 1000./i;
11end
12
13figure(4);
14plot(x,test1)  %横軸は波長λの逆数(単位μm-1)
15title('スペクトル分布(横軸はλの逆数.単位μm-1)')
16
17
18for i = 1:1:41;
19    
20    for j = 1:1:8192;
21    si(i,j) = test1(1,i).*cos(2.*pi./(0.78 + (i .* 0.001)).*(j.*0.1 - 400));
22    %41本の干渉信号作成
23    %サンプリング距離は0.1μm,データ点数は8192
24    end
25    
26end
27
28for i = 1:1:8192;
29    y(1,i) = i.*0.1 - 400;  %干渉信号の横軸計算(実際の実験では参照ミラーの移動距離.μm)
30end
31
32S = sum(si)./41;
33figure(1); %干渉信号
34plot(y,S);
35xlim([-100,100]);%本当は横軸が±400μmほどのデータがあるが、見やすさのために切り取り
36
37z = fft(S);
38m = abs(z)./41;
39p = (angle(z)); %位相(FFTの絶対値最大の点での位相が0になってほしい…)
40for i = 1:1:8192;
41    la(1,i) = 10./8192.*i;  %横軸
42end
43
44figure(2);
45plot(la,m);
46title('FFT(横軸は1/0.82~1/0.78までを切り取り)')
47xlim([1.221001221,1.280409731]); 
48
49
50figure(3);
51plot(la,p);
52title('位相分布(横軸は1/0.82~1/0.78までを切り取り)')
53xlim([1.221001221,1.280409731]);
54