Main Content

高速フーリエ変換 (FFT) から固定小数点への変換

この例では、高速フーリエ変換 (FFT) アルゴリズムのテキストブック バージョンを固定小数点の MATLAB® コードに変換する方法を説明します。

以下のコードを実行して現在の状態を取得し、グローバル状態をリセットします。

FIPREF_STATE = get(fipref);
reset(fipref)

テキストブック FFT アルゴリズム

FFT は、時間領域から周波数領域への複素数の線形変換です。たとえば、2 つの正弦波の合計としてベクトルを作成し、それを FFT で変換すると、FFT 振幅プロットで周波数のピークを確認できます。

n = 64;                                     % Number of points
Fs = 4;                                     % Sampling frequency in Hz
t  = (0:(n-1))/Fs;                          % Time vector
f  = linspace(0,Fs,n);                      % Frequency vector
f0 = 0.2; f1 = 0.5;                         % Frequencies, in Hz
x0 = cos(2*pi*f0*t) + 0.55*cos(2*pi*f1*t);  % Time-domain signal
x0 = complex(x0);                           % The textbook algorithm requires the input to be complex
y0  = fft(x0);                              % Frequency-domain transformation fft() is a MATLAB built-in function

fi_fft_demo_ini_plot(t,x0,f,y0);            % Plot the results from fft and time-domain signal

Figure contains 2 axes objects. Axes object 1 contains an object of type line. This object represents x0. Axes object 2 contains an object of type line. This object represents abs(fft(x0)).

周波数プロットの 0.2 および 0.5 Hz でのピークは、それらの周波数での時間領域信号の 2 つの正弦波に対応します。

3.5 および 3.8 Hz の反射ピークに注目してください。この場合のように FFT への入力が実数値である場合、出力 y は共役対称になります。

y(k)=conj(y(n-k))

FFT のさまざまな実装があり、それぞれに独自のコストと利点があります。アプリケーションによっては、ここで指定したアルゴリズムではない別のアルゴリズムの方が適していることがあります。このアルゴリズムは、独自の調査を開始する方法の例として使用しています。

この例は、Charles Van Loan 著『Computational Frameworks for the Fast Fourier Transform』の「Algorithm 1.6.2」(45 ページ) に示されている時間間引きのユニット幅 FFT を使用しています。

疑似コードでのテキストブックのアルゴリズムは次のとおりです。

アルゴリズム 1.6.2。x が長さ nn=2t の複素数ベクトルの場合、次のアルゴリズムによって xFnx に書き換えられます。

x=Pnx

w=wn(long) (Van Loan のセクション 1.4.11 を参照)

for q=1:t

L=2q;r=n/L;L*=L/2;

for k=0:r-1

for j=0:L*-1

τ=w(L*-1+j)x(kL+j+L*)

x(kL+j+L*)=x(kL+j)-τ

x(kL+j)=x(kL+j)+τ

end

end

end

テキストブックのアルゴリズムは、ゼロから始まるインデックスを使用します。Fn は n 行 n 列のフーリエ変換行列、Pn は n 行 n 列のビット反転の置換行列、w は回転因子の複素数ベクトルです。回転因子 w は、以下のアルゴリズムによって計算される 1 の根となる複素数です。

function w = fi_radix2twiddles(n) 
%FI_RADIX2TWIDDLES  Twiddle factors for radix-2 FFT example.
%   W = FI_RADIX2TWIDDLES(N) computes the length N-1 vector W of
%   twiddle factors to be used in the FI_M_RADIX2FFT example code.
%
%   See also FI_RADIX2FFT_DEMO.

%   Reference:
%
%   Twiddle factors for Algorithm 1.6.2, p. 45, Charles Van Loan,
%   Computational Frameworks for the Fast Fourier Transform, SIAM,
%   Philadelphia, 1992.
%
%   Copyright 2003-2011 The MathWorks, Inc.
%     

t = log2(n);
if floor(t) ~= t
  error('N must be an exact power of two.');
end

w = zeros(n-1,1);
k = 1;
L = 2;
% Equation 1.4.11, p. 34
while L <= n
  theta = 2*pi/L;
  % Algorithm 1.4.1, p. 23
  for j = 0:(L/2 - 1)
    w(k) = complex(cos(j*theta),-sin(j*theta));
    k = k + 1;
  end
  L = L*2;
