%-------------------------------------------------------------------------%
%--     固定縮小比例的FFT算法Matlab驗證程序。                             --%
%--         作者：Xudong                                                --%
%--         時間：2004年5月6日                                          --%
%------------------------------------------------------------------------%
function [g1out,g2out,g3out,g4out,g5out]=simfft(indata);
%--輸出數據：g1out...g5out都是8位有符號數（-128～+127）；
%--輸入數據：indata也是8位有符號數（-128～+127）；

PI=3.14159265357;
DataWidth=8; %--旋轉因子的精度（數據位寬）。
PointDot=32;   %--時間抽取基－2FFT的點數。

WnDot=PointDot/2; %--旋轉因子數組元素個數。
theta=[0:1/WnDot*PI:(WnDot-1)/WnDot*PI];
cos_ary=cos(theta);
sin_ary=sin(theta);

wn_r_t=round(cos_ary*2^(DataWidth-1)); %--四舍五入取整。
wn_i_t=-round(sin_ary*2^(DataWidth-1));%--此處要加“－”號，為便于后面處理。

for i=1:WnDot %--產生wn_r和wn_i數組（不用補碼來表示）。
    if wn_r_t(i)==2^(DataWidth-1)  %--避免對“+1”歸一化時的錯誤。
        wn_r(i)=wn_r_t(i)-1;
    else
        wn_r(i)=wn_r_t(i);
    end;
    if wn_i_t(i)==2^(DataWidth-1)
        wn_i(i)=wn_i_t(i)-1;
    else
        wn_i(i)=wn_i_t(i);
    end;    
end;
wn=wn_r+wn_i*sqrt(-1);%--構建復數形式的旋轉因子。

for i=1:PointDot %--產生仿真波形中的輸入數組（不要用補碼來表示）。
    if indata(i)==2^(DataWidth-1)  %--避免對“+1”歸一化時的錯誤。
        indata_comp(i)=indata(i)-1;
    else
        indata_comp(i)=indata(i);
    end;
end;
indata_bitrev=bitrevorder(indata_comp);%--輸入數據位反轉。
%--第一級輸出。
for i=1:2:PointDot
        [g1out(i),g1out(i+1)]=butfly(indata_bitrev(i),indata_bitrev(i+1),wn(1));
end;
%--第二級輸出。
for i=1:2%--外循環：每組蝶形結的個數。
    for j=1:4:PointDot%--內循環：每組蝶形結需要循環的次數。
        [g2out(i-1+j),g2out(i-1+j+2)]=butfly(g1out(i-1+j),g1out(i-1+j+2),wn(8*(i-1)+1));
    end;
end;
%--第三級輸出。
for i=1:4
    for j=1:8:PointDot
        [g3out(i-1+j),g3out(i-1+j+4)]=butfly(g2out(i-1+j),g2out(i-1+j+4),wn(4*(i-1)+1));
    end;
end;
%--第四級輸出。
for i=1:8
    for j=1:16:PointDot
        [g4out(i-1+j),g4out(i-1+j+8)]=butfly(g3out(i-1+j),g3out(i-1+j+8),wn(2*(i-1)+1));
    end;
end;
%--第五級（最后一級）輸出。
for i=1:16
        [g5out(i),g5out(i+16)]=butfly(g4out(i),g4out(i+16),wn(i));
end;
%--經過驗證第一級和第五級的輸出完全符合ISE與Modelsim中的結果；
%--幅度誤差與Matlab中算出的理論值也相差很小，可是相位誤差卻很大，原因有待進一步分析。
%--模值小的點相位誤差大，而模值大的點，例如仿真正弦信號頻率點處的相位誤差居然也很大。

%-------------------------------------------------------------------------%
%--幅度理論值：
%Columns 1 through 8 
%
%         0    1.3212    4.3405   61.9818    7.7471    4.0671    2.8353    2.3614
%
%  Columns 9 through 16 
%
%    1.9445    1.8409    1.5481    1.4261    1.3912    1.2839    1.2585    1.3036
%
%  Columns 17 through 24 
%
%    1.2500    1.3036    1.2585    1.2839    1.3912    1.4261    1.5481    1.8409
%
%  Columns 25 through 32 
%
%    1.9445    2.3614    2.8353    4.0671    7.7471   61.9818    4.3405    1.3212

