编写一维快速傅里叶变换算法,并与matlab的fft函数进行对比
代码:
function [RET] = myfft(V)
%% 填充n=2^l
n = size(V,1);
originaln = n;
l = floor(log2(n));
VV = V;
if( power(2,l) ~= n)
l = floor(log2(n)) + 1;
n = power(2,l);
VV = zeros(n,1);
VV(1:size(V,1)) = V;
end
%% calculate inverse bit order
ibo = int16([1:n]');
for loop=1:n
ibo(loop) = bitinverse( loop-1, l );
end
%% assign inverse bit order
VV = VV(ibo+1);
%% 计算 FFT
Q = zeros(size(VV));
% 层数
for loop=1:l
%扫描
for index = 0:n-1
gap = power(2,loop);
block = floor(index / gap);
bindex = mod(index, gap/2);
half = mod(floor(index*2/gap),2);
W = exp(-i*2*pi*bindex/gap);
if( half==0)
Q(index+1) = VV(block*gap+bindex+1) + W.*VV(block*gap+bindex+1+gap/2);
else
Q(index+1) = VV(block*gap+bindex+1) - W.*VV(block*gap+bindex+1+gap/2);
end
end
VV=Q;
end
% return RET;
RET = VV(1:originaln);
end
%%% 比特倒序函数
function [ib] = bitinverse(bits, length)
ib=0;
for loop = 0:length-1
index = loop;
bitval = bitand( bits, bitshift(1,index));
bitval = bitshift( bitval, -index);
invindex = length-1-loop;
ib = bitor(ib, bitshift(bitval, invindex));
end
% returen ib
end