2N点实序列的FFT, MATLAB CODE

2N点实序列的FFT

• 空间

function X = fft_2N(x)

% 计算长度为2N的实序列的FFT

% John G.Proakis,数字信号处理,第四版

% Define variables:

%   x       --- Input real vector of length 2*N(power of 2)

%   X       --- Output vector,FFT result

% Record of revisions:

%       Date          Programmer        Description of change

%       ====          ==========        =====================

%      2014/0421        GOU DH           Original code

%

% test

clear,clc

x = randn(64,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x = x(:);

nx = length(x);

t = reshape(x,2,nx/2);

x1 = t(1,:)';

x2 = t(2,:)';            % 时域抽取,奇偶分离

[X1,X2] = fft_2_N(x1,x2,nx/2);

k = [0:nx/2-1]';

W2Nk = exp(-j*2*pi*k/(nx));

T = X2.*W2Nk;

G1 = X1 + T;

G2 = X1 - T;

X = [G1;G2];

% test

sum(abs(X - fft(x,nx)))

return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Subfunction

function [X1,X2] = fft_2_N(x1,x2,N)

% 计算2个N点实序列的FFT

if ~isreal(x1) | ~isreal(x2)

    return;

end

x = x1 + j*x2;

X = fft_dif2(x,-1);

T = [X(1); flipud(X(2:end))];

T = conj(T);

X1 = (X + T)/2;

X2 = (X - T)/(2*j);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Subfunction

function [X]=fft_dif2(x,f)

% 频域抽取基2-FFT(基2-DIF-FFT)

nx = length(x);

M = fix(log2(nx));

if 2^M < nx, M = M + 1; end       % 级数M

x = [x(:);zeros(2^M-nx,1)];

N = length(x);

%

C = j*pi*f;

for m = 0:1:M-1          % 级循环序号

    B = 2^(M-1-m);       % 蝶形因子序号间距=蝶形因子个数=2^(M-1-m)

    for r = 0:1:B-1      % 蝶形因子循环序号

        Wnp = exp(C*r/B);

        for k = r:2*B:N-1

            tmp1 = x(k+1);

            tmp2 = x(k+B+1);

            x(k+1) = tmp1+tmp2;

            x(k+B+1) = (tmp1-tmp2)*Wnp;

        end

    end

end

% 倒位序

% X = x;

LH = N/2;

J = LH;

N1 = N-2;

for I = 1:1:N1

    if I < J

        T = x(I+1);

        x(I+1) = x(J+1);

        x(J+1) = T;

    end

    K = LH;

    while J >= K

        J = J-K;

        K = K/2;

    end

    J = J+K;

end

if f == 1

    x = x/N;

end

X = x;