ドキュメンテーション

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

fftfilt

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

構文

y = fftfilt(b,x)
y = fftfilt(b,x,n)
y = fftfilt(d,x)
y = fftfilt(d,x,n)
y = fftfilt(gpuArrayb,gpuArrayX,n)

説明

fftfilt では、FIR フィルターに対してのみ機能する周波数領域のフィルター処理手法である、効率的な FFT ベースの "オーバーラップ加算" 法を使用して、データをフィルター処理します。

y = fftfilt(b,x) では、ベクトル x 内のデータが係数ベクトル b の記述するフィルターでフィルター処理されます。これはデータ ベクトル y を返します。fftfilt で実行される演算は、差分方程式により "時間領域" 内で、以下のように表されます。

y(n)=b(1)x(n)+b(2)x(n1)++b(nb+1)x(nnb)

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

Y(z)=(b(1)+b(2)z1++b(nb+1)znb)X(z)

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

x が行列の場合、fftfilt では各列がフィルター処理されます。b が行列の場合、fftfilt では、b の各列にあるフィルターが信号ベクトル x に適用されます。bx が両方とも同じ列数をもつ行列の場合は、b の i 列目を使用して、x の i 列目がフィルター処理されます。

y = fftfilt(b,x,n)n を使用して FFT の長さを求めます。詳細は、アルゴリズムを参照してください。

y = fftfilt(d,x)digitalFilter オブジェクト d を使用してベクトル x のデータをフィルター処理します。d を周波数応答仕様に基づいて生成するには、関数 designfilt を使用します。

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

y = fftfilt(gpuArrayb,gpuArrayX,n) は、gpuArraygpuArrayb 内の FIR フィルター係数で gpuArray オブジェクト、gpuArrayX 内のデータをフィルター処理します。gpuArray オブジェクトの詳細については、GPU での MATLAB 関数の実行 (Parallel Computing Toolbox)を参照してください。fftfiltgpuArray オブジェクトと使用するには、Parallel Computing Toolbox™ ソフトウェアが必要です。どの GPU がサポートされているかについては、リリース別の GPU サポート (Parallel Computing Toolbox)を参照してください。フィルター処理されたデータ y は、gpuArray オブジェクトです。GPU でのオーバーラップ加算フィルター処理の例は、GPU でのオーバーラップ加算フィルター処理を参照してください。

fftfilt は、実数入力および複素数入力の両方に使用できます。

filter 関数との比較

入力信号が比較的大きい場合は、関数 filter よりも、x の各サンプルに対して N 回 (N はフィルター長) の乗算を行う fftfilt を使用する方が効率的です。fftfilt では、以下のコストで 2 つの FFT 演算 (信号ブロック長 L の FFT と、FFT 同士の積の逆 FFT) が実行されます。

½Llog2L

ここで、L はブロック長です。次に、L ポイントごとの乗算を実行します。総コストは次のようになります。

L + Llog2L = L(1 + log2L)

したがって、2 つの関数のコスト比は、以下のようになります。

L(1+log2L) / (NL) = (1 + log2L)/N

これは log2L / N に近似できます。

つまり、fftfilt は、log2L が N 未満である場合に効果的になります。

すべて折りたたむ

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

rng default

N = 100;

shrt = 20;
long = 2000;

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

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

    bshrt = rand(1,shrt);

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

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

    blong = rand(1,long);

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

end

平均時間を比較します。

fprintf('%4d-tap filter averages: fftfilt: %f s; filter: %f s\n',shrt,tfs,tls)
  20-tap filter averages: fftfilt: 0.412923 s; filter: 0.021208 s
fprintf('%4d-tap filter averages: fftfilt: %f s; filter: %f s\n',long,tfl,tll)
2000-tap filter averages: fftfilt: 0.104964 s; filter: 0.221990 s

この例では、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.Coefficients;

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

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

アルゴリズム

fftfiltfft を使用して、"オーバーラップ加算法" (参考文献 [1]) を実装します。この手法は、連続的な周波数領域のフィルター処理された入力シーケンスのブロックを組み合わせるものです。fftfilt は、入力シーケンス x を長さ L のデータ ブロックに分割します。このとき、L はフィルター長 N よりも長くなければなりません。

さらに、以下の操作により、フィルター b で各ブロックを畳み込みます。

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

ここで、nfft は FFT の長さです。fftfilt は、連続的な出力セクションを n-1 点ずつオーバーラップさせ、その和を計算します。ここで、n はフィルターの長さです。

fftfilt では、FFT の長さ n を指定するかどうか、またフィルターや信号の長さによって、さまざまな方法で主要パラメーターの Lnfft が選択されます。FFT の長さを決定する n の値を指定しない場合は、fftfilt によってはこれらの主要パラメーターが自動的に選択されます。

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

  • length(b)length(x) に等しいまたは大きい場合、fftfilt では、以下の長さの単一 FFT が使用されます。

    2^nextpow2(length(b) + length(x) - 1)
    

    これは基本的に、次の計算を行うものです。

    y = ifft(fft(B,nfft).*fft(X,nfft))
    

n に値を指定すると、fftfilt では FFT の長さ nfft2^nextpow2(n)、データ ブロック長に nfft - length(b) + 1 が選択されます。nlength(b) 未満の場合、fftfilt では nlength(b) が設定されます。

参照

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

拡張機能

R2006a より前に導入