DSP- FFT経由の周波数ドメインでのフィルタリング

https://stackoverflow.com/questions/2929401

05-10-2019
|

質問

私はFFTのExocortexの実装で少し遊んでいますが、いくつかの問題があります。

IFFTを呼び出す前に周波数ビンの振幅を変更するたびに、結果の信号には、特に信号に低い周波数が存在する場合（ドラムやベースなど）。ただし、これはすべてのビンを同じ要因で減衰させた場合には発生しません。

4サンプルFFTの出力バッファーの例を挙げましょう。

// Bin 0 (DC)
FFTOut[0] = 0.0000610351563
FFTOut[1] = 0.0

// Bin 1
FFTOut[2] = 0.000331878662
FFTOut[3] = 0.000629425049

// Bin 2
FFTOut[4] = -0.0000381469727
FFTOut[5] =  0.0

// Bin 3, this is the first and only negative frequency bin.
FFTOut[6] =  0.000331878662
FFTOut[7] = -0.000629425049

出力は、フロートのペアで構成されており、それぞれが単一のビンの実際と想像力の部分を表しています。したがって、bin 0（配列インデックス0、1）は、DC周波数の実際の部分と想像上の部分を表します。ご覧のとおり、ビン1と3の両方が同じ値を持っています（IM部分の兆候を除く）。周波数ビン。

次に、周波数bin 1を減衰させる1これが私がしていることです：

// Attenuate the 'positive' bin
FFTOut[2] *= 0.5;
FFTOut[3] *= 0.5;

// Attenuate its corresponding negative bin.
FFTOut[6] *= 0.5;
FFTOut[7] *= 0.5;

実際のテストでは、1024レングスのFFTを使用していますが、常にすべてのサンプルを提供しているため、0パディングは必要ありません。

// Attenuate
var halfSize = fftWindowLength / 2;
float leftFreq = 0f;
float rightFreq = 22050f; 
for( var c = 1; c < halfSize; c++ )
{
    var freq = c * (44100d / halfSize);

    // Calc. positive and negative frequency indexes.
    var k = c * 2;
    var nk = (fftWindowLength - c) * 2;

    // This kind of attenuation corresponds to a high-pass filter.
    // The attenuation at the transition band is linearly applied, could
    // this be the cause of the distortion of low frequencies?
    var attn = (freq < leftFreq) ? 
                    0 : 
                    (freq < rightFreq) ? 
                        ((freq - leftFreq) / (rightFreq - leftFreq)) :
                        1;

    // Attenuate positive and negative bins.
    mFFTOut[ k ] *= (float)attn;
    mFFTOut[ k + 1 ] *= (float)attn;
    mFFTOut[ nk ] *= (float)attn;
    mFFTOut[ nk + 1 ] *= (float)attn;
}

明らかに私は何か間違ったことをしていますが、何を理解することはできません。

非常に基本的なダイナミックイコライザーを実装しようとしているため、FFT出力をFFT係数を生成する手段として使用したくありません。

周波数ドメインでフィルタリングする正しい方法は何ですか？何が足りないの？

また、負の周波数も減衰する必要がありますか？ NEGがあるFFT実装を見てきました。周波数値は合成前にゼロになります。

前もって感謝します。

解決

2つの問題があります。FFTと特定のフィルターの使用方法です。

フィルタリングは、従来、時間領域の畳み込みとして実装されています。入力信号とフィルター信号のスペクトルを掛けることは同等です。ただし、離散フーリエ変換（DFT）（高速フーリエ変換アルゴリズムで実装された速度で実装）を使用すると、実際に真のスペクトルのサンプリングバージョンを計算します。これには多くの意味がありますが、フィルタリングに最も関連するものは、時間ドメイン信号が周期的であるという意味です。

これが例です。正弦波入力信号を検討してください x 期間に1.5サイクル、および単純なローパスフィルターがあります h. 。 Matlab/Octave構文で：

