function x1 = FFT(x)
%FFT 自己实现的快速傅里叶变换程序
% Input:Vector x on time domain, Output: Vector y on frequency domain.
% 2022.11.17 复旦大学 张叶昊
N = length(x);
n = ceil(log2(N));
len = 2^n;
%% 补零操作,这里采用的是一个复制,一定程度上连累了速度。
x1 = zeros(len,1);
% y1 = zeros(len,1);
x1(1:N) = x;
% y1(1:N) = y;
%% 倒位序的实现
% 我们希望找到x(k)被挪到了X(K)的K与k的对应关系,以及实现算法
% 给定一个十进制数k,k有一个二进制表示,以五位为例,不妨记为(abcde),
% 那么K的二进制表示就是(edcba),然后将它换回十进制。
% 在k的数组里面是00000,00001,00010,00011这样进位的,是从右往左进位的。
% 在K的数组里面是00000,10000,01000,11000这样进位的,是从左往右进位的。
% 下面这个叫Radar算法
half_len=len/2;%数组半长
j = half_len;
for i=1:N-2%这里实现了奇偶前后分开排序
%比较前半部分序数【0 1 2 3 】,对每对中的后一个偶数进行交换1换4,3换6
if i<j
t=x1(j+1);
x1(j+1)=x1(i+1);%交换
x1(i+1)=t;
end
%求下一个倒序数
k =half_len;
while(j>=k)%j为下一个倒序数,比较100的最高位1,如果为1,置0
j=j-k;
k=k/2;%变为010,准备比较第二位的1,循环
% j=j+k;%找到为0 的一位,补成1,j就是下一个倒序数
end
% if j<k
j=j+k;%找到为0 的一位,补成1,j就是下一个倒序数
% end
end
%% 下一步是蝶形运算
% 先把所有要用到的矩阵都存储起来
% m是运算级数,q是蝶形运算的组数。
% coef = -ones(len,1);
% coef(1:half_len-1) = exp(-2*1j*pi *((1:half_len-1)./ len));
% coef(half_len+1:end) = -coef(1:half_len);
coef = ones(half_len,1);
coef(2:half_len) = exp(-2*1j*pi *((1:half_len-1)./ len));
%% 蝶形算法简介:
% 比如一个16bit,已经置换好的情况下,
% 我们将进行蝶形运算时使用同一个旋转因子的称为一组,
% 其中不用乘旋转因子的称为第一部分,要乘蝶形因子的成为第二部分。
% 我们现在一共有N = 2^n个数,在step = k,也就是第k级运算的时候
% 我们有group_num = 2^(k-1)组,其中每组里面有两个部分,所以一共有2^k个部分
% 每个部分里面有2^(n-k)个元素,每个部分内部元素相隔ingroup_dist = 2^k
% 下面给一个例子,对于16比特而言,是step = 1,只有一组,ingroup_dist = 2
% 13579111315 是第一部分,246810121416是第二部分
% 进行第一次操作的时候,第一组不动,第二组要乘旋转因子W_{ingroup_dist}^0
% 然后蝶形运算一次
% step = 2
% group_num = 2^(2-1) = 2
% (15913,371115) 是第一组,(261014,481216)是第二组。
% ingroup_dist = 4
% 每一组的第二部分要乘旋转因子W_{ingroup_dist}^{group_idx}
% 然后进行蝶形运算
% step =3
% group_num = 2^(3-1) = 4
% (19,513)(26,1014)(37,1115)(48,1216)有四组
% 每一组的第二部分要乘旋转因子W_{ingroup_dist}^{group_idx}
% 然后进行蝶形运算
% 最后要注意,当step = k的时候,N = ingroup_dist*2^(n-k),所以
% W_{ingroup_dist}^{group_idx} = W_{N}^{group_idx*2^(n-k)}
% = W_{N}^{group_idx*part_vol} = coef(group_idx*part_vol)
% 回顾变量:ingroup_dist,同一组的一个部分里面两个元素索引差
% group_vol = 一个组里面两个部分加起来的总容量
for step = 1:n
ingroup_dist = 2^step;
group_num = ingroup_dist/2;%2^(step-1)
group_vol = 2^(n-step);
% part_dist = group_num;% 同一组内相互作用的两个元素的差
x1 = reshape(x1,ingroup_dist,group_vol);
%对coef的向量化并没有起到非常好的效果。
% coef = reshape(coef,group_vol,group_num);
temp = x1(group_num+1:end,:).*coef(1:group_vol:end);
% temp = x1(group_num+1:end,:).*coef(1,:).';
x1(group_num+1:end,:) = x1(1:group_num,:)-temp;
x1(1:group_num,:) = x1(1:group_num,:)+temp;
% 本来应该写成下面这个样子,但是从数值实验上面看这一段是真的慢
% 慢的原因是索引太慢了,所以干脆改掉,通过重排,改成上面那个样子
% 并且引入向量化运算,速度提升了很多。
% for group_idx= 1:group_num
% I = group_idx:ingroup_dist:N-ingroup_dist+group_idx;
% II = I+group_num;%实际上加的是part_dist,这里就不多引入变量了。
% temp = coef((group_idx-1)*group_vol+1)*x1(II);
% x1(II) = x1(I)-temp;
% x1(I) = x1(I)+temp;
% end
end
end