[i=s] 本帖最后由 fanyan861110 于 2011-1-26 15:33 编辑

基于FPGA实现FFT算法的思路，利用Matlab完成算法各级输出仿真，并且模拟了定点数计算，考虑有限字长。
对FFT算法的实现有真正的理论价值。


时域单周期采样序列图.JPG


理论摸值图.JPG


自定义函数simfft计算的模值图.JPG

Matlab代码如下
%----------------------------------------------------------------
% 作者：樊彦
% 功能：时间抽取POINTS(2^n)点基-2固定缩小比例(2/级)fft的Matlab验证程序
% 时间：2010年12月20日
% 说明：本代码主要用来分析，此种FFT算法中相位误差的来源
%       主函数主要用来产生仿真所需的输入数据，以及将理论输出数据(无截断)变换成
%   FPGA软件显示的数据格式，以便于进行比较。主函数调用了simfft函数。
% 版本：1.0
%----------------------------------------------------------------

function fft_theroy_test  %建立主函数
clc;
clear all;
close all;

PI=3.14159265357;
DWIDTH=8;   %数据位宽
%--下面的程序用于产生仿真所需的输入信号(正弦信号)
%  |
%  |          *
%  |      *       *
%  |   *             *
%  | *                 *
%--|*-------------------*---------------------→                  
%  |0 1 ...              *                 *31
%  |                       *             *
%  |                          *       *
%  |                              *

%-----------------------------------------------------------------
%采样点数压缩于一个模拟周期
% f=31.25;Um=1;nT=1;  %输入信号频率，振幅，显示周期数
% POINTS=32;T=1/f; %POINTS为一个周期的采样点数，T为信号周期
% dt=T/POINTS;     %dt为采样间隔
% Fs=POINTS*f;     %fs为采样频率    
% N=0:nT*POINTS-1;     %nt个周期的样点n从0开始
% tN=N*dt;        %确定时间序列样点在时间轴上的位置
% indata_signal=Um*cos(2*PI*f*tN);
%-----------------------------------------------------------------
%将输入信号恰好完整从满整个点数
Fs=1000;Um=100;
POINTS=32;
MULT_F=3;
OMIGA=2*PI*MULT_F/(POINTS-1);
tN=0:POINTS-1;
N=tN;
indata_signal=Um*cos(OMIGA*N);
%-----------------------------------------------------------------
%利用正弦型号的2规定：1→Fs=m*f0,m>2的正整数，2→N=k*m,k>1的整数
%推导：sin(2*pi*f0*t)→采样:sin(2*pi*f0/Fs*n),n=0,1,...,N-1
%→代入规定数据得：sin(2*pi**n/m)=sin(2*pi*n*k/N)
% POINTS=32;m=4;k=POINTS/m;
% f=250;Um=100;Fs=m*f;
% MIGA=2*PI*k/POINTS;
% tN=0:POINTS-1;
% N=tN;
% indata_signal=Um*cos(OMIGA*N);
%-----------------------------------------------------------------
%图形-数据可视化
set(figure(1),'Name','时域单周期采样序列图')
stem(tN,indata_signal,'.');grid on;
axis([min(tN) max(tN)  -1.5*max(abs(indata_signal)) 1.5*max(abs(indata_signal))]);  %限定横坐标和纵坐标的显示范围
text(tN,indata_signal,num2cell(N));
title('时域单周期采样序列图');
%semilogy(w,Hr);grid on;   %x轴线性,y轴对数
ylabel(' Amplitude value ');
xlabel(' Discrete Time (s)');
%-----------------------------------------------------------------
%归一化后对数据进行截断取整→否则随着运算级数的增加，位宽扩展成指数增加
CONSTANT=0;
indata_unify=(indata_signal+CONSTANT)/max(indata_signal+CONSTANT);  %归一化数据-1～１
%数据截断取整到有符号8位2进制表示的十进制范围(bit7-bit0,其中bit7为符号位)-128～127
indata=round(indata_unify*2^(DWIDTH-1));        

