DSP-通過FFT在頻域進行濾波

Question

我一直在玩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 

輸出由對浮體，每個代表一個單個倉的實部和imaginay零件。因此，bin 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; 
} 

顯然我做錯了什麼，但不知道是什麼。

我不想使用FFT輸出作爲生成一組FIR係數的手段，因爲我試圖實現一個非常基本的動態均衡器。

在頻域中過濾的正確方法是什麼？我錯過了什麼？

另外，是否真的需要衰減負頻率呢？我見過一個FFT實現，其中neg。頻率值在合成之前歸零。

在此先感謝。

Answer 1

有兩個問題：您使用FFT的方式和特定的過濾器。

過濾傳統上是作爲時域卷積來實現的。你說得對，輸入信號和濾波器信號的頻譜是相等的。但是，當您使用離散傅立葉變換（DFT）（使用快速傅立葉變換算法實現速度）時，實際上可以計算真實頻譜的採樣版本。這有很多含義，但與濾波最相關的一個含義是時域信號是週期性的。

下面是一個例子。考慮在該週期中具有1.5個週期的正弦輸入信號x以及簡單的低通濾波器h。在Matlab /倍頻語法：

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); 

而這裏的圖：

在該塊開頭的故障是不是我們所期望的。但如果你考慮fft(x)，這是有道理的。離散傅立葉變換假定信號在變換塊內是週期性的。至於DFT知道，我們要求變換的這一時期：

與DFT的過濾時，這導致第一重要的考慮因素：你實際上是實現circular convolution，而不是線性卷積。所以當考慮數學時，第一個圖中的「故障」並不是真正的小故障。那麼問題就變成了：是否有辦法解決週期性問題？答案是肯定的：使用overlap-save processing。基本上，你計算如上的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))。

現在到第二個問題：特定的過濾器。在你的循環，你創建一個過濾器誰的譜實質上是H = [0 (1:511)/512 1 (511:-1:1)/512]';如果你這樣做hraw = real(ifft(H)); plot(hraw)，您可以：

很難看，但也有在圖表的最左邊一堆非零點，然後在最右邊的一堆。使用倍頻的內置freqz功能看，我們看到（這樣做freqz(hraw)）頻率響應：

幅度響應有很多來自高通絡的漣漪下降到零。同樣，DFT固有的週期性正在起作用。就DFT而言，hraw一遍又一遍地重複。但如果你的hraw一個時期，隨着freqz呢，它的頻譜是週期性版本的完全不同。

所以，讓我們定義一個新的信號：hrot = [hraw(513:end) ; hraw(1:512)];我們只需旋轉的原始DFT輸出，使之在塊內連續。現在讓我們看看使用freqz(hrot)的頻率響應：

好得多。期望的信封在那裏，沒有所有的漣漪。當然，執行不那麼簡單，現在，你必須fft(hrot)而不僅僅是縮放每個復斌做了充分的複數乘法，但至少你會得到正確的答案。

請注意，對於速度，您通常會預先計算填充的h的DFT，我將其單獨留在循環中以便於與原始文件進行比較。

Answer 2

首先，關於規範化：這是一個已知（非）問題。 DFT/IDFT需要在每一個因子中（除了標準餘弦/正弦因子外）（直接相反）使其相互正比並且真正可逆。另一種可能性是將它們中的一個（直接或倒數）除以N，這是方便和品味的問題。通常FFT例程不執行這種標準化，用戶需要意識到這一點並按照他的喜好進行標準化。 See

二：在（比方說）16點DFT，你叫什麼斌0將對應於零頻率（DC），紙槽1個低頻率... 紙槽4中等頻率， bin 8到最高頻率和分檔9 ... 15到「負頻率」。在你的例子中，bin 1實際上都是低頻和中頻。除了這個考慮之外，你的「均衡」在概念上沒有任何錯誤。我不明白你的意思是「信號在低頻時失真」。你怎麼看？

Answer 3

您是否將DC頻率採樣值減小到零？在你的例子中，你似乎並沒有削弱它。由於您正在實施高通濾波器，因此還需要將DC值設置爲零。

這將解釋低頻失真。如果DC值由於較大的轉換而不爲零，則在低頻時頻率響應會出現大量波動。

這裏是MATLAB /八度爲例來演示一下可能發生的：

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'); 

frequency response http://www.freeimagehosting.net/uploads/e10109e535.png

