學期開始時，我們希望能實作出使用 EQ 時常見的六種 Filter， 因此便開始研究音訊處理的基本原理， 過程中，也認識到 Filter 分為兩種類型：FIR 與 IIR 由其各自的特性，在做 DSP 時，通常使用 IIR 會更有效率，其中最常見的是 Biquad Filter。 但經過老師的提點，我們發現在現今的設備上實現高階的 FIR filter 是非常 OK 的，
因此，我們最後規劃以：
- FIR 實現 LPF、HPF
- IIR 實現 Notch、Peak、LSF、HSF
- 並期待能提供 GUI，畫出訊號的頻譜 & Filter 的頻率響應
剛開始，我們只希望能先看到一些成果，因此，查到 這個網站 介紹如何產生 LPF 的 impulse response 時，
/* Construct impulse response of Low Pass Filter (FIR) */
std::array <float, 560> h = { 0.0 };
for (int i = -(order - 1) / 2, j = 0; i <= order / 2; i++, j++) {
if (i != 0)
h[j] = (sin(2 * fc * i * MY_PI) / (2 * fc * i * MY_PI));
else
h[j] = (1.0);
}
/* Set blackman window */
w.clear();
for (int i = 0; i < order; i++)
w.push_back(0.42 - 0.5 * cos(2 * MY_PI * i / (order - 1)) + 0.08 * os(4 * MY_PI * i / (order - 1)));
/* Windowed-Sinc Filter */
for (int i = 0; i < order; i++) {
h[i] = h[i] * w.at(i);
impulse_response_sum += h[i];
}
/* Normalized windowed-sinc filter */
for (int i = 0; i < order; i++)
h[i] /= impulse_response_sum;
我們就直接將
因此，再往下研究之後，辣個男人 出現了。
下面是針對某段訊號進行濾波後的結果， 可以看見，比起直接在時域作 convolution， 先將 Signal 及 Filter 傅立葉轉換至頻域，做點乘再逆轉換回來的運算速度，快了將近 25 倍。 Image Source
時間複雜度分析 ↓
/* signal FFT */
forwardFFT.performRealOnlyForwardTransform(&fftData.at(0));
/* Dot product with the frequency response of the filter */
for (int i = 0; i < fftSize * 2; i += 2)
{
float sr = fftData[i];
float si = fftData[i + 1];
float fr = FIR_freq_response_c[i];
float fi = FIR_freq_response_c[i + 1];
fftData[i] = sr * fr - si * fi;
fftData[i + 1] = sr * fi + si * fr;
}
/* signal iFFT */
forwardFFT.performRealOnlyInverseTransform(&fftData.at(0));
/* filter FFT */
for (int i = 0; i < fftSize * 2; i++)
{
H[i] = (i < order) ? h[i] : 0; // 將 filter 0-padding 至 fft size
}
forwardFFT.performRealOnlyForwardTransform(&H.at(0));
:::warning 但做完之後，聲音聽起來還是怪怪的? :::
一個長度
因此，每個 Process Block 做完 Filtering 都會多出一段不完整的輸出， 下次的結果須再加上此段訊號。
動圖為例 GIF Source ↓
/* Overlap-Add method */
L = buffer.getNumSamples(); // process block size
M = numTap; // filter size
for (int i = 0; i < M - 1; i++)
{
fftData[i] += overlap[i];
}
for (int i = 0; i<buffer.getNumSamples(); i++) {
buffer.setSample(0, i, fftData[i]);
buffer.setSample(1, i, fftData[i]);
}
for (int i = 0; i < M - 1; i++)
{
overlap[i] = fftData[L + i];
}
在修正這些訊號之後，聲音聽起來就比較沒有雜音了！
這個部分，我們使用 juce 內建的 IIR API，取出特定 filter 的
/* initialize the first few input and output signal to 0 */
x = {0, 0, 0}; //input signal
y = {0, 0}; //output signal
/* construct impulse response */
coef = juce::dsp::IIR::ArrayCoefficients<float>::makeNotch(getSampleRate(), f, Q);
h1 = { coef[2] / coef[3], coef[1] / coef[3], coef[0] / coef[3] };
h2 = { coef[5] / coef[3], coef[4] / coef[3] };
/* apply IIR filter */
for (int i = 0; i < buffer.getNumSamples(); ++i) {
x.erase(x.begin());
x.push_back(buffer.getSample(0, i));
value = 0.0;
for (int j = 0; j < 3; j++)
value += x.at(j) * h1.at(j);
for (int j = 0; j < 2; j++)
value -= y.at(j) * h2.at(j);
y.erase(y.begin());
y.push_back(value);
buffer.setSample(0, i, value);
buffer.setSample(1, i, value);
}
- 串接多個 IIR filters 時會產生雜音，目前尚不確定雜音產生的原因
- 首先，我們先看到套用多個 IIR filters 的計算方法，如下：
void EQAudioProcessor::applyIIRFilter(juce::AudioBuffer<float> &buffer) { float value; std::vector <float> x; // input signal std::vector <float> y; // output signal std::vector <float> h1; // impulse response for input std::vector <float> h2; // impulse response for output for (int j = 0; j < coefs.size(); ++j) { /* initialize the first few input and output signal to 0 */ x = {0, 0, 0}; y = {0, 0}; /* construct impulse response */ auto coef = coefs[j]; h1 = { coef[2] / coef[3], coef[1] / coef[3], coef[0] / coef[3] }; h2 = { coef[5] / coef[3], coef[4] / coef[3] }; /* apply IIR filter */ for (int i = 0; i < buffer.getNumSamples(); ++i) { x.erase(x.begin()); x.push_back(buffer.getSample(0, i)); value = 0.0; for (int j = 0; j < 3; j++) value += x.at(j) * h1.at(j); for (int j = 0; j < 2; j++) value -= y.at(j) * h2.at(j); y.erase(y.begin()); y.push_back(value); buffer.setSample(0, i, value); buffer.setSample(1, i, value); } } }- 這裡使用
for迴圈依次套用各個 IIR filter
- 這裡使用
- 根據此演算法，我們對於雜音產生的原因有兩種猜想
- 計算時間的問題：若套用 N 個 filters，則計算的時間就會變成 N 倍
- 計算上產生誤差：二階 IIR filter 的前 2 個音訊因為沒有前面的資料作為參考，便將前面的資料直接當作 0 使用，因而產生些許誤差。在套用多個 IIR filters 後，該誤差就放大到可以被人耳辨識出來了
將 Filter 的脈衝響應做 FFT 之後，得到頻率響應。
/* FFT */
for (int i = 0; i < fftSize * 2; i++)
{
H_freq[i] = (i < order) ? h[i] : 0;
}
forwardFFT.performFrequencyOnlyForwardTransform(&H_freq.at(0));
利用 JUCE 提供的 API，給定 IIR 的
juce::dsp::IIR::Coefficients<float> IIR_filter(coef[0], coef[1], coef[2], coef[3], coef[4], coef[5]);
IIR_filter.getMagnitudeForFrequencyArray(&frequencies[0], &response[0], fftSize/2, getSampleRate());
為了待會與 FIR 的頻率響應相結合，我們將 IIR 的頻率響應也切為 fftSize$/2$ 個 bin。
std::vector<double> response;
response.resize(fftSize/2);
juce::FloatVectorOperations::fill(&response[0], 0, fftSize/2);
設定欲求的特定頻率組合：
for (int i = 0; i < fftSize/2; ++i)
{
frequencies.push_back(binToFreq(i,sampleRate));
}
float binToFreq(int binNum, double sr)
{
float binWidth = sr / fftSize;
auto binFreq = binNum * binWidth;
return binFreq;
}
// 0 ~ 511 bin
// ↓
// 0 ~ 24000 Hz
//
// frequencies[] = { 0Hz, 46.875Hz, ... , 24000Hz }*/
首先將剛剛求得的 FIR、IIR 各自的頻率響應結合：
for (int i = 0; i < fftSize / 2; i++)
{
freqResponse[i] = IIR_Response[i] * FIR_freq_response[i];
}
接著我們沿用助教提供的範例程式碼，實作音訊的頻譜和 filter 的響應曲線。
/* filter response curve generator */
void generateResponsePath(const std::vector<double>& renderData,
juce::Rectangle<float> fftBounds,
int fftSize,
float binWidth,
float negativeInfinity)
{
auto left = fftBounds.getX();
auto top = fftBounds.getY();
auto bottom = fftBounds.getHeight();
auto width = fftBounds.getWidth();
int numBins = (int)fftSize/2; // FFT block 的後半段與前半段對稱，僅取前半
PathType p;
p.preallocateSpace(3 * (int)fftBounds.getWidth());
/* 利用 lambda expression 將頻率響應 map 到 GUI 上 */
auto map = [bottom, top, negativeInfinity](float v)
{
return juce::jmap(v,
0.f, 2.5f,
float(bottom), top);
};
auto y = map(renderData[0]);
if (std::isnan(y) || std::isinf(y))
y = bottom;
p.startNewSubPath(left, y);
const int pathResolution = 1;
for (int binNum = 1; binNum < numBins; binNum += pathResolution)
{
y = map(renderData[binNum]);
if (!std::isnan(y) && !std::isinf(y))
{
auto binFreq = binNum * binWidth;
auto normalizedBinX = juce::mapFromLog10(binFreq, 20.f, 20000.f);
int binX = std::floor(left + normalizedBinX * width);
p.lineTo(binX, y);
}
}
pathFifo.push(p);
}
每個 band 的左上角都有 activate 的選項，可以設定該 band 的 filter 是否要套用
- 程式流程設計的問題：目前設定每當旋鈕的數值改變時，filter 就會重新生成，因此當快速轉動旋扭時，會產生大量的計算，進而導致雜音產生