Jump to content
  • Advertisement

Archived

This topic is now archived and is closed to further replies.

Crispy

noise from a low-pass filter

This topic is 5274 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

If you intended to correct an error in the post then please contact us.

Recommended Posts

When doing FFT convolution (multiplying the frequency spectra of a signal and a lowpass filter), the lower the cutoff frequency is moved, the more noise (in the form of crackles) the convolution produces. If fc = 0.5, eg the filter is an all-pass filter, the sound is left unchanged. Whenever fc becomes other than 0 or 0.5, very strong noise becomes an issue. To provide at least some insight, here's the general outline of my code:
size = some_power_of_two;

FillWithSignal(&signal, size);
CreateFilter(&filter, cutoff, size);

FFT(signal, size);
FFT(filter, size);

for(int i = 0; i < size * 2; i++)
  signal[i] *= filter[i];

////////////////////////////////////////////////

//the same effect is produced by simply doing://

//                                            //

//for(int i = 50; i < size * 2; i++)          //

//  signal[i] = 0;                            //

////////////////////////////////////////////////


IFFT(signal, size);

for(int i = 0; i < size * 2; i++)
  signal[i] /= size;
Can anyone pinpoint a mistake or put forth a suggestion? PS: and, yes - FillWithSignal() and CreateFilter() function properly. Cheers edit: stoopid forum parser...

TCP - Transmission by Carrier Pigeons
   
[edited by - crispy on April 13, 2004 5:54:01 PM]

Share this post


Link to post
Share on other sites
Advertisement
What do you mean by crackles? Is this some kind of block based filtering for a continuous signal or do you filter a single finite singal in one FFT/mult/IFFT? Are you sure FFT/IFFT works correct?

What are you actually trying to do?

Share this post


Link to post
Share on other sites
>> What do you mean by crackles?

A "pounding repetitive noise" as if every n''th sample came out way louder than all all other samples. It sounds like crackling in the directest of senses.

>> Is this some kind of block based filtering for a continuous signal or do you filter a single finite singal in one FFT/mult/IFFT?

I''m filtering by blocks of 512 samples right now, taken from a continuously generated signal. Unless you''re saying I need to use the overlap-add in this case, I don''t really understand the question. The filtering is done in real-time. If I do need to use the overlap-add method, I''d appreciate some insight as to why.

>> Are you sure FFT/IFFT works correct?

Yes. FFT output corresponds to what I expected and so does IFFT (the priginal signal with mild alterations, which I write off to be round-off noise).

>> What are you actually trying to do?

Apply a low-pass filter via FFT convolution.




TCP - Transmission by Carrier Pigeons

Share this post


Link to post
Share on other sites
quote:
Original post by Crispy
>> What do you mean by crackles?

A "pounding repetitive noise" as if every n''th sample came out way louder than all all other samples. It sounds like crackling in the directest of senses.


I could almost guess it was something like that. You''re not alone having a similar problem when filtering a time signal in non-time domain.

quote:
Original post by Crispy
>> Is this some kind of block based filtering for a continuous signal or do you filter a single finite singal in one FFT/mult/IFFT?

I''m filtering by blocks of 512 samples right now, taken from a continuously generated signal. Unless you''re saying I need to use the overlap-add in this case, I don''t really understand the question. The filtering is done in real-time. If I do need to use the overlap-add method, I''d appreciate some insight as to why.


You understood it correct. That''s a block based filtering of a continuous signal. The other thing I mentioned was if you had a, say, prerecorded N-sample (N is quite large) signal, and performed the convolution offline using N-point FFT/mult/IFFT.

quote:
Original post by Crispy
>> Are you sure FFT/IFFT works correct?

Yes. FFT output corresponds to what I expected and so does IFFT (the priginal signal with mild alterations, which I write off to be round-off noise).


Sound good to me.

quote:
Original post by Crispy
>> What are you actually trying to do?

Apply a low-pass filter via FFT convolution.


Ok, asked just in case this was a part of a bigger project or so.


I think I know what''s causing the noise. If you can save the signal and watch it offline, you will probably see that there are sudden jumps in the signal every N sample, where N is the block length. More specifically, i happens exactly between two blocks. It''s because the signal in each block does not line up perfectly when going from one block to the following block.

I think the overlap-add method, which you mentioned, would work here. What you want is linear convolution. In time domain that is the same as performing a convolution with an FIR-filter or whatever you use, and just let the signal pass through it without interruption.

You''ve probably heard about the saying convolution in time domain is multiplication in frequency domain, or a variant of it. This is of course true, but the result of multiplication in frequency domain is a circular convolution in time domain, not a linear convolution, and that''s bad in this case. The blockwise circular convolution instead of linear convolution (which involves overlapping blocks) is what causes the discontinuities in the signal between the blocks.