for i=1:POINTS         %产生仿真波形中的输入数组(要用补码来表示)补码表示范围0～255
    if indata(i)<0
        indata_comp(i)=2^DWIDTH+indata(i);  %取负数的补码
    elseif 2^(DWIDTH-1)==indata(i)   %避免对数值“+128”归一化时的错误
        indata_comp(i)=indata(i)-1;  %用最大整数“+127”代替
    else
        indata_comp(i)=indata(i);
    end
end

%-----------------------------------------------------------------
%幅频响应&&相频响应--理论输出
fftout_mat=fftshift(fft(indata));  %Matlab计算出的FFT理论值
fftout_xil=floor(fftout_mat/2^log2(POINTS));  %因为在设计中为了防止每级蝶形运算的数据溢出，所以每级输出数据缩小了2倍
absout_xil=abs(fftout_xil);           %FFT结果的理论摸值
%absout_xil_dB=20*log10(absout_xil);   %理论摸值(dB)
thetaout=angle(fftout_xil);           %FFT结果的理论相位值(弧度表示)
thetaout=unwrap(thetaout);            % 解卷绕
thetaout=thetaout/PI*180;              %度数表示的FFT相位理论输出结果
thetaout_xil=thetaout*2^(DWIDTH-3);  %按Xilinx的CORDIC核中定义的数据格式缩放后的数值(波形中的理论数值，便于比较用)

sN=length(fftout_xil);
F_axis=(-Fs/2:Fs/sN:Fs/2-Fs/sN);           %FFT程序的实际频率
%-----------------------------------------------------------------
%理论数值图像可视化
set(figure(2),'Name','幅频响应&&相频响应图')
subplot(211)
stem(F_axis,absout_xil,'.');grid on;               %画归一化的幅频率响应曲线
%axis([-1.1*max(abs(F_axis)) 1.1*max(abs(F_axis))  -1.1*max(abs(absout_xil)) 1.1*max(abs(absout_xil))]);  %限定横坐标和纵坐标的显示范围
text(F_axis,absout_xil,num2cell(N));
title('理论摸值图');
%semilogy(w,Hr);grid on;   %x轴线性,y轴对数
ylabel('Amplitude Freq. Res.');
xlabel('Digital Freq (Hz)');
subplot(212)
stem(F_axis,thetaout,'.');grid on;             %画相频响应虚线
text(F_axis,thetaout,num2cell(N));
title('理论相位值图');
ylabel(' Phase Freq. Res.(Degree)');
xlabel(' Digital Freq (Hz)');
%-----------------------------------------------------------------
%结果输出，调用自定义函数simfft（转成列向量便于观察）
[g0out,g1out,g2out,g3out,g4out]=simfft(indata.')  %输出FFT运算每一级的运算结果。
%-----------------------------------------------------------------
%幅频响应&&相频响应--理论输出
fftout_simfft_xil=fftshift(g4out);  %调用自定义函数simfft计算出的FFT值
%fftout_simfft_xil=floor(fftout_simfft_xil/2^log2(POINTS));  %因为在设计中为了防止每级蝶形运算的数据溢出，所以每级输出数据缩小了2倍
absout_simfft_xil=abs(fftout_simfft_xil);  %FFT结果的摸值
%absout_xil_dB=20*log10(absout_xil);      %摸值(dB)
thetaout_simfft=angle(fftout_simfft_xil);        %FFT结果的相位值(弧度表示)
thetaout_simfft=unwrap(thetaout_simfft);         % 解卷绕
thetaout_simfft=thetaout_simfft/PI*180;            %度数表示的FFT相位输出结果
%thetaout_simfft_xil=thetaout_simfft*2^(DWIDTH-3);  %按Xilinx的CORDIC核中定义的数据格式缩放后的数值(波形中的理论数值，便于比较用)

sN=length(fftout_simfft_xil);
F_axis=(-Fs/2:Fs/sN:Fs/2-Fs/sN);           %FFT程序的实际频率
%-----------------------------------------------------------------
%数值图像可视化
set(figure(3),'Name','幅频响应&&相频响应图')
subplot(211)
stem(F_axis,absout_simfft_xil,'.');grid on;               %画归一化的幅频率响应曲线
%axis([-1.1*max(abs(F_axis)) 1.1*max(abs(F_axis))  -1.1*max(abs(absout_xil)) 1.1*max(abs(absout_xil))]);  %限定横坐标和纵坐标的显示范围
text(F_axis,absout_simfft_xil,num2cell(N));
title('自定义函数simfft计算的摸值图');
%semilogy(w,Hr);grid on;   %x轴线性,y轴对数
ylabel('Amplitude Freq. Res.');
xlabel('Digital Freq (Hz)');
subplot(212)
stem(F_axis,thetaout_simfft,'.');grid on;             %画相频响应虚线
text(F_axis,thetaout_simfft,num2cell(N));
title('自定义函数simfft计算的相位值图');
ylabel(' Phase Freq. Res.(Degree)');
xlabel(' Digital Freq (Hz)');
%-----------------------------------------------------------------

%---------------------------------------------------------------%
%--      32点的基-2FFT每级数据输出函数（32点→log2(32)=5级）    --%
%--            作者：樊彦                                      --%
%--            时间：2010年12月21日                            --%
%---------------------------------------------------------------%
function [g0out,g1out,g2out,g3out,g4out]=simfft(indata);
%--输出数据：g1out...g5out都是8位有符号数(-128～+127);
%--输入数据：indata也是8位有符号数(-128～+127);
PI=3.14159265357;
DataWidth=8;   %旋转因子的精度（数据宽度）
PointDot=32;   %时间抽取基-2FFT的点数

WnDot=PointDot/2;  %旋转因子数组元素个数
%最高级（M-1）的旋转因子
%wn=cos(2*PI/PointDot*r)-j*sin(2*PI/PointDot*r),r=0,1,...,PointDot/2-1
%最高级的旋转因子包含了各级的旋转因子，并呈规律的连续取值
%故旋转因子矩阵表示如下：
%注意：旋转因子需要扩大方便后续计算
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));  %四舍五入取整（旋转因子扩大了2^(DataWidth-1)倍）
wn_i_t=-round(sin_ary*2^(DataWidth-1)); %此处要加“-”号，以便于后面处理

