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