I think I skip the explanation of the overlap-add method. It sounds too much like a schoolwork assignment so you have to do the reasearch on it yourself, but if you have any problems, feel free to ask for pointers. I don''t have any deep practical experience in the overlap-add method of filtering, but I will try not to say too much I have no clue about

Share this post


Link to post
Share on other sites
>> I could almost guess it was something like that. You''re not alone having a similar problem when filtering a time signal in non-time domain.

Figures .

>> Ok, asked just in case this was a part of a bigger project or so.

Well, at the moment it''s not. And this particular thing isn''t even related to my school stuff (although I am currently writing a thesis on DSP and filtering, but I mainly deal with images and PSF convolution in the spatial domain). Nevertheless, I don''t need to write this stuff for reasons related to schoolwork, if that''s your worry - the reason I''m doing this is because it kinda interests me and I''ve been wanting to write a proper real-time sudio signal filter application for 2-3 years now. Right now I feel I have the background I always lacked, but a few things aren''t going as smoothly as I''d hoped for, things get forgotten and overall progress is slow due to time constraints.

>> I think I know what''s causing the noise. If you can save the signal and watch it offline, you will probably see that there are sudden jumps in the signal every N sample, where N is the block length. More specifically, i happens exactly between two blocks. It''s because the signal in each block does not line up perfectly when going from one block to the following block.

That''s probably it - as I said the noise is repetitive and periodic. I went over the overlap-add method again (oh the things you read about and forget about just as quickly!) and here''s the basic outline as to how it''s supposed to work (for me it doesn''t work - in fact for me the output doesn''t change much at all):


filter = CreateFilter(N_samples_long, append_N_samples_of_zeroes_to_end);
signal = GetSignal(N_samples_long, append_N_samples_of_zeroes_to_end);

FFT(filter, N * 2);
FFT(signal, N * 2);

for(int i = 0; i < N; i++)
{
//handle both real and imaginary parts

Multiply(signal[i], filter[i]);
}

IFFT(signal, N * 2);

//overlap first N samples from previous block

for(int i = 0; i < N; i++)
{
//real and imaginary parts are interleaved in both arrays,

//so, in rectangular co-ordinates, only N / 2 samples are added

signal[i] += overlap[i];
}

//store last N samples for next block

for(int i = N, j = 0; i < N * 2; i++, j++)
{
//same deal with the interleaved samples

overlap[j] = signal[i];
}


The reason there''s a Multiply() function to multiply the signals is to denote polar multiplication of the two components (s* for signal, f* for filter, *re for real, *im for imaginary):


tmp = sre[k] * fre[k] - sim[k] * fim[k];
sim[k] = sre[k] * fre[k] + sim[k] * fim[k];
sre[k] = tmp;


Anyway, I can''t get the darn thing to work properly, although the whole deal seems very straightforward. I suppose you have more experience with these things so maybe you could suggest something?

PS: cheers for your feedback so far!




TCP - Transmission by Carrier Pigeons

Share this post


Link to post
Share on other sites
I''m a little bit confused about what is what your code. You use the same size for the block length and overlap length, and you use the same arrays for both time and frequency domain data, and store both imaginary and real data in them. Usually the block length is larger than the overlap. The overlap is the length of the filter in time domain. But as I understand it, the general idea is correct.

I made a MATLAB-script for a quick look at how well the overlap-add method works, and it works great. If you have access to MATLAB you can try it and step through it to see what happens. If you know MATLAB, it shouldn''t be too difficult to translate it into some other language.

function overlap_add;

M = 21; % Filterlength (also one more than the amount of overlap)
L = 512; % Block length
B = 10; % Number of blocks of testdata
fc = 0.2; % Cutoff frequency

% Adjust block length to make the size of the FFT a power of two.
L = L-M+1;

% Generate random input data.
x = randn(L*B, 1);

% Create a lowpass filter.
% H should be an L+M-1-tap DFT of h, which is an M-tap filter!
h = fc * sinc(fc * [(-M+1)/2: (M-1)/2]'');
H = fft(h, L+M-1);

% Fill y with M-1 zeros for first overlap.
y = zeros(M-1, 1);

I = [1: L]'';
J = [1: L+M-1]'';
K = [1: M-1]'';

for n = 0: L: L*B-1

% Get block from input, append zeros and transform.
% Multiply by the filter to get the output.
X = fft(x(I), L+M-1);
Y = X .* H;

