fftfilt

オーバーラップ加算法を使用した FFT ベースの FIR フィルター処理

構文

y = fftfilt(b,x)

y = fftfilt(b,x,nfft)

y = fftfilt(d,x)

y = fftfilt(d,x,nfft)

説明

y = fftfilt(b,x) は、ベクトル x で指定されたデータをフィルター処理します。この関数は、係数ベクトル b によって記述されるフィルターを使用します。

y = fftfilt(b,x,nfft) は nfft を使用して FFT 長を求めます。

y = fftfilt(d,x) は digitalFilter オブジェクト d を使用してベクトル x のデータをフィルター処理します。

y = fftfilt(d,x,nfft) は nfft を使用して FFT 長を求めます。

例

短いフィルターおよび長いフィルター用の fftfilt および filter

小さなオペランドにはfilterの方が、大きなオペランドにはfftfiltの方が効率的であることを確認します。20 タップの短いフィルターと 2,000 タップの長いフィルターの 2 つの乱数フィルターを使用して、${10}^6$ 個の乱数をフィルター処理します。tic と toc を使用して実行時間を測定します。実験を 100 回繰り返して統計を向上させます。

rng("default")

N = 100;

s = 20;
l = 2000;

tfs = 0;
tls = 0;
tfl = 0;
tll = 0;

for kj = 1:N
    
    x = rand(1,1e6);

    bshrt = rand(1,s);

    tic
    sfs = fftfilt(bshrt,x);
    tfs = tfs+toc/N;

    tic
    sls = filter(bshrt,1,x);
    tls = tls+toc/N;

    blong = rand(1,l);

    tic
    sfl = fftfilt(blong,x);
    tfl = tfl+toc/N;
    
    tic
    sll = filter(blong,1,x);
    tll = tll+toc/N;

end

平均時間を比較して表示します。

table(table(1000*[tfs;tls],1000*[tfl;tll], ...
    RowNames=["fftfilt" "filter"],VariableNames=[s;l]+"-tap"), ...
    VariableNames="Filter averages (milliseconds)")

ans=2×1 table
            Filter averages (milliseconds)        
    ______________________________________________

               20-tap    2000-tap
               ______    ________
                                                  
    fftfilt    142.65     60.027 
    filter     10.064     149.01 

GPU でのオーバーラップ加算フィルター処理

この例では、Parallel Computing Toolbox™ ソフトウェアが必要です。サポートされている GPU の一覧については、GPU 計算の要件 (Parallel Computing Toolbox)を参照してください。

加法性ホワイト ガウス ノイズの正弦波の合計から構成される信号を作成します。正弦波の周波数は、2.5、5、10 および 15 kHz です。サンプリング周波数は 50 kHz です。

Fs = 50e3;
t = 0:1/Fs:10-(1/Fs);
x = cos(2*pi*2500*t) + ...
    0.5*sin(2*pi*5000*t) + ...
    0.25*cos(2*pi*10000*t)+ ...
    0.125*sin(2*pi*15000*t) + randn(size(t));

designfiltを使用してローパス FIR 等リップル フィルターを設計します。

d = designfilt("lowpassfir",SampleRate=Fs, ...
    PassbandFrequency=5500,StopbandFrequency=6000, ...
    PassbandRipple=0.5,StopbandAttenuation=50);
B = d.Numerator;

オーバーラップ加算法によって、GPUでデータをフィルター処理します。gpuArray を使用してデータを GPU に入力します。gather を使用して出力を MATLAB® ワークスペースに返し、フィルターされたデータのパワー スペクトル密度推定をプロットします。

y = fftfilt(gpuArray(B),gpuArray(x));
periodogram(gather(y),rectwin(length(y)),length(y),50e3)

入力引数

b — フィルター係数

ベクトル | 行列

フィルター係数。ベクトルとして指定します。b が行列の場合、fftfilt では、b の各列にあるフィルターが信号ベクトル x に適用されます。

x — 入力データ

ベクトル | 行列

入力データ。ベクトルとして指定します。x が行列の場合、fftfilt では各列がフィルター処理されます。b と x が同じ列数を持つ行列の場合、b の i 列目を使用して x の i 列目がフィルター処理されます。fftfilt は、実数入力および複素数入力の両方に使用できます。

nfft — FFT 長

正の整数

FFT 長。フィルターの長さ以上の 2 のべき乗の正の整数として指定します。

• 既定では、fftfilt では効率的な実行時間になるように FFT の長さとデータ ブロック長が選択されます。

• nfft を 2 のべき乗の正の整数として指定しない場合、または nfft がフィルターの長さより小さい場合、fftfilt は 2^nextpow2(max(nfft,N)) の FFT 長を選択します。ここで、N はフィルター長です。

  • b がベクトルの場合、N = numel(b)。

  • b が行列の場合、N = size(b,1)。

d — デジタル フィルター

digitalFilter オブジェクト