for i=1:WnDot  %产生wn_r和wn_i数组（不用补码来表示）
    if 2^(DataWidth-1)==wn_r_t(i)  %避免对“128”归一化时的错误
        wn_r(i)=wn_r_t(i)-1;
    else
        wn_r(i)=wn_r_t(i);
    end
    if 2^(DataWidth-1)==wn_i_t(i)
        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 2^(DataWidth-1)==indata(i)  %避免对“128”归一化时的错误
        indata_comp(i)=indata(i)-1;
%     elseif indata(i)<0
%         indata_comp(i)=2^DataWidth+indata(i);       
    else
        indata_comp(i)=indata(i);
    end
end
indata_comp=indata_comp.';   %转成列向量便于观察
indata_bitrev=bitrevorder(indata_comp);  %输入数据位反转（位码倒置）注意：输入数据必须为2^n

%-----------------------------------
%以“组”和“级”的概念进行循环运算输出
%“组”和“级”的关系：
%第m级的组数=PointDot/2^(m+1),m=0,1,...,M-1(log2(PointDot)-1)
%旋转因子的分布于“级”的关系如下：
%第m级的旋转因子矩阵：规律如下
%m=0级，W（r,2）,r=0
%m=1级，W（r,4）,r=0,1
%m=2级，W（r,8）,r=0,1,2,3
%...
%m=M-1级，W（r,PointDot）,r=0,1,...,PointDot/2-1
%归结起来得：
%第m级，W(r,2^(m+1)),r=0,1,...,2^m-1
%-----------------------------------
%旋转因子位码的确定：
%首先，在调用butfly时,输入参数wn的位码需要把理论旋转因子W(r,2^(m+1))的2^(m+1)换算成PointDot，
%即放大PointDot/2^(m+1)=WnDot/2^m倍。
%然后，r值等价变化，最后根据每组蝶形的个数确定旋转因子矩阵的位码，别忘记在确定的位码上+1！！
%-----------------------------------
%第m级每蝶形单元两输入数据位码的关系：
%|p-q|=2^m
%-----------------------------------
%运算的实现如下：分外,内循环实现
%外：每组蝶形个数（1：WnDot(每级蝶形总数)/组数）由于WnDot/PointDot*2^(m+1)=2^m，所以化简为(1:2^m)
%内：每组蝶形需要的循环的次数（1:2*每组蝶形个数(每组蝶形的点数):PointDot(总点数)）
%输入输出数据按位码顺序排列的取法：q,q+2^m
%注意：~！！！！
%第0级蝶形输入为时序序列x(n)的位码倒序值，即indata_bitrev
%第1级蝶形输入为第0级输出g0out的同位码值，即g0out
%第2级蝶形输入为第1级输出g1out的同位码值，即g1out
%...
%以此类推
%-----------------------------------
%--第0级输出
for i=1:2:PointDot  %内循环（外循环略，因为每组蝶形数位1）
 [g0out(i),g0out(i+1)]=butfly(indata_bitrev(i),indata_bitrev(i+1),wn(1));
