%-------------------------------------------------------------------------%
%--     時間抽取32點基－2fft的matlab驗證程序。                          --%
%--         作者：Xudong                                                --%
%--         時間：2004年5月4日                                          --%
%-------------------------------------------------------------------------%
clear;
clc;

PI=3.14159265357;
DWIDTH=8;
POINTS=32;

%i=sqrt(-1);
%indata=[39+58*i, 37+58*i, 36+5*i,  33+5*i,  36+69*i, 22+99*i, 22+65*i, 22+67*i,...%--仿真輸入數據（波形文件中的數據）。
%        22+58*i, 33+53*i, 32+51*i, 32+54*i, 32+23*i, 35+21*i, 52+14*i, 96+25*i,...
%        55+22*i, 25+34*i, -20+21*i,22+20*i, 22+21*i, 22+58*i, 58+98*i, 55+56*i,...
%        69+68*i, 5+57*i,  57+59*i, 58+63*i, 5+33*i,  6+32*i,  9+31*i,  7+38*i];

%indata=[indata_inv(1) indata_inv(17) indata_inv(9)  indata_inv(25) indata_inv(5) indata_inv(21) indata_inv(13) indata_inv(29)...
%        indata_inv(3) indata_inv(19) indata_inv(11) indata_inv(27) indata_inv(7) indata_inv(23) indata_inv(15) indata_inv(31)...
%        indata_inv(2) indata_inv(18) indata_inv(10) indata_inv(26) indata_inv(6) indata_inv(22) indata_inv(14) indata_inv(30)...
%        indata_inv(4) indata_inv(20) indata_inv(12) indata_inv(28) indata_inv(8) indata_inv(24) indata_inv(16) indata_inv(32)];
%        indata_inv是前面的indata。

%--下面的程序用于產生仿真所需的輸入信號（正弦信號）。
%                  ^
%                  |       *
%                  |    *     *
%                  |  *         *
%                  | *            *               
%                --|*--------------*--------------*----->
%                  |0 1 ...         *            *31
%                  |                  *         *
%                  |                    *     *
%                  |                       *

MULT_F=3;
OMIGA=2*PI/(POINTS-1)*MULT_F;
%indata=round(sin(0:OMIGA*1:2*PI*MULT_F)*2^(DWIDTH-1));
N=[0:1:31];
%indata=round(sin(OMIGA*N)*2^(DWIDTH-1));
indata=round(sin(OMIGA*N)*2^(6-1));%--考慮A/D動態范圍。
for i=1:POINTS %--產生仿真波形中的輸入數組（要用補碼來表示）。
    if indata(i)<0
        indata_comp(i)=2^DWIDTH+indata(i); %--取負數的補碼。
    elseif indata(i)==2^(DWIDTH-1)  %--避免對“+1”歸一化時的錯誤。
        indata_comp(i)=indata(i)-1;
    else
        indata_comp(i)=indata(i);
    end;
end;

fftout_mat=fft(indata);%--Malab計算出的FFT理論輸出。
fftout_xil=fftout_mat/(2^log2(POINTS));%--因為在設計中為了防止每級蝶形運算的數據逸出，
                                           %--所以每級輸出數據縮小了2倍。
%fftout_xil_sample=abs([fftout_xil(1) fftout_xil(17) fftout_xil(9)  fftout_xil(25) fftout_xil(5) fftout_xil(21) fftout_xil(13) fftout_xil(29) ...
%                   fftout_xil(3) fftout_xil(19) fftout_xil(11) fftout_xil(27) fftout_xil(7) fftout_xil(23) fftout_xil(15) fftout_xil(31)]);   
           
absout=abs(fftout_xil);%--理論FFT結果的幅度值。
thetaout=angle(fftout_xil);%--理論FFT結果相位值（弧度表示）。
%--與thetaout=angle(fftout_mat)效果是一樣的。

thetaout_xil=thetaout*2^(DWIDTH-3);%--按Xilin的CORDIC核中定義的數據格式縮放后的數值（波形中的理論數值，便于比較用）。

%thetaout_xil_sample=[thetaout_xil(1) thetaout_xil(17) thetaout_xil(9)  thetaout_xil(25) thetaout_xil(5) thetaout_xil(21) thetaout_xil(13) thetaout_xil(29) ...
%                   thetaout_xil(3) thetaout_xil(19) thetaout_xil(11) thetaout_xil(27) thetaout_xil(7) thetaout_xil(23) thetaout_xil(15) thetaout_xil(31)];

thetaout_du=thetaout/PI*180;%--度數表示的FFT相位理論輸出結果。

%--下面的程序是為了將Modelsim中的仿真結果變換成角度的表示，以便更好的分析誤差。
%--xilsimtheta_r=[  4850 -15678 -16190 -18165 -20950 -23310 -24708 -23941 -21972 ...
%--               -19181 -19927 -21461 -22926 -22415 -20438 -19181 -20182 -21972 ...
%--               -21972 -22415 -21972 -20438 -20183 -20950 -21972 -21460 -20438 ...
%--               -20182 -20182 -21972 -21972 -21972 -19181 -20182 -20182 -21460 ...
%--               -22415 -20438 -20438 -19181 -20182 -20950 -21460 -20438 -20438 ...
%--               -18421 -20182 -20950 -21460 -21972 -19181 -18421    0      0   ...
%--                  0     0      0       0      0      0      0      0      0   ... 
%--                  0];%--輸出結果原始值（從波形中得來）。
            
% xilsimtheta_du=xilsimtheta_r/8192/PI*180; %--輸出結果的度數表示。     
% abs_xilsimtheta_du=abs(xilsimtheta_du);
% abs_thetaout_du=abs(thetaout_du);%--理論值的絕對值。
% det_du=abs(abs_xilsimtheta_du-abs_thetaout_xil_du);%--仿真結果（FFT輸出相位）的絕對值誤差。

%--圖形化輸出。
subplot(2,1,1);
stem(indata,'*-r');
title('輸入數據時域值','color','b','FontSize',8);

subplot(2,1,2);
stem(abs(fftout_xil),'*-r');
title('FFT輸出理論值','color','b','FontSize',8);