结果
代码:
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点的结果: