基-4FFT实施

Question

的倍频radix-4 FFT代码下面工作正常值逐案。

$ octave fft4.m 

ans = 1.4198e-015 

但是，如果我去掉循环代码我收到以下错误

$ octave fft4.m 

error: `stage' undefined near line 48 column 68 
error: evaluating argument list element number 1 
error: evaluating argument list element number 2 
error: called from: 
error: r4fftN at line 48, column 22 
error: c:\Users\david\Documents\Visual Studio 2010\Projects\mv_fft\fft4.m at line 80, column 7 

the "error" refers to a line the in fft function code which otherwise works correctly when xp is not set by a loop ... very strange. 


    function Z = radix4bfly(x,segment,stageFlag,W) 
     % For the last stage of a radix-4 FFT all the ABCD multiplers are 1. 
     % Use the stageFlag variable to indicate the last stage 
     % stageFlag = 0 indicates last FFT stage, set to 1 otherwise 

     % Initialize variables and scale to 1/4 
     a=x(1)*.25; 
     b=x(2)*.25; 
     c=x(3)*.25; 
     d=x(4)*.25; 

     % Radix-4 Algorithm 
     A=a+b+c+d; 
     B=(a-b+c-d)*W(2*segment*stageFlag + 1); 
     C=(a-b*j-c+d*j)*W(segment*stageFlag + 1); 
     D=(a+b*j-c-d*j)*W(3*segment*stageFlag + 1); 

     % assemble output 
     Z = [A B C D]; 
    end % radix4bfly() 


    % radix-4 DIF FFT, input signal must be floating point, real or complex 
    % 
    function S = r4fftN(s) 

     % Initialize variables and signals: length of input signal is a power of 4 
     N = length(s); 
     M = log2(N)/2; 

     % Initialize variables for floating point sim 
     W=exp(-j*2*pi*(0:N-1)/N); 
     S = complex(zeros(1,N)); 
     sTemp = complex(zeros(1,N)); 

     % FFT algorithm 
     % Calculate butterflies for first M-1 stages 
     sTemp = s; 
     for stage = 0:M-2 
     for n=1:N/4 
      S((1:4)+(n-1)*4) = radix4bfly(sTemp(n:N/4:end), floor((n-1)/(4^stage)) *(4^stage), 1, W); 
     end 
     sTemp = S; 
     end 

     % Calculate butterflies for last stage 
     for n=1:N/4 
     S((1:4)+(n-1)*4) = radix4bfly(sTemp(n:N/4:end), floor((n-1)/(4^stage)) * (4^ 
    stage), 0, W); 
     end 
     sTemp = S; 

     % Rescale the final output 
     S = S*N; 
    end % r4fftN(s) 


    % test FFT code 
    % 
    xp = 2; 

    % ERROR if I uncomment loop! 

    %for xp=1:8 
     N = 4^xp; % must be power of: 4 16 64 256 1024 4086 .... 
     x = 2*pi/N * (0:N-1); 
     x = cos(x); 
     Y_ref = fft(x); 
     Y = r4fftN(x); 
     Y = digitrevorder(Y,4); 
     %Y = bitrevorder(Y,4); 
     abs(sum(Y_ref-Y)) % compare fft4 to built-in fft 
    %end 

Answer 1

的问题是环开往指数XP作为fft4代码假定至少2个由2应该开始基数-4蝴蝶阶段

对不起人:(

-David

Answer 2

请在下面找到一个完全工作的Matlab实现的基数-4抽取频率为的FFT算法。我还提供了复杂矩阵乘法和加法的总体操作计数。可以确实证明，每个基4蝶形运算包括复数乘法和8复数加法运算。由于有log_4(N) = log_2(N)/2阶段，并且每个阶段涉及N/4蝴蝶，以使操作数是

complex multiplications = (3/8) * N * log2(N) 
complex additions  = N * log2(N) 

下面是代码：

% --- Radix-2 Decimation In Frequency - Iterative approach 

clear all 
close all 
clc 

% --- N should be a power of 4 
N = 1024; 

% x = randn(1, N); 
x = zeros(1, N); 
x(1 : 10) = 1; 
xoriginal = x; 
xhat = zeros(1, N); 

numStages = log2(N)/2; 

W = exp(-1i * 2 * pi * (0 : N - 1)/N); 
omegaa = exp(-1i * 2 * pi/N); 

mulCount = 0; 
sumCount = 0; 

M = N/4; 
for p = 1 : numStages; 
    for index = 0 : (N/(4^(p - 1))) : (N - 1); 
     for n = 0 : M - 1; 
      a = x(n + index + 1) +  x(n + index + M + 1) + x(n + index + 2 * M + 1) +  x(n + index + 3 * M + 1); 
      b = (x(n + index + 1) -  x(n + index + M + 1) + x(n + index + 2 * M + 1) -  x(n + index + 3 * M + 1)) .* omegaa^(2 * (4^(p - 1) * n)); 
      c = (x(n + index + 1) - 1i * x(n + index + M + 1) - x(n + index + 2 * M + 1) + 1i * x(n + index + 3 * M + 1)) .* omegaa^(1 * (4^(p - 1) * n)); 
      d = (x(n + index + 1) + 1i * x(n + index + M + 1) - x(n + index + 2 * M + 1) - 1i * x(n + index + 3 * M + 1)) .* omegaa^(3 * (4^(p - 1) * n)); 
      x(n + 1 + index) = a; 
      x(n + M + 1 + index) = b; 
      x(n + 2 * M + 1 + index) = c; 
      x(n + 3 * M + 1 + index) = d; 
      mulCount = mulCount + 3; 
      sumCount = sumCount + 8; 
     end; 
    end; 
    M = M/4; 
end 

xhat = bitrevorder(x); 

tic 
xhatcheck = fft(xoriginal); 
timeFFTW = toc; 

rms = 100 * sqrt(sum(sum(abs(xhat - xhatcheck).^2))/sum(sum(abs(xhat).^2))); 

fprintf('Theoretical multiplications count \t = %i; \t Actual multiplications count \t = %i\n', ... 
     (3/8) * N * log2(N), mulCount); 
fprintf('Theoretical additions count \t\t = %i; \t Actual additions count \t\t = %i\n\n', ... 
     N * log2(N), sumCount); 
fprintf('Root mean square with FFTW implementation = %.10e\n', rms);