20

私は 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

出力は float のペアで構成され、それぞれが 1 つのビンの実部と虚部を表します。したがって、ビン 0 (配列インデックス 0、1) は、DC 周波数の実部と虚部を表します。ご覧のとおり、ビン 1 と 3 は両方とも同じ値 (Im 部分の符号を除く) を持っているため、ビン 3 が最初の負の周波数であり、最後にインデックス (4, 5) が最後の正の値であると推測します周波数ビン。

次に、周波数ビン 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;
}

明らかに私は何か間違ったことをしていますが、何がわかりません。

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

周波数領域でフィルタリングする正しい方法は何ですか? 私は何が欠けていますか？

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

前もって感謝します。

4

4 に答える 4

36

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

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

これが例です。周期が 1.5 サイクルの正弦波入力信号xと、単純なローパス フィルターを考えます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

ブロックの最初のグリッチは、私たちが期待するものではありません。しかし、あなたが考えるならばfft(x)、それは理にかなっています. 離散フーリエ変換は、信号が変換ブロック内で周期的であると仮定します。DFT が知る限り、次の 1 周期の変換を求めました。 DFT への非周期的な入力

これは、DFT でフィルタリングする際の最初の重要な考慮事項につながります。実際には、線形畳み込みではなく循環畳み込みを実装しています。したがって、最初のグラフの「グリッチ」は、数学を考えると、実際にはグリッチではありません。では、問題は次のようになります。周期性を回避する方法はありますか? 答えはイエスです。オーバーラップ保存処理を使用してください。基本的に、上記のように 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ます： 正しい

この図は理にかなっています - フィルタからのスタートアップ トランジェントが予想され、結果は定常状態の正弦波応答に落ち着きます。xnowは任意の長さにできることに注意してください。制限はNproc > 2*min(length(x),length(h))です。

2 つ目の問題は、特定のフィルターです。あなたのループでは、スペクトルが本質的にあるフィルターを作成します。そうすると、次のようH = [0 (1:511)/512 1 (511:-1:1)/512]'; にhraw = real(ifft(H)); plot(hraw)なります。 ひらめく

見にくいですが、グラフの左端にゼロ以外の点が多数あり、右端にさらに多くの点があります。Octave の組み込み関数を使用して、( を実行して) freqz周波数応答を確認します。freqz(hraw)freqz(hraw)

振幅応答には、ハイパス エンベロープからゼロまでの多くのリップルがあります。ここでも、DFT 固有の周期性が機能しています。DFTに関する限り、hraw何度も繰り返します。hrawしかし、の 1 周期を取るとfreqz、そのスペクトルは周期的なものとはまったく異なります。

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

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

速度のために、通常は padded の DFT を事前に計算することに注意してくださいh。元のデータと簡単に比較できるように、ループ内にそのまま残しました。

于 2010-06-01T11:12:49.303 に答える
11

あなたの主な問題は、短い時間間隔で頻度が明確に定義されていないことです。これは特に低周波に当てはまります。そのため、そこで最も問題に気付くのです。

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

たとえば、いくつかの妥当な数値を取ります: 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 までの短いセグメントに可能なすべての周波数を含めて初めて、元の正弦波の適切な表現が得られます。

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

于 2010-06-06T04:42:12.880 に答える
2

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 を実行します (入力信号よりも大きな 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 を実行して時間ドメインに戻すことにより、各入力オーディオ バッファーに変換します。これは、サンプルごとに行っているすべての分岐よりもはるかに効率的です。

于 2010-06-01T07:58:42.363 に答える
1

まず、正規化について: これは既知の (非) 問題です。DFT/IDFTでは、それらを対称的で真に可逆にするために、それぞれ (標準のコサイン/サイン ファクターとは別に) 1/sqrt(N)のファクター (逆を指示する) が必要です。もう 1 つの可能性は、それらの 1 つ (直接または逆) をNで割ることです。これは、利便性と好みの問題です。多くの場合、FFT ルーチンはこの正規化を実行しません。ユーザーはそれを認識して、好みに応じて正規化することが期待されます。見る

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

于 2010-05-28T13:47:57.443 に答える