手で微調整する際にアーティファクトを防ぐために FFT スペクトルを滑らかにする方法の経験則はありますか?

Question

FFT 振幅スペクトルがあり、それから周期的なノイズ源 (例: ノイズ源) を選択的に通過させるフィルターを作成したいと考えています。正弦波スプリアス）と、ランダムな背景ノイズに関連付けられた周波数ビンをゼロにします。このフィルタが時間領域に IFFT されると、周波数領域の急激な遷移によってリンギング アーティファクトが生成されることは理解しています...そこで、そのようなリンギングを避けるために、そのようなフィルターでの遷移を滑らかにする方法の経験則はあるのかどうか疑問に思っています。

たとえば、FFT に 1M の周波数ビンがあり、バックグラウンド ノイズ フロアから突き出ている 5 つのスプリアスがある場合、5 つのスプリアスそれぞれに関連付けられたピーク ビンを除くすべてのビンをゼロにしたいとします。問題は、時間領域でのアーティファクトを防ぐために隣接するスプリアス ビンをどのように処理するかです。たとえば、スプリアス ビンの両側のビンを 50% の振幅に設定する必要がありますか?拍車ビンの両側に 2 つのビンを使用する必要がありますか (50% で最も近いビン、25% で次に近いビンなど)。ご意見をお待ちしております。ありがとう！

Answer 1

私は次の方法が好きです。

• 理想的な振幅スペクトルを作成します (DC に関して対称になるように注意してください)。
• 時間領域への逆変換
• ブロックサイズの半分だけブロックを回転します
• Hann ウィンドウを適用する

あなたが提案しているほどシャープなものでそれを試したことはありませんが、かなり滑らかな周波数領域の結果が作成されることがわかりました。Kaiser-Bessel ウィンドウを使用すると、よりシャープなフィルターを作成できる可能性がありますが、パラメーターを適切に選択する必要があります。シャープにすることで、おそらくサイドローブを 6 dB 程度削減できると思います。

以下に Matlab/Octave コードのサンプルをいくつか示します。結果をテストするために使用しました freqz(h, 1, length(h)*10);.

function [ht, htrot, htwin] = ArbBandPass(N, freqs)
%# N = desired filter length
%# freqs = array of frequencies, normalized by pi, to turn into passbands
%# returns raw, rotated, and rotated+windowed coeffs in time domain

if any(freqs >= 1) || any(freqs <= 0)
    error('0 < passband frequency < 1.0 required to fit within (DC,pi)')
end

hf = zeros(N,1); %# magnitude spectrum from DC to 2*pi is intialized to 0
%# In Matlabs FFT, idx 1 -> DC, idx 2 -> bin 1, idx N/2 -> Fs/2 - 1, idx N/2 + 1 -> Fs/2, idx N -> bin -1
idxs = round(freqs * N/2)+1; %# indeces of passband freqs between DC and pi
hf(idxs) = 1; %# set desired positive frequencies to 1
hf(N - (idxs-2)) = 1; %# make sure 2-sided spectrum is symmetric, guarantees real filter coeffs in time domain
ht = ifft(hf); %# this will have a small imaginary part due to numerical error
if any(abs(imag(ht)) > 2*eps(max(abs(real(ht)))))
    warning('Imaginary part of time domain signal surprisingly large - is the spectrum symmetric?')
end
ht = real(ht); %# discard tiny imag part from numerical error
htrot = [ht((N/2 + 1):end) ; ht(1:(N/2))]; %# circularly rotate time domain block by N/2 points
win = hann(N, 'periodic'); %# might want to use a window with a flatter mainlobe
htwin = htrot .* win;
htwin = htwin .* (N/sum(win)); %# normalize peak amplitude by compensating for width of window lineshape