end
figure(gcf)
clf
w0 = fi_radix2twiddles(n);
polarplot(angle(w0),abs(w0),'o')
title('Twiddle Factors: Complex roots of unity')

Figure contains an axes object. The axes object contains an object of type line.

浮動小数点コードの確認

MATLAB® でアルゴリズムを実装するために、関数 fi_bitreverse を使用して入力シーケンスをビット反転します。それらを 0 ベースから 1 ベースに変換するために、インデックスに 1 を加えなければなりません。

function x = fi_m_radix2fft_algorithm1_6_2(x, w)
%FI_M_RADIX2FFT_ALGORITHM1_6_2  Radix-2 FFT example.
%   Y = FI_M_RADIX2FFT_ALGORITHM1_6_2(X, W) computes the radix-2 FFT of
%   input vector X with twiddle-factors W.  Input X is assumed to be
%   complex.
%
%   The length of vector X must be an exact power of two.
%   Twiddle-factors W are computed via
%      W = fi_radix2twiddles(N)
%   where N = length(X).
%
%   This version of the algorithm has no scaling before the stages.
%
%   See also FI_RADIX2FFT_DEMO, FI_M_RADIX2FFT_WITHSCALING.

%   Reference:
%     Charles Van Loan, Computational Frameworks for the Fast Fourier
%     Transform, SIAM, Philadelphia, 1992, Algorithm 1.6.2, p. 45.
%
%   Copyright 2004-2015 The MathWorks, Inc.

    n = length(x);  t = log2(n);
    x = fi_bitreverse(x,n);
    for q=1:t
        L = 2^q; r = n/L; L2 = L/2;
        for k = 0:(r-1)
            for j = 0:(L2-1)
                temp          = w(L2-1+j+1) * x(k*L+j+L2+1);
                x(k*L+j+L2+1) = x(k*L+j+1)  - temp;
                x(k*L+j+1)    = x(k*L+j+1)  + temp;
            end
        end
    end
end

可視化

MATLAB® にアルゴリズムが正しく実装されたことを確認するには、既知の信号を実行し、その結果と MATLAB® FFT 関数によって生成される結果を比較します。

下のプロットに示すように、誤差が MATLAB® 組み込み FFT 関数の許容誤差内にあるので、アルゴリズムが正しく実装されたことを確認できます。

y = fi_m_radix2fft_algorithm1_6_2(x0, w0);

fi_fft_demo_plot(real(x0),y,y0,Fs,'Double data', ...
    {'FFT Algorithm 1.6.2','Built-in FFT'});

型テーブルを使用するように関数を変換

アルゴリズムからデータ型を分離する方法は以下のとおりです。

  1. データ型定義のテーブルを作成します。

  2. そのテーブルのデータ型を使用するようアルゴリズム コードを変更します。

この例では、さまざまなファイルを作成することによって反復的な手順について説明します。実際には、同じファイルに反復的な変更を加えることができます。

元の型テーブル

変数のプロトタイプが元の型に設定された構造体を使用して、型テーブルを作成します。ベースライン型を使用して、初期変換が正しく行われたことを検証すると共に、プログラムによって浮動小数点型と固定小数点型の間で関数を切り替えます。インデックス変数は MATLAB® Coder™ によって自動的に整数に変換されるため、この table で型を指定する必要はありません。

プロトタイプ値を空白 ([ ]) として指定します。これはこの値ではなくデータ型を使用するためです。

function T = fi_m_radix2fft_original_types()
%FI_M_RADIX2FFT_ORIGINAL_TYPES Types Table Example

%   Copyright 2015 The MathWorks, Inc.

    T.x = double([]);
    T.w = double([]);
    T.n = double([]);

end

型認識アルゴリズム関数

関数への入力として型 table T を追加し、これを使用することでアルゴリズム本体を変更せずに変数を特定の型にキャストします。

