【发布时间】：2014-07-29 15:58:46
【问题描述】：
我实现了我的过滤器，其中使用了重叠添加方法来防止循环抽搐。
输入 - 有噪音的文件，输出应该是过滤后的文件。
我的结果： 稍微修改了输出，没有削减频率
我的猜测是我错误地在滤波器内核上的频域输入信号中相乘 （我的意图是切断不在 [300,3700] 范围内的频率）。乘法应该怎么做？
我使用 blackmanwindow 构建内核 - 我的理解正确吗？ （我计算每个滤波器样本的频率量，然后遍历样本，看看它是否在我想要截断的范围内，我使用布莱克曼窗口的公式计算频率。）
我刚开始学习 DSP。
这是我的实现（有什么问题？？？）：
void DeleteFrequencies(char* fileWithNoise, char* resultFile, const int bufferSize, int lowestFrequency, int highestFrequency, int sampleRate )
{
// |1|. files
std::fstream in;
std::fstream out;
in.open (fileWithNoise, std::ios::in | std::ios::binary);
out.open(resultFile, std::ios::out | std::ios::binary);
// |2|. Filter kernel design. I shall use blackman window
// fundamental params
const int filterKernelLength = 200; // 512
const int filterMaxFrequency = sampleRate / 2; // 8000 / 2
const int frequencyPerSamle = filterMaxFrequency / filterKernelLength;
double *RealFilterResp = new double [bufferSize / 2];
double *ImmFilterResp = new double [bufferSize / 2];
// coefficients for Blackman window
const double a0 = 0.42659;
const double a1 = 0.49656;
const double a2 = 0.076849;
// construct filter kernel
for (int i = 0 ; i < bufferSize / 2; ++i)
{
if ( i >= filterKernelLength ) // padd filter kernel with zeroes
{
RealFilterResp[i] = 0;
ImmFilterResp[i] = 0;
}
else if (i * frequencyPerSamle < lowestFrequency || i * frequencyPerSamle > highestFrequency)
{
// apply blackman window (to eleminate frequencies < 300 hz and > 3700 hz)
RealFilterResp[i] = a0 - a1 * cos (2 * M_PI * i / (bufferSize / 2 - 1)) + a2 * cos (4 * M_PI / (bufferSize / 2 - 1));
ImmFilterResp[i] = a0 - a1 * cos (2 * M_PI * i / (bufferSize / 2 - 1)) + a2 * cos (4 * M_PI / (bufferSize / 2 - 1));
}
else
{
RealFilterResp[i] = 1;
ImmFilterResp[i] = 1;
}
}
// |3|. overlap add method
// calculate parameters for overlap add method (we use it to prevent circular convultion)
const int FFT_length = pow (2.0 ,(int)(log(bufferSize + filterKernelLength - 1.0)/log(2.0)) + 1.0);
double *OLAP = new double[bufferSize / 2 ]; // holds the overlapping samples from segment to segment
memset(OLAP,0, bufferSize / 2 * sizeof (double));
double *RealX = new double[bufferSize];
memset(RealX, 0, bufferSize * sizeof(double));
double *ImmX = new double[bufferSize];
memset(ImmX, 0, bufferSize * sizeof(double));
short* audioDataBuffer = new short[bufferSize];
memset(audioDataBuffer, 0 , sizeof(short) * bufferSize);
// start reading from file by chunks of bufferSize
while (in.good())
{
// get proper chunk of data
FillBufferFromFile(audioDataBuffer, bufferSize, in); // read chunk from file
ShortArrayToDoubleArray(audioDataBuffer, RealX, bufferSize); // fill RealPart
ForwardRealFFT(RealX, ImmX, bufferSize); // go to frequency domain
// perform convultion as multiplication in frequency domain
for (int j = 0; j < bufferSize / 2; ++j)
{
double tmp = RealX[j] * RealFilterResp[j] - ImmX[j] * ImmFilterResp[j];
ImmX[j] = RealX[j] * ImmFilterResp[j] + ImmX[j] * RealFilterResp[j];
RealX[j] = tmp;
}
// Inverse FFT
ReverseRealFFT(RealX, ImmX, bufferSize); // go to time domain
// add last segment overlap to this segment
for (int j = 0; j < filterKernelLength - 2; ++j )
{
RealX[j] += OLAP[j];
}
// save samples that will overlap the next segment
for (int j = bufferSize/2 + 1; j < bufferSize; ++j )
{
OLAP[j - bufferSize/2 - 1] = audioDataBuffer[j];
}
// write results
DoubleArrayToShortArray(RealX, audioDataBuffer, bufferSize);
FillFileFromBuffer(audioDataBuffer, bufferSize, out);
}
/*ReverseRealFFT(RealX, ImmX, bufferSize
);
DoubleArrayToShortArray(RealX, audioDataBuffer, bufferSize);*/
delete [] audioDataBuffer;
delete [] RealFilterResp;
delete [] ImmFilterResp;
delete [] OLAP;
delete [] RealX;
delete [] ImmX;
in.close();
out.close();
}
【问题讨论】：
-
'我的猜测是我乘错了......'你不能用调试器消除你的疑虑吗？
-
@ πάντα ῥεῖ - 我认为存在逻辑错误，调试器无法帮助我发现，我已经连续调试了 11 个小时的代码，之前已经检查过数千次发布到 SO。
标签： c++ audio filtering signal-processing fft