结果

代码：

clear,clc,close all;
%自己编程实现时间抽取基2fft算法
f=50;%频率
fs=10000;%采样频率
Ts=1/fs;%采样间隔
N=128;
n=1:N;
data1=sin(2*pi*f*n*Ts);
% figure;
% plot(n*Ts,data1);title('原始信号');
data1;%原始信号
%%将数据长度变为2的幂次方
if(mod(log2(N),1)~=0)
    temp=ceil(log2(N));
    data1=[data1,zeros(1,2^temp-N)];
    N=N+2^temp-N;
end
%%对数据进行重新排列
for index1=1:N
    index2(index1)=bin2dec(fliplr(dec2bin(index1-1,log2(N))));
    data_p(index1)=data1(index2(index1)+1);
end
%%进行运算
L1=log2(N);%L1级运算
data=zeros(N,L1);
data(:,1)=data_p';


%m级运算
for m=0:L1-1
    %k组,k=N/(2^(m+1))
    for k=1:N/(2^(m+1))
        
        %计算鲽形对
        %%%%为什么这里用exp(j*2*pi*r/2^(m+1))才正确，按理来说应该用exp(-j*2*pi*r/2^(m+1))？
        %用exp(j*2*pi*r/2^(m+1))和fft计算结果一致
        for r=0:2^(m)-1
            w(r+1)=exp(j*2*pi*r/2^(m+1));
        end
        
        %每组n个鲽形单元，n=2^m
        for n=1:2^(m+1)
            if(n<=2^m)
                flag=1;
                data((k-1)*(2^(m+1))+n,m+2)=data((k-1)*(2^(m+1))+n,m+1)+flag*w(mod((k-1)*(2^(m+1))+n-1,2^m)+1)*data((k-1)*(2^(m+1))+n+2^m,m+1);
            end
            if(n>2^m)
                flag=-1;
                data((k-1)*(2^(m+1))+n,m+2)=data((k-1)*(2^(m+1))+n-2^m,m+1)+flag*w(mod((k-1)*(2^(m+1))+n-1,2^m)+1)*data((k-1)*(2^(m+1))+n,m+1);
            end
            
        end
    end
end
fft1=fft(data1)';
figure;
subplot(211)
plot(abs(data(:,end)))
title('自己写的代码计算的');
subplot(212)
plot(abs(fft1));
title('fft计算的');
answer1=data(:,end)-fft1;

1024点的结果：