end
g0out=g0out.';
%--第1级输出
for i=1:2  %外循环:每组蝶形个数
    for j=1:4:PointDot  %内循环：每组蝶形需要的循环的次数
        [g1out(i-1+j),g1out(i-1+j+2)]=butfly(g0out(i-1+j),g0out(i-1+j+2),wn(8*(i-1)+1));
    end
end
g1out=g1out.';
%--第2级输出
for i=1:4
    for j=1:8:PointDot
        [g2out(i-1+j),g2out(i-1+j+4)]=butfly(g1out(i-1+j),g1out(i-1+j+4),wn(4*(i-1)+1));
    end
end
g2out=g2out.';
%--第3级输出
for i=1:8
    for j=1:16:PointDot
        [g3out(i-1+j),g3out(i-1+j+8)]=butfly(g2out(i-1+j),g2out(i-1+j+8),wn(2*(i-1)+1));
    end
end
g3out=g3out.';
%--第4级输出
for i=1:16 %外循环（内循环略）
    [g4out(i),g4out(i+16)]=butfly(g3out(i),g3out(i+16),wn(i));
end
g4out=g4out.';

%---------------------------------------------------------------%
%--          单个蝶形运算函数butfly                            --%
%--            作者：樊彦                                      --%
%--            时间：2010年12月21日                            --%
%---------------------------------------------------------------%
function [xk1,xk2]=butfly(x1,x2,wn);
%输入：x1,x2,wn都是8位有符号数(-128～127),复数形式
%输出：xk1,xk2也是8位有符号数(-128～127),复数形式
%注意：wn应是先放大2^(DataWidth-1)倍后作为输入数据输入
DataWidth=8;
PointDot=32;
xr1=real(x1);
xi1=imag(x1);
xr2=real(x2);
xi2=imag(x2);
wr=real(wn);
wi=imag(wn);

%蝶形运算求积运算（均是复数运算）：
%(xr2+j*xi2)*(wr+j*wi)=midsubout+j*midaddout=(xr2*wr-xi2*wi)+j*(wr2*wi+xi2*wr)
mul1out=xr2*wr;   %4个乘法器的输出
mul2out=xi2*wi;
mul3out=xr2*wi;
mul4out=xi2*wr;
midsubout=mul1out-mul2out;  %中间的加法器和减法器
midaddout=mul3out+mul4out;
%蝶形运算求和运算（均是复数运算）：
%xk1=x1+wn*x2=xkr1+j*xki1=(xr1+midsubout)+j*(xi1+midaddout); 
%xk2=x1-wn*x2=xkr2+j*xki2=(xr1-midsubout)+j*(xi1-midaddout);
%注意：输入的wn事先已经扩大了2^(DataWidth-1)倍
xkr1=floor(xr1/2)+floor(midsubout/2^DataWidth);  %最后输出前的加法和减法器
xkr2=floor(xr1/2)-floor(midsubout/2^DataWidth);  %为防止数据溢出输出数据固定缩小2倍
xki1=floor(xi1/2)+floor(midaddout/2^DataWidth);  
xki2=floor(xi1/2)-floor(midaddout/2^DataWidth);

xk1=xkr1+j*xki1;  %组合成复数后输出
xk2=xkr2+j*xki2;