デジタル フィルター。digitalFilter オブジェクトで指定します。d を周波数応答仕様に基づいて生成するには、関数 designfilt を使用します。

出力引数

y — 出力データ

ベクトル | 行列

出力データ。ベクトルまたは行列として返されます。

詳細

filter 関数との比較

入力信号が比較的大きい場合は、fftfilt が filter よりも高速です。

filter は x の各サンプルに対して N 回 (N はフィルター長) の乗算を行います。fftfilt では、$$\frac{1}{2}L{\log }_{2}L$$ のコストで 2 つの FFT 演算 (信号ブロック長 L の FFT と、FFT どうしの積の逆 FT) が実行されます。ここで、L はブロック長です。次に、L ポイントごとの乗算を実行します。総コストは $$L+L{\log }_{2}L=L(1+{\log }_{2}L)$$ 回の乗算です。そのため、コスト比は $$L(1+{\log }_{2}L)/(NL)=(1+{\log }_{2}L)/N$$ になります。これは log₂L / N に近似できます。

つまり、fftfilt は、log₂L が N 未満である場合により高速になります。

アルゴリズム

fftfilt では、FIR フィルターに対してのみ機能する周波数領域のフィルター処理手法である、効率的な FFT ベースの "オーバーラップ加算" 法[1]を使用して、連続的な周波数領域のフィルター処理された入力シーケンスのブロックを組み合わせることで、データをフィルター処理します。

N_x 個の要素をもつ入力信号ベクトル x と N 個の要素をもつフィルター ベクトル b があると仮定します。ここで、b = [b₁ b₂ ⋯ b_N] です。fftfilt で実行される演算は、差分方程式により時間領域内で、以下のように表されます。

$$y[n]=b_{\text{}1}\text{}x[n]+b_{\text{}2}\text{}x[n-1]+\cdots +b_{\text{}N}\text{}x[n-(N-1)]$$

これと等価な表現は、Z 変換、すなわち "周波数領域" 表現で、以下のように表されます。

$$Y(z)=\left[b_{\text{}\text{}1}+b_{\text{}2}\text{}{z}^{-1}+\cdots +b_{\text{}N}\text{}{z}^{-(N-1)}\right]X(z)$$

fftfilt 関数は、fft を使用してオーバーラップ加算法を実装します。fftfilt は、入力シーケンス x を長さが L の k 個のデータ ブロックに分割します。ここで、L はフィルターの長さ N よりも大きくなければなりません。ここで、k = ⌈N_x/L⌉ および ⌈⌉ の記号は天井関数を表します。

関数は、次を使用して各ブロックをフィルター b で畳み込みます。

y = ifft(fft(x(i:i+L-1),nfft).*fft(b,nfft));

ここで、i = 1, L+1, 2L+1, ⋯ と nfft は FFT 長です。fftfilt は、連続する出力セクションを N–1 点ずつオーバーラップさせ、その和を計算します。

fftfilt は、フィルターと信号の FFT 長 nfft が指定されているかどうかに応じ、異なる方法で主要パラメーターの L と nfft を選択します。

• nfft (これにより FFT 長が決まる) の値が指定されていない場合、fftfilt は次の主要パラメーターを自動的に選択します。

  • length(x) が length(b) より大きい場合、fftfilt では、ブロック数と FFT ごとのフロップ数を乗算した結果が最小になる値が選択されます。

  • length(b) が length(x) 以上の場合、fftfilt は、長さ 2^nextpow2(length(b) + length(x) - 1) の単一 FFT を使用します。

    これらの仮定により、y が次のように得られます。

    y = ifft(fft(b,nfft).*fft(x,nfft))
    

• nfft に値が指定されている場合、fftfilt は、FFT 長として 2^nextpow2(nfft) を、データ ブロック長として nfft-length(b)+1 を選択します。nfft が length(b) 未満の場合、fftfilt は、2^nextpow(length(b)) の FFT 長を選択します。

参照

[1] Oppenheim, Alan V., Ronald W. Schafer, and John R. Buck. Discrete-Time Signal Processing. 2nd Ed. Upper Saddle River, NJ: Prentice Hall, 1999.

拡張機能

C/C++ コード生成

MATLAB® Coder™ を使用して C および C++ コードを生成します。

使用上の注意および制限:

digitalFilter オブジェクトは、コード生成でサポートされません。

GPU コード生成

GPU Coder™ を使用して NVIDIA® GPU のための CUDA® コードを生成します。

使用上の注意および制限:

digitalFilter オブジェクトは、コード生成でサポートされません。

GPU 配列

Parallel Computing Toolbox™ を使用してグラフィックス処理装置 (GPU) 上で実行することにより、コードを高速化します。

この関数は、GPU 配列を完全にサポートします。詳細については、GPU での MATLAB 関数の実行 (Parallel Computing Toolbox)を参照してください。

バージョン履歴

R2006a より前に導入

参考

conv | designfilt | digitalFilter | filter | filtfilt