www.gusucode.com > GPS仿真Matlab编程源码程序 > GPS仿真Matlab编程源码程序/myfft.m

    % x=ones(1,128); %输入的信号，自己可以改变 
clear
clc
Nt=200;
t=1:Nt;
fc1=1.25;
% fc1=4.30975427e6;
fs=7.5;
x=exp(-1i*2*pi*fc1/fs*t);%输入的信号，自己可以改变 
%整体运用原位计算 
m=nextpow2(x);N=2^m; % 求x的长度对应的2的最低幂次m 
if length(x)<N 
x=[x,zeros(1,N-length(x))]; % 若x的长度不是2的幂，补零到2的整数幂 
end 
nxd=bin2dec(fliplr(dec2bin([1:N]-1,m)))+1; % 求1:2^m数列序号的倒序 
y=x(nxd); % 将x倒序排列作为y的初始值 
for mm=1:m % 将DFT作m次基2分解,从左到右，对每次分解作DFT运算,共做m级蝶形运算，每一级都有2^(mm-1)个蝶形结 
Nz=2^mm;u=1; % 旋转因子u初始化为WN^0=1 
WN=exp(-1i*2*pi/Nz); % 本次分解的基本DFT因子WN=exp(-i*2*pi/Nz) 
for j=1:Nz/2 % 本次跨越间隔内的各次蝶形运算,在进行第mm级运算时需要2^(mm-1)个 蝶形 
for k=j:Nz:N % 本次蝶形运算的跨越间隔为Nz=2^mm 
kp=k+Nz/2; % 蝶形运算的两个因子对应单元下标的关系 
t=y(kp)*u; % 蝶形运算的乘积项 
y(kp)=y(k)-t; % 蝶形运算 
y(k)=y(k)+t; % 蝶形运算 
end 
u=u*WN; % 修改旋转因子,多乘一个基本DFT因子WN 
end 
end 
y;
y1=fft(x,N); %自己编的FFT跟直接调用的函数运算以后的结果进行对比