N = 1024;
n = (1:N)'-1; %'# define the time index
x = sin(2*pi*1.5*n/N); %# input with 1.5 cycles per 1024 points
h = hanning(129) .* sinc(0.25*(-64:1:64)'); %'# windowed sinc LPF, Fc = pi/4
h = [h./sum(h)]; %# normalize DC gain

y = ifft(fft(x) .* fft(h,N)); %# inverse FT of product of sampled spectra
y = real(y); %# due to numerical error, y has a tiny imaginary part
%# Depending on your FT/IFT implementation, might have to scale by N or 1/N here
plot(y);

そして、これがグラフです： IFFT of product

ブロックの先頭にあるグリッチは、私たちが期待するものではありません。しかし、あなたが考慮するなら fft(x), 、それは理にかなっている。離散フーリエ変換は、信号が変換ブロック内で周期的であると仮定します。 DFTが知っている限りでは、これの1つの期間の変換を求めました。 Aperiodic input to DFT

これは、DFTSでフィルタリングするときに最初の重要な考慮事項につながります：あなたは実際に実装しています円形の畳み込み, 、線形畳み込みではありません。したがって、最初のグラフの「グリッチ」は、数学を考慮するとき、本当にグリッチではありません。それで、質問は次のとおりです。周期性を回避する方法はありますか？答えはイエスです：使用オーバーラップ節約処理. 。基本的に、上記のようにNロング製品を計算しますが、N/2ポイントのみを保持します。

Nproc = 512;
xproc = zeros(2*Nproc,1); %# initialize temp buffer
idx = 1:Nproc; %# initialize half-buffer index
ycorrect = zeros(2*Nproc,1); %# initialize destination
for ctr = 1:(length(x)/Nproc) %# iterate over x 512 points at a time
    xproc(1:Nproc) = xproc((Nproc+1):end); %# shift 2nd half of last iteration to 1st half of this iteration
    xproc((Nproc+1):end) = x(idx); %# fill 2nd half of this iteration with new data
    yproc = ifft(fft(xproc) .* fft(h,2*Nproc)); %# calculate new buffer
    ycorrect(idx) = real(yproc((Nproc+1):end)); %# keep 2nd half of new buffer
    idx = idx + Nproc; %# step half-buffer index
end

そして、これがのグラフです ycorrect:

この写真は理にかなっています - フィルターから一時的なスタートアップが期待され、結果は定常状態の正弦波応答に落ち着きます。今それに注意してください x 任意に長くすることができます。制限はです Nproc > 2*min(length(x),length(h)).

次に、2番目の問題について：特定のフィルター。ループでは、スペクトルが本質的にあるフィルターを作成します H = [0 (1:511)/512 1 (511:-1:1)/512]'; もし、するなら hraw = real(ifft(H)); plot(hraw), 、あなたは得ます： hraw

見るのは難しいですが、グラフの左端にはゼロ以外のポイントがたくさんあり、その後、右端にさらに束があります。 Octaveの組み込みを使用します freqz 私たちが見る周波数応答を見るための関数（することによって freqz(hraw)):

マグニチュード応答には、ハイパスエンベロープからゼロまで多くの波紋があります。繰り返しますが、DFTに固有の周期性は機能しています。 DFTに関する限り、 hraw 何度も何度も繰り返します。しかし、あなたが1つの期間をとる場合 hraw, 、なので freqz そのスペクトルは、周期バージョンとはまったく異なります。

それでは、新しい信号を定義しましょう。 hrot = [hraw(513:end) ; hraw(1:512)]; 生のDFT出力を単純に回転させて、ブロック内で連続するようにします。次に、使用を使用して周波数応答を見てみましょう freqz(hrot):

ずっといい。望ましい封筒は、すべての波紋がなく、そこにあります。もちろん、実装はそれほど単純ではありません。完全に複雑な乗算を行う必要があります fft(hrot) 各複雑なビンを単にスケーリングするのではなく、少なくとも正しい答えが得られます。

