Back History
数字图像处理——作业6

作业内容：编写一维快速傅里叶变换算法，并与matlab的fft函数进行对比

代码：
{{{{<nowiki>
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
</nowiki>
}}}}