注意，在我的代碼，我創建了DC值的一個例子是非零，然後突然變爲零，然後斜升。然後我將IFFT轉換到時域。然後，我在該時域信號上執行零填充fft（當您通過比輸入信號更大的fft尺寸時，由MATLAB自動完成）。時域中的零填充會導致頻域中的插值。使用這個，我們可以看到濾波器如何在濾波器樣本之間做出響應。

一個要記住最重要的事情是，即使你被衰減DFT的輸出，在給定的頻率設置濾波器響應值，這保證沒有爲樣本點之間發生的頻率。這意味着你的變化越突然，樣品之間的過沖和振盪就會發生。

現在回答你的問題，如何完成這個過濾。有很多方法，但最容易實現和理解的是窗口設計方法。目前設計的問題是轉換寬度很大。大多數情況下，儘可能快地進行轉換，儘可能減少波動。

在接下來的代碼，我將創建一個理想的濾波器並顯示響應：

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'); 

frequency response http://www.freeimagehosting.net/uploads/c86f5f1700.png

注意，有很多引起的突然變化振盪。

FFT或離散傅立葉變換是傅立葉變換的採樣版本。傅立葉變換應用於連續範圍的信號 - 無窮大到無窮大，同時將DFT應用於有限數量的採樣。由於我們只處理有限數量的採樣，因此這在使用DFT時會導致時域中的正方形窗口（截斷）。不幸的是，方波的DFT是sinc型函數（sin（x）/ x）。

與具有在過濾器的尖銳過渡（快速跳轉從0到1在一個樣品中）的問題是，這具有在時域中，它是由一個方形窗截斷很長的響應。因此，爲了最大限度地減少這個問題，我們可以將時域信號乘以更漸進的窗口。如果我們通過增加線乘以寧窗口：

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

服用IFFT後，我們得到這樣的迴應：

frequency response http://www.freeimagehosting.net/uploads/944da9dd93.png

所以我會建議您嘗試實現，因爲它的櫥窗設計方法相當簡單（有更好的方法，但它們變得更復雜）。由於您正在實現均衡器，我假設您希望能夠即時更改衰減，所以我建議每當參數發生變化時都要在頻域中計算和存儲濾波器，然後您可以將其應用於頻域通過將輸入緩衝區的fft乘以您的頻域濾波器樣本，然後執行ifft以返回到時域，從而對每個輸入音頻緩衝區進行調整。這將比您爲每個樣本進行的所有分支更有效率。

Answer 4

您的主要問題是：頻率不超過短的時間間隔定義。這對於低頻尤其如此，這就是爲什麼你最關注的問題。因此，當您從聲音系列中取出非常短的片段，然後過濾這些片段時，過濾後的片段不會以產生連續波形的方式進行過濾，並且您會聽到片段之間的跳躍，這就是生成點擊你在這裏。例如，採取一些合理的數字：我以27.5Hz（鋼琴上的A0）的波形開始，以44100Hz數字化，它將如下所示（其中紅色部分是1024個樣本長）：

alt text http://i48.tinypic.com/zim802.png

所以首先我們先從低通40Hz開始。因此，由於原始頻率小於40Hz，截止頻率爲40Hz的低通濾波器應該不會產生任何效果，我們將得到與輸入幾乎完全匹配的輸出。對？ 錯誤，錯誤，錯誤 - 這基本上是你的問題的核心。的問題是，對於短的部分的想法27.5的赫茲沒有明確的規定，且不能在DFT來表示良好。

即27.5赫茲不在短段特別有意義可以通過在下面的圖看DFT可以看出。請注意，雖然較長段的DFT（黑點）在27.5 Hz處顯示峯值，但短點（紅點）不會。

alt text http://i50.tinypic.com/14w6luw.png

顯然，然後在下方40Hz的過濾，將只捕獲DC偏移和40Hz的低通濾波器的結果在下面的綠色示出。

alt text http://i48.tinypic.com/2vao21w.png

藍色曲線（用200赫茲的截止截取）的開始匹配好得多。但請注意，並不是讓它匹配得好的低頻，而是包含了高頻。直到我們在短段中包含每個可能的頻率，高達22KHz，我們才能最終得到原始正弦波的良好表示。

所有這一切的原因是27.5 Hz正弦波的一小部分是而不是是27.5 Hz正弦波，而它的DFT與27.5 Hz沒有太大關係。