速度では、通常、パッド入りのDFTを事前に計算することに注意してください h, 、私はそれをループに残して、オリジナルとより簡単に比較しました。

他のヒント

あなたの主な問題はそれです 周波数は短い時間間隔ではあまり定義されていません. 。これは、低周波数に特に当てはまります。そのため、問題に最も問題に気付くのはそのためです。

したがって、サウンドトレインから本当に短いセグメントを取り出して、これらをフィルタリングすると、フィルタリングされたセグメントが連続波形を生成する方法でフィルターをかけず、セグメント間のジャンプが聞こえます。これがここでクリックを生成するものです。

たとえば、いくつかの合理的な数字を取る：27.5 Hz（ピアノ上のA0）の波形から始め、44100 Hzでデジタル化され、このように見えます（赤い部分は1024サンプルの長さです）：

それで、最初に40Hzのローパスから始めます。したがって、元の周波数は40Hz未満であるため、40Hzのカットオフを備えたローパスフィルターは実際には効果がありません。入力とほぼ一致する出力を取得します。右？ 間違っている、間違っている、間違っている - そして、これは基本的にあなたの問題の中核です。問題は、短いセクションの場合です アイディア 27.5 Hzは明確に定義されておらず、DFTでうまく表現することはできません。

27.5 Hzは、以下の図のDFTを見ることで、短いセグメントでは特に意味がありません。長いセグメントのDFT（黒い点）は27.5 Hzのピークを示していますが、短いもの（赤い点）はそうではないことに注意してください。

明らかに、40Hz以下のフィルタリングはDCオフセットをキャプチャするだけで、40Hzのローパスフィルターの結果を以下の緑色に示します。

青い曲線（200 Hzのカットオフで撮影）は、はるかに良く一致し始めています。しかし、それをうまく一致させているのは低周波数ではなく、高周波数が含まれていることに注意してください。最終的な22kHzまでのすべての周波数を最大22kHzに含めるまで、最終的に元の正弦波の適切な表現を得ることができません。

このすべての理由は、27.5 Hzの正弦波の小さなセグメントが いいえ 27.5 Hzの正弦波、DFTは27.5 Hzとはあまり関係がありません。

DC周波数サンプルの値をゼロに減衰させていますか？あなたはあなたの例でそれをまったく減衰させていないようです。ハイパスフィルターを実装しているため、DC値をゼロに設定する必要があります。

これにより、低周波数の歪みが説明されます。そのDC値が大きな遷移のために非ゼロである場合、低周波数での周波数応答には多くのリップルがあります。

ここにMatlab/Octaveの例があり、何が起こっているのかを示します。

N = 32;
os = 4;
Fs = 1000;
X = [ones(1,4) linspace(1,0,8) zeros(1,3) 1 zeros(1,4) linspace(0,1,8) ones(1,4)];
x = ifftshift(ifft(X));
Xos = fft(x, N*os);
f1 = linspace(-Fs/2, Fs/2-Fs/N, N);
f2 = linspace(-Fs/2, Fs/2-Fs/(N*os), N*os);

hold off;
plot(f2, abs(Xos), '-o');
hold on;
grid on;
plot(f1, abs(X), '-ro');
hold off;
xlabel('Frequency (Hz)');
ylabel('Magnitude');

私のコードでは、DC値がゼロ以外であるという例を作成しており、その後にゼロに急激に変更され、その後ランプアップが行われていることに注意してください。次に、IFFTを使用して時間領域に変換します。次に、その時間領域信号で、入力信号よりも大きいFFTサイズを渡すときにMATLABによって自動的に行われるゼロパッドFFT（これはMATLABによって自動的に行われます）を実行します。時間領域のゼロパッジングは、周波数領域に補間されます。これを使用して、フィルターがフィルターサンプル間でどのように応答するかを確認できます。