%--“abs(g5out)”幅度值（與Modelsim中的輸出值不完全吻合，主要是CORDIC模塊的誤差，
    %但是g5out與Modelsim中算的FFT值是完全吻合的）：
%Columns 1 through 8 
%
%    4.0000    2.8284    2.0000   58.5235    8.0623    6.4031    4.4721    4.4721
%
%  Columns 9 through 16 
%
%    2.8284    3.6056    2.8284    2.8284    1.4142    2.2361    2.2361    2.2361
%
%  Columns 17 through 24 
%
%         0    2.0000    2.0000    2.2361    1.0000    2.2361    2.8284    2.8284
%
%  Columns 25 through 32 
%
%    2.8284    2.2361    2.8284    4.4721    7.6158   61.6847    3.6056    1.0000

%--“g5out”與“fftout_xil”比較：
%--g5out
%Columns 1 through 8 
%
%  -4.0000            -2.0000 - 2.0000i        0 + 2.0000i  17.0000 +56.0000i  -4.0000 - 7.0000i  -4.0000 - 5.0000i  -2.0000 - 4.0000i  -2.0000 - 4.0000i
%
%  Columns 9 through 16 
%
%  -2.0000 - 2.0000i  -3.0000 - 2.0000i  -2.0000 - 2.0000i  -2.0000 - 2.0000i  -1.0000 - 1.0000i  -2.0000 - 1.0000i  -2.0000 - 1.0000i  -2.0000 - 1.0000i
%
%  Columns 17 through 24 
%
%        0            -2.0000            -2.0000            -1.0000 + 2.0000i        0 + 1.0000i  -2.0000 + 1.0000i  -2.0000 + 2.0000i  -2.0000 + 2.0000i
%
%  Columns 25 through 32 
%
%  -2.0000 + 2.0000i  -1.0000 + 2.0000i  -2.0000 + 2.0000i  -2.0000 + 4.0000i  -3.0000 + 7.0000i  18.0000 -59.0000i   2.0000 - 3.0000i        0 - 1.0000i
%--fftout_xil
% Columns 1 through 8 
%
%        0             0.1295 - 1.3148i   0.8468 - 4.2571i  17.9924 -59.3129i  -2.9647 + 7.1574i  -1.9172 + 3.5868i  -1.5752 + 2.3575i  -1.4980 + 1.8254i
%
%  Columns 9 through 16 
%
%  -1.3750 + 1.3750i  -1.4230 + 1.1678i  -1.2872 + 0.8601i  -1.2577 + 0.6723i  -1.2853 + 0.5324i  -1.2286 + 0.3727i  -1.2343 + 0.2455i  -1.2973 + 0.1278i
%
%  Columns 17 through 24 
%
%  -1.2500            -1.2973 - 0.1278i  -1.2343 - 0.2455i  -1.2286 - 0.3727i  -1.2853 - 0.5324i  -1.2577 - 0.6723i  -1.2872 - 0.8601i  -1.4230 - 1.1678i
%
%  Columns 25 through 32 
%
%  -1.3750 - 1.3750i  -1.4980 - 1.8254i  -1.5752 - 2.3575i  -1.9172 - 3.5868i  -2.9647 - 7.1574i  17.9924 +59.3129i   0.8468 + 4.2571i   0.1295 + 1.3148i

%--由此可見這是算法本身的缺陷導致了后面計算相位時的誤差是如此的大。下面要分析產生如此大的相位誤差的原因所在，以及能夠采用的改進方法。
   %ffttheorytest.m文件用于分析相位誤差的原因所在。2004年5月7日上午。
%--原因已經查明：是Wi中的符號差了一個負號的原因。修正后的結果在信號頻點處的幅度和相位誤差都已經達到了很好的精度，在沒有信號的頻點處
  %的相位誤差仍然很大，這是因為這些點的模值太小，相位受截斷誤差的影響是不可避免的。2004年5月7日14:16。