function x = fi_m_radix2fft_algorithm1_6_2_typed(x, w, T)
%FI_M_RADIX2FFT_ORIGINAL_TYPED Radix-2 FFT example.
%   Y = FI_M_RADIX2FFT_ALGORITHM1_6_2_TYPED(X, W, T) computes the radix-2
%   FFT of input vector X with twiddle-factors W.  Input X is assumed to be
%   complex.
%
%   The length of vector X must be an exact power of two.
%   Twiddle-factors W are computed via
%      W = fi_radix2twiddles(N)
%   where N = length(X).
%
%   T is a types table to cast variables to a particular type, while keeping
%   the body of the algorithm unchanged.
%
%   This version of the algorithm has no scaling before the stages.
%
%   See also FI_RADIX2FFT_DEMO, FI_M_RADIX2FFT_WITHSCALING.
%
%   Reference:
%     Charles Van Loan, Computational Frameworks for the Fast Fourier
%     Transform, SIAM, Philadelphia, 1992, Algorithm 1.6.2, p. 45.

%   Copyright 2015 The MathWorks, Inc.
%
%#codegen

    n = length(x);
    t = log2(n);
    x = fi_bitreverse_typed(x,n,T);
    LL = cast(2.^(1:t),'like',T.n);
    rr = cast(n./LL,'like',T.n);
    LL2 = cast(LL./2,'like',T.n);
    for q=1:t
        L = LL(q);
        r = rr(q);
        L2 = LL2(q);
        for k = 0:(r-1)
            for j = 0:(L2-1)
                temp          = w(L2-1+j+1) * x(k*L+j+L2+1);
                x(k*L+j+L2+1) = x(k*L+j+1)  - temp;
                x(k*L+j+1)    = x(k*L+j+1)  + temp;
            end
        end
    end
end

型認識ビット反転関数

関数への入力として型 table T を追加し、これを使用することでアルゴリズム本体を変更せずに変数を特定の型にキャストします。

function x = fi_bitreverse_typed(x,n0,T)
%FI_BITREVERSE_TYPED  Bit-reverse the input.
%   X = FI_BITREVERSE_TYPED(x,n,T) bit-reverse the input sequence X, where
%   N=length(X).
%
%   T is a types table to cast variables to a particular type, while keeping
%   the body of the algorithm unchanged.
%
%   See also FI_RADIX2FFT_DEMO.

%   Copyright 2004-2015 The MathWorks, Inc.
%
%#codegen
n = cast(n0,'like',T.n);
nv2 = bitsra(n,1);
j = cast(1,'like',T.n);
for i = 1:(n-1)
  if i < j
    temp = x(j);
    x(j) = x(i);
    x(i) = temp;
  end
  k = nv2;
  while k < j
    j(:) = j-k;
    k = bitsra(k,1);
  end
  j(:) = j+k;
end

変更後の関数の検証

関数を変更するたびに、関数の結果がベースラインと引き続き一致することを検証します。型テーブルで元の型を使用したため、出力は同一になるはずです。これにより、アルゴリズムから型を正しく分離するための変換が行われたことが検証されます。

T1 = fi_m_radix2fft_original_types(); % Get original data types declared in table

x = cast(x0,'like',T1.x);
w = cast(w0,'like',T1.w);

y = fi_m_radix2fft_algorithm1_6_2_typed(x, w, T1);

fi_fft_demo_plot(real(x),y,y0,Fs,'Double data', ...
    {'FFT Algorithm 1.6.2','Built-in FFT'});

Figure contains 3 axes objects. Axes object 1 contains an object of type line. This object represents Double data. Axes object 2 contains 2 objects of type line. These objects represent FFT Algorithm 1.6.2, Built-in FFT. Axes object 3 contains an object of type line. This object represents abs(error).

固定小数点型 table の作成

変数のプロトタイプをもつ構造体を使用して固定小数点型 table を作成します。プロトタイプ値を空白 ([ ]) として指定します。これはこの値ではなくデータ型を使用するためです。

function T = fi_m_radix2fft_fixed_types()
%FI_M_RADIX2FFT_FIXED_TYPES Example function

%   Copyright 2015 The MathWorks, Inc.

    T.x = fi([],1,16,14); % Picked the following types to ensure that the
    T.w = fi([],1,16,14); % inputs have maximum precision and will not
                          % overflow
    T.n = int32([]);      % Picked int32 as n is an index

end

固定小数点問題の特定

