FIR filtering using window function: implementation problem
I want to make a FIR filter using a window function. I have some sample data and size variable is a count of samples. The windowSize variable is a size of the window function. At first I create the window function (blackman window): the variable window Then I need to multiply it by sin(x) / x function and convolve with real data (variable data):
for (int i = 0; i < size; ++i) {
for (j = 0; j < windowSize; ++j) {
double arg = 2.0 * PI * ((double)j - (double)windowSize / 2.0) / (double)windowSize;
if (i + j - windowSize / 2 < 0)
continue;
if (arg == 0) {
filteredData[i] += data[i + j - windowSize / 2] * window[j] * 1.0 / (double)windowSize;
} else
filteredData[i] += data[i + j - windowSize / 2] * window[j] * (sin(arg) / arg) / (double)windowSize;
}
}
The problem:
As a result I get a filtered data with an average which very different than the average of the original data. Where is a mistake?
In the DSP book it is written that in order to make a FIR filter we should multiply the function sin(x) / x by a window function and then perform the convolution, but nothing 开发者_StackOverflow中文版is written about x in the sin(x) / x, so I used the:
double arg = 2.0 * PI * ((double)j - (double)windowSize / 2.0) / (double)windowSize;
for the x value, the argument of sine, is it correct?
The sin(x)/x
filter is a lowpass filter. That is, it suppresses all frequencies above a certain cutoff frequency.
If the sampling frequency is Fs
(Hertz) and you want a cutoff frequency of fc
(Hertz), You should be using x = 2*PI*fc/(2*Fs)*n
where n
goes from -N
to +N
and N
is large enough that the sin(x)/x function is close to zero. Don't forget that sin(x)/x is 1 when x is zero.
To maintain the average of the signal you have to normalize the filter coefficients by their sum. I.e., set f_norm[k] = f[k] / sum(f[k], k=...)
That's all I have to say at this point. It seems like you have a lot to learn. I suggest a good book on signal processing.
As far as the implementation is concerned it looks like you need to initialise filteredData[i]
, e.g.
for (int i = 0; i < size; ++i) {
filteredData[i] = 0;
for (j = 0; j < windowSize; ++j) {
...
精彩评论