覚えておくべき最も重要なことの1つは、DFTの出力を減衰させることにより、指定された周波数でフィルター応答値を設定しているにもかかわらず、サンプルポイント間で発生する周波数に対して何も保証されないことです。これは、変化が急激に急激になるほど、サンプル間のオーバーシュートと振動が発生することを意味します。

ここで、このフィルタリングをどのように行うべきかについての質問に答えてください。いくつかの方法がありますが、実装して理解するのが最も簡単な方法の1つは、ウィンドウ設計方法です。現在の設計の問題は、遷移幅が巨大であることです。ほとんどの場合、可能な限り波紋が少なく、できるだけ迅速に移行する必要があります。

次のコードでは、理想的なフィルターを作成し、応答を表示します。

N = 32;
os = 4;
Fs = 1000;
X = [ones(1,8) zeros(1,16) ones(1,8)];
x = ifftshift(ifft(X));
Xos = fft(x, N*os);
f1 = linspace(-Fs/2, Fs/2-Fs/N, N);
f2 = linspace(-Fs/2, Fs/2-Fs/(N*os), N*os);

hold off;
plot(f2, abs(Xos), '-o');
hold on;
grid on;
plot(f1, abs(X), '-ro');
hold off;
xlabel('Frequency (Hz)');
ylabel('Magnitude');

突然の変化によって引き起こされる多くの振動があることに注意してください。

FFTまたは離散フーリエ変換は、フーリエ変換のサンプリングバージョンです。フーリエ変換は、DFTが有限数のサンプルに適用される一方で、連続範囲 - インフィニティを介して信号に適用されます。これにより、DFTを使用するときの時間領域では、四角いウィンドウ（切り捨て）が得られます。残念ながら、正方形の波のdftはsincタイプ関数（sin（x）/x）です。

フィルターに急激な遷移があること（1つのサンプルで0から1までのクイックジャンプ）の問題は、これが正方形のウィンドウによって切り捨てられている時間領域で非常に長い応答があることです。したがって、この問題を最小限に抑えるために、タイムドメイン信号をより段階的なウィンドウで掛けることができます。ラインを追加してハニングウィンドウを掛けた場合：

x = x .* hanning(1,N).';

IFFTを取得した後、この応答を取得します。

したがって、ウィンドウ設計方法を実装することをお勧めします。これはかなり簡単だからです（より良い方法がありますが、より複雑になります）。あなたはイコライザーを実装しているので、私はあなたがその場で減衰を変更できるようにしたいと思いますので、パラメーターに変更があるときはいつでも周波数領域にフィルターを計算して保存することをお勧めします、そしてあなたはそれを適用することができます入力バッファーのFFTを取得し、周波数ドメインフィルターサンプルを乗算し、IFFTを実行して時間ドメインに戻ることにより、各入力オーディオバッファーに。これは、各サンプルに対してあなたがしている分岐のすべてよりもはるかに効率的です。

まず、正規化について：それは既知の（非）問題です。 DFT/IDFTには要因が必要です 1/sqrt（n） （標準的なコサイン/サイン因子を除く）それぞれ（逆に指示）して、それらを模倣して真に反転させるものにします。別の可能性は、それらの1つ（直接または逆）を分割することです n, 、これは便利さと味の問題です。多くの場合、FFTルーチンはこの正規化を実行しないため、ユーザーはそれを認識し、好むように正規化することが期待されます。見る

2番目：（たとえば）16ポイントdftで、あなたが呼ぶもの ビン0 ゼロ周波数（DC）に対応し、 ビン1 低いfreq ... ビン4 ミディアムフレック、 ビン8 最高の周波数に ビン9 ... 15 「負の周波数」へ。あなたの例では、 ビン1 実際、低周波と中周波数の両方です。この考慮事項とは別に、「イコライゼーション」には概念的に間違っていることはありません。私はあなたが何を意味するのかわかりません 「信号は低周波数で歪んでいます」. 。それをどのように観察しますか？

ライセンス： CC-BY-SA と帰属

所属していません StackOverflow