ここで、入力データを固定小数点に変換してもアルゴリズムに問題がないかを確かめます。この第 1 段階では、fi コンストラクターを使用して、符号付き固定小数点データに対してすべての既定値を使用します。

T2 = fi_m_radix2fft_fixed_types(); % Get fixed point data types declared in table

x = cast(x0,'like',T2.x);
w = cast(w0,'like',T2.w);

固定小数点入力を使用して同じアルゴリズムを再実行します。

y  = fi_m_radix2fft_algorithm1_6_2_typed(x,w,T2);
fi_fft_demo_plot(real(x),y,y0,Fs,'Fixed-point data', ...
    {'Fixed-point FFT Algorithm 1.6.2','Built-in FFT'});

Figure contains 3 axes objects. Axes object 1 contains an object of type line. This object represents Double data. Axes object 2 contains 2 objects of type line. These objects represent FFT Algorithm 1.6.2, Built-in FFT. Axes object 3 contains an object of type line. This object represents abs(error).

固定小数点の FFT の振幅プロット (中央) は、組み込み FFT のプロットとは似ていないことに注目してください。誤差 (下部のプロット) は、丸めの誤差に予想される誤差より大きいので、オーバーフローが起こる可能性が高くなります。

最小/最大計測機能を使用したオーバーフローの特定

MATLAB® コードをインストルメント化するには、buildInstrumentedMex コマンドを使用して MATLAB® 関数から MEX 関数を作成します。buildInstrumentedMex への入力は fiaccel への入力と同じですが、buildInstrumentedMex には fi オブジェクトの制限がありません。buildInstrumentedMex の出力は、計測機能を挿入された MEX 関数です。この MEX 関数を実行すると、すべての名前付き変数および中間値について、シミュレーションされた最小値と最大値が記録されます。

生成される MEX 関数に名前を付けるには、'-o' オプションを使用します。'-o' オプションを使用しない場合は、元の MATLAB® 関数の名前の最後に '_mex' を付けた名前が MEX 関数に付けられます。MEX 関数に MATLAB® 関数と同じ名前を付けることもできますが、MEX 関数の方が MATLAB® 関数より優先されるため、同じ名前の MEX 関数が再生成されるか削除されてクリアされるまでは MATLAB® 関数に対する変更が実行されない点に注意する必要があります。

値の最大範囲を保持して潜在的なオーバーフローを特定できるように、スケーリングされた double データ型の入力を作成します。

function T = fi_m_radix2fft_scaled_fixed_types()
%FI_M_RADIX2FFT_SCALED_FIXED_TYPES Example function

%   Copyright 2015 The MathWorks, Inc.

    DT = 'ScaledDouble';                % Data type to be used for fi
                                        % constructor
    T.x = fi([],1,16,14,'DataType',DT); % Picked the following types to
    T.w = fi([],1,16,14,'DataType',DT); % ensure that the inputs have
                                        % maximum precision and will not
                                        % overflow
    T.n = int32([]);                    % Picked int32 as n is an index

end
T3 = fi_m_radix2fft_scaled_fixed_types(); % Get fixed point data types declared in table

x_scaled_double = cast(x0,'like',T3.x);
w_scaled_double = cast(w0,'like',T3.w);

buildInstrumentedMex fi_m_radix2fft_algorithm1_6_2_typed ...
    -o fft_instrumented -args {x_scaled_double w_scaled_double T3}

インストルメント化した MEX 関数を実行して、最小値/最大値を記録します。

y_scaled_double = fft_instrumented(x_scaled_double,w_scaled_double,T3);

計測結果を表示します。

showInstrumentationResults fft_instrumented

InstrumentationResults.png

計測結果から、変数 x に割り当てるときにオーバーフローが発生したことがわかります。

固定小数点問題に対処するためのアルゴリズムの修正

FFT の個別ビンの振幅は、多くても n 倍 (n は FFT の長さ) までしか大きくなりません。そのため、データを 1/n 倍スケーリングすることで、入力でのオーバーフローの発生を防ぐことができます。長さ n の FFT の第 1 ステージへの入力のみを 1/n 倍スケーリングすると、n^2 に比例した S/N 比が取得されます [Oppenheim & Schafer 1989、方程式 9.101]、[Welch 1969]。ただし、FFT の各ステージへの入力を 1/2 倍スケーリングすると、1/n の全体スケーリングが得られ、n に比例した S/N 比を生成できます [Oppenheim & Schafer 1989、方程式 9.105]、[Welch 1969]。