% Append calculated output almost at the end of y,
% adding the two where they overlap.
y(J) = [y(K); zeros(L, 1)] + real(ifft(Y));

I = I + L;
J = J + L;
K = K + L;

end

% Remove last M-1 samples from output buffer to
% make y exactly B blocks long.
y = y(1: L*B);

% As reference, filter the input using time-domain convolution.
yy = filter(h, 1, x);

n = [0: L*B-1]'';
plot(n, y, ''b-'', n, yy, ''r-'');
axis([0, L*B-1 -1.5 1.5]);

return;

Share this post


Link to post
Share on other sites
Well, a few of the details in your MatLab script are a bit vague since I''ve never used MatLab for any constructive purposes and can therefore only rely on my understanding of general programming concepts. Besides, I don''t own MatLab, so I can''t test it...

I''ll just bring a numeric example because this just refuses to add up (probably just goes to prove how much of an idiot I am ):

I''m currently using 512 sample blocks and a 512 sample filter kernel, which I''ve substituted with a perfect hard-coded lowpass filter (everything above sample N is simply zeroed).

1) the 512 sample data block is interlaced with zeroes, bringing the total length of the array to 1024 units
2) Taking a 512 sample FFT of it produces a 1024 unit frequency spectrum (one sample translates to one real and one imaginary value)
3) The two spectra are multiplied (everything above N Hz zeroed)
4) A 512 sample IFFT is taken of the result
5) The first 256 samples (256 real + 256 imaginary values = 512 units of interlaced data*) are added to the result from the overlap array, which is initially zeroed
6) The last 256 samples (256 real + 256 imaginary values) of the result are copied to the overlap array to be used with the next block

* the imaginary part is quite useless at this point, but it''s just slighly simpler to code and keeps all numbers in the same ballpark

It works as it is, but the result isn''t at all that different from circular convolution.

Anyway, I''m pretty confused about the number of samples I have to use and where to use how many. Common sense tells me that if I zero everything above, say, 50 samples, the last 512 samples are bound to be all zeroes, which means the overlap for the next block will be all zeroes, which in turn means the whole overlapping doesn''t change a thing. Damn - it''s weird how a few numbers can confuse you out of your wits and 3 years of university strikes you as crystal clear...




TCP - Transmission by Carrier Pigeons

Share this post


Link to post
Share on other sites
Ok, I will try to describe the steps using the numbers you use. Unless I misunderstood you, pretty much all of the steps you posted, except the first, are wrong in one way or another

First, precalculate the DFT of the filter kernel. No need to do it every frame at least.

  1. Add 511 zeros after the filter kernel and perform DFT. What you get is a 1023-tap DFT.



For each block.

  1. Extract 512 samples from the input signal.

  2. Add 511 zeros after the block and perform 1023-tap DFT.

  3. Multiply the 1023 complex values of the signal block and filter.

  4. Perform 1023-tap IDFT and get 1023 time samples as new output.

  5. Add 511 overlap samples from previous block to the first 511 samples of the output.

  6. Save last 511 samples for next block.


As you can see, with the values you use, the size of the DFT/IDFT are not a power of two, so you can''t really use the FFT. Either change the block length or the size of the filter kernel to make it a power of two. I suggest you increase the block length to 513 samples.

Oh, and for the first block, the overlap from "previous" block should of course be all zeros.

Share this post


Link to post
Share on other sites
DAMMMMMMMMMMMIIIIIIIIIIIIIIIIIIITTTTTTT. I'm about to go on a violent rampage right now. REALLY REALLY violent.

Well, the problem turned out to be the following (if you can spot it):

for(int k = 0; k < 2048; k++)
{
float temp = signal_fft[k] * filter_fft[k] - signal_fft[k + 1] * filter_fft[k + 1];
signal_fft[k + 1] = signal_fft[k] * filter_fft[k + 1] + signal_fft[k + 1] * filter_fft[k];
signal_fft[k] = temp;
}

Anyway, I need something sweet to eat (or, alternatively, do something suicidal) because I've been breaking my skull over this for 4-5 days now.

edit: in all the anger towards myself I completely forgot about thanking you for being a babysitter and talking me through this . Cheers!





TCP - Transmission by Carrier Pigeons



[edited by - crispy on April 16, 2004 12:30:14 PM]

Share this post


Link to post
Share on other sites

  • Advertisement
×

Important Information

By using GameDev.net, you agree to our community Guidelines, Terms of Use, and Privacy Policy.

We are the game development community.

Whether you are an indie, hobbyist, AAA developer, or just trying to learn, GameDev.net is the place for you to learn, share, and connect with the games industry. Learn more About Us or sign up!

Sign me up!