Why Ideal band pass filter not working as expected?
Here is latest version that produce effect close to the desired
void DeleteFrequencies(short *audioDataBuffer, const int bufferSize, int lowestFrequency, int highestFrequency, int sampleRate )
{
int frequencyInHzPerSample = sampleRate / bufferSize;
/* __________________________
/* ___________ __________________________ filter kernel */
int nOfPointsInFilterKernel = (lowestFrequency / frequencyInHzPerSample) + ( bufferSize - highestFrequency / frequencyInHzPerSample);
U u;
double *RealX = new double[bufferSize];
double *ImmX = new double[bufferSize];
ShortArrayToDoubleArray(audioDataBuffer, RealX, bufferSize);
// padd with zeroes, so that inputSignalSamplesNumber + kernelLength - 1 = bufferSize
// convert to frequency domain
ForwardRealFFT(RealX, ImmX, bufferSize);
// cut frequences < 300 && > 3400
int Multiplyer = 1;
for (int i = 0; i < 512; ++i)
{
if (i * 8000 / 1024 > 3400 || i * 8000 / bufferSize < 300 )
{
RealX[i] = 0;
ImmX[i] = 0;
}
if (i < lowestFrequency / frequencyInHzPerSample || i > highestFrequency / frequencyInHzPerSample )
Multiplyer = 0;
else
Multiplyer = 1;
RealX[i] = RealX[i] * Multiplyer /*ReH[f]*/ - ImmX[i] * Multiplyer;
ImmX[i] = ImmX[i] * Multiplyer + RealX[i] * Multiplyer;
}
ReverseRealFFT(RealX, ImmX, bufferSize);
DoubleArrayToShortArray(RealX, audioDataBuffer, bufferSize);
delete [] RealX;
delete [] ImmX;
}
but why it works this way???
Important that I just started learning DSP , so I can be unaware of some important ideas (i appologise for that, but I have task which I need to solve: i need to reduce background noise in the recorder speeach, I try to approach that by cuting off from recorded speech frequencies in ranges <300 && > 3700 (as human voice in [300;3700] range) I started from that method as it is simple, but I found out - it can`t be applied (please see - https://dsp.stackexchange.com/questions/6220/why-is-it-a-bad-idea-to-filter-by-zeroing-out-fft-bins/6224#6224 - thanks to @SleuthEye for reference).
So can you please suggest me simple solution based on the FFT usage that will allow me at least remove given ranges of frequneces ?
I am trying to implement ideal band pass filter. But it isn't working as I expect - only high frequencies are cut.
Here is my implementation description:
union U
{
char ch[2];
short sh;
};
std::fstream in;
std::fstream out;
short audioDataBuffer[1024];
in.open ("mySound.pcm", std::ios::in | std::ios::binary);
out.open("mySoundFilteres.pcm", std::ios::out | std::ios::binary);
int i = 0;
bool isDataInBuffer = true;
U u;
while (in.good())
{
int j = 0;
for (int i = 0; i < 1024 * 2; i+=2)
{
if (false == in.good() && j < 1024) // padd with zeroes
{
audioDataBuffer[j] = 0;
}
in.read((char*)&audioDataBuffer[j], 2);
cout << audioDataBuffer[j];
++j;
}
// Algorithm
double RealX [1024] = {0};
double ImmX [1024] = {0};
ShortArrayToDoubleArray(audioDataBuffer, RealX, 1024);
// convert to frequency domain
ForwardRealFFT(RealX, ImmX, 1024);
// cut frequences < 300 && > 3400
for (int i = 0; i < 512; ++i)
{
if (i * 8000 / 1024 > 3400 || i * 8000 / 1024 < 300 )
{
RealX[i] = 0;
ImmX[i] = 0;
}
}
ReverseRealFFT(RealX, ImmX, 1024);
DoubleArrayToShortArray(RealX, audioDataBuffer, 1024);
for (int i = 0; i < 1024; ++i) // 7 6 5 4 3 2 1 0 - byte order hence we write ch[1] then ch[0]
{
u.sh = audioDataBuffer[i];
out.write(&u.ch[1], 1);
out.write(&u.ch[0], 1);
}
}
in.close();
out.close();
when I write result to a file, open it audacity and check spectr analysis, and see that high frequences are cut, but low still remains (they starts from 0)
What I am doing wrong?
Here is sound frequency spectr before 
Here is sound frequency after I zeroed needed values 
Please help!
Update:
Here is code I came up with, what I should padd with Zeroes???
void DeleteFrequencies(short *audioDataBuffer, const int bufferSize, int lowestFrequency, int highestFrequency, int sampleRate )
{
// FFT must be the same length as output segment - to prevent circular convultion
//
int frequencyInHzPerSample = sampleRate / bufferSize;
/* __________________________
/* ___________ __________________________ filter kernel */
int nOfPointsInFilterKernel = (lowestFrequency / frequencyInHzPerSample) + ( bufferSize - highestFrequency / frequencyInHzPerSample);
U u;
double *RealX = new double[bufferSize];
double *ImmX = new double[bufferSize];
ShortArrayToDoubleArray(audioDataBuffer, RealX, bufferSize);
// padd with zeroes, so that inputSignalSamplesNumber + kernelLength - 1 = bufferSize
// convert to frequency domain
ForwardRealFFT(RealX, ImmX, bufferSize);
// cut frequences < 300 && > 3400
int Multiplyer = 1;
for (int i = 0; i < 512; ++i)
{
/*if (i * 8000 / 1024 > 3400 || i * 8000 / bufferSize < 300 )
{
RealX[i] = 0;
ImmX[i] = 0;
}*/
if (i < lowestFrequency / frequencyInHzPerSample || i > highestFrequency / frequencyInHzPerSample )
Multiplyer = 0;
else
Multiplyer = 1;
RealX[i] = RealX[i] * Multiplyer /*ReH[f]*/ - ImmX[i] * Multiplyer;
ImmX[i] = ImmX[i] * Multiplyer + RealX[i] * Multiplyer;
}
ReverseRealFFT(RealX, ImmX, bufferSize);
DoubleArrayToShortArray(RealX, audioDataBuffer, bufferSize);
delete [] RealX;
delete [] ImmX;
}
it produce the following spectrum (low frequencies are cut, but high not) 
void ForwardRealFFT(double* RealX, double* ImmX, int nOfSamples)
{
short nh, i, j, nMinus1, nDiv2, nDiv4Minus1, im, ip, ip2, ipm, nOfCompositionSteps, LE, LE2, jm1;
double ur, ui, sr, si, tr, ti;
// Step 1 : separate even from odd points
nh = nOfSamples / 2 - 1;
for (i = 0; i <= nh; ++i)
{
RealX[i] = RealX[2*i];
ImmX[i] = RealX[2*i + 1];
}
// Step 2: calculate nOfSamples/2 points using complex FFT
// advantage in efficiency, as nOfSamples/2 requires 1/2 of the time as nOfSamples point FFT
nOfSamples /= 2;
ForwardDiscreteFT(RealX, ImmX, nOfSamples );
nOfSamples *= 2;
// Step 3: even/odd frequency domain decomposition
nMinus1 = nOfSamples - 1;
nDiv2 = nOfSamples / 2;
nDiv4Minus1 = nOfSamples / 4 - 1;
for (i = 1; i <= nDiv4Minus1; ++i)
{
im = nDiv2 - i;
ip2 = i + nDiv2;
ipm = im + nDiv2;
RealX[ip2] = (ImmX[i] + ImmX[im]) / 2;
RealX[ipm] = RealX[ip2];
ImmX[ip2] = -(RealX[i] - RealX[im]) / 2;
ImmX[ipm] = - ImmX[ip2];
RealX[i] = (RealX[i] + RealX[im]) / 2;
RealX[im] = RealX[i];
ImmX[i] = (ImmX[i] - ImmX[im]) / 2;
ImmX[im] = - ImmX[i];
}
RealX[nOfSamples * 3 / 4] = ImmX[nOfSamples / 4];
RealX[nDiv2] = ImmX[0];
ImmX[nOfSamples * 3 / 4] = 0;
ImmX[nDiv2] = 0;
ImmX[nOfSamples / 4] = 0;
ImmX[0] = 0;
// 3-rd step: combine the nOfSamples frequency spectra in the exact reverse order
// that the time domain decomposition took place
nOfCompositionSteps = log((double)nOfSamples) / log(2.0);
LE = pow(2.0,nOfCompositionSteps);
LE2 = LE / 2;
ur = 1;
ui = 0;
sr = cos(M_PI/LE2);
si = -sin(M_PI/LE2);
for (j = 1; j <= LE2; ++j)
{
jm1 = j - 1;
for (i = jm1; i <= nMinus1; i += LE)
{
ip = i + LE2;
tr = RealX[ip] * ur - ImmX[ip] * ui;
ti = RealX[ip] * ui + ImmX[ip] * ur;
RealX[ip] = RealX[i] - tr;
ImmX[ip] = ImmX[i] - ti;
RealX[i] = RealX[i] + tr;
ImmX[i] = ImmX[i] + ti;
}
tr = ur;
ur = tr * sr - ui * si;
ui = tr * si + ui * sr;
}
}
You may want to have a look at this answer for an explanation for the effects you are observing.
Otherwise, the 'ideal' filter you are trying to achieve is more a mathematical tool than a practical implementation since the rectangular function in the frequency domain (with a zero-transition and infinite stopband attenuation) corresponds to an infinite length impulse response in the time domain.
To obtain a more practical filter, you must first define desired filter characteristics such as transition width and stopband attenuation based on your specific application needs. Based on these specifications, the filter coefficients can be derived using one of various filter design methods such as:
Perhaps the closest to what you're doing is the Window method. Using that method, something as simple as a triangular window can help increase the stopband attenuation, but you may want to experiment with other window choices (many available from the same link). Increasing the window length would help reduce the transition width.
Once you have completed your filter design, you can apply the filter in the frequency-domain using the overlap-add method or the overlap-save method. Using either of these methods, you would split your input signal in chunks of length L, and pad to some convenient size N >= L+M-1, where M is the number of filter coefficients (for example if you have a filter with 42 coefficients, you might choose N = 128, from which L = N-M+1 = 87).
使用FFT / IFFT进行快速卷积滤波需要零填充至少为滤波器长度的两倍（出于性能原因，通常为2的下幂），然后使用重叠添加或重叠保存方法来消除循环卷积伪像。
After doing the real FFT you get your spectral data twice: Once in the bins from 0 to 512, and a mirror spectrum in the bins 513 to 1024. Your code however only clears the lower spectrum.
Try this:
for (int i = 0; i < 512; ++i)
{
if (i * 8000 / 1024 > 3400 || i * 8000 / 1024 < 300 )
{
RealX[i] = 0;
ImmX[i] = 0;
// clear mirror spectrum as well:
RealX[1023-i] = 0;
ImmX[1023-i] = 0;
}
}
That may help unless your FFT implementation does this step automatically.
Btw, just zeroing out frequency bins like you did is not a good way to do such filters. Expect a very nasty phase response and a lot of ringing in your signal.
链接地址: http://www.djcxy.com/p/33842.html上一篇: 在Javascript中使用FFT计算音频文件的平均幅度
下一篇: 为什么理想带通滤波器不能按预期工作？