固定小数点で 1/2 倍スケーリングする効率的な方法は、データを右シフトすることです。これを行うには、右にビット シフトする算術関数 bitsra を使用します。FFT の各ステージをスケーリングし、インデックス変数の計算を最適化すると、アルゴリズムは以下のようになります。

function x = fi_m_radix2fft_withscaling_typed(x, w, T)
%FI_M_RADIX2FFT_WITHSCALING  Radix-2 FFT example with scaling at each stage.
%   Y = FI_M_RADIX2FFT_WITHSCALING_TYPED(X, W, T) computes the radix-2 FFT of
%   input vector X with twiddle-factors W with scaling by 1/2 at each stage.
%   Input X is assumed to be complex.
%
%   The length of vector X must be an exact power of two.
%   Twiddle-factors W are computed via
%      W = fi_radix2twiddles(N)
%   where N = length(X).
%
%   T is a types table to cast variables to a particular type, while keeping
%   the body of the algorithm unchanged.
%
%   This version of the algorithm has no scaling before the stages.
%
%   See also FI_RADIX2FFT_DEMO.

%   Reference:
%     Charles Van Loan, Computational Frameworks for the Fast Fourier
%     Transform, SIAM, Philadelphia, 1992, Algorithm 1.6.2, p. 45.
%
%   Copyright 2004-2015 The MathWorks, Inc.
%
%#codegen

    n = length(x);  t = log2(n);
    x = fi_bitreverse(x,n);

    % Generate index variables as integer constants so they are not computed in
    % the loop.
    LL = cast(2.^(1:t),'like',T.n);
    rr = cast(n./LL,'like',T.n);
    LL2 = cast(LL./2,'like',T.n);
    for q = 1:t
        L = LL(q); r = rr(q); L2 = LL2(q);
        for k = 0:(r-1)
            for j = 0:(L2-1)
                temp          = w(L2-1+j+1) * x(k*L+j+L2+1);
                x(k*L+j+L2+1) = bitsra(x(k*L+j+1) - temp, 1);
                x(k*L+j+1)    = bitsra(x(k*L+j+1) + temp, 1);
            end
        end
    end
end

固定小数点データでスケーリングされたアルゴリズムを実行します。

x = cast(x0,'like',T3.x);
w = cast(w0,'like',T3.w);

y = fi_m_radix2fft_withscaling_typed(x,w,T3);

fi_fft_demo_plot(real(x), y, y0/n, Fs, 'Fixed-point data', ...
    {'Fixed-point FFT with scaling','Built-in FFT'});

Figure contains 3 axes objects. Axes object 1 contains an object of type line. This object represents Fixed-point data. Axes object 2 contains 2 objects of type line. These objects represent Fixed-point FFT Algorithm 1.6.2, Built-in FFT. Axes object 3 contains an object of type line. This object represents abs(error).

Figure contains 3 axes objects. Axes object 1 contains an object of type line. This object represents Fixed-point data. Axes object 2 contains 2 objects of type line. These objects represent Fixed-point FFT with scaling, Built-in FFT. Axes object 3 contains an object of type line. This object represents abs(error).

スケーリングされた固定小数点の FFT アルゴリズムでは、組み込み FFT と 16 ビットの固定小数点データに必要な許容誤差が一致することがわかります。

クリーン アップ

以下のコードを実行してグローバル状態を復元します。

fipref(FIPREF_STATE);

参考文献

Charles Van Loan, Computational Frameworks for the Fast Fourier Transform, SIAM, 1992.

Cleve Moler, Numerical Computing with MATLAB, SIAM, 2004, Chapter 8 Fourier Analysis.

Alan V. Oppenheim and Ronald W. Schafer, Discrete-Time Signal Processing, Prentice Hall, 1989.

Peter D. Welch, "A Fixed-Point Fast Fourier Transform Error Analysis," IEEE® Transactions on Audio and Electroacoustics, Vol. AU-17, No. 2, June 1969, pp. 151-157.