25 Nov
2009
25 Nov
'09
1:52 p.m.
Hi,
I'm trying to do a simple effect in the frequency domain, nothing extraordinary.
I need to apply various gains to specific frequencies.
For this purpose, I want to process a 16bit audio input stream, block by block,
apply a hanning window, do a fft, now apply my effect, and then an inverse fft,
followed by overlap and add buffering.
For now, I'm simply trying to do all of this without applying any effect, just
to see if my time to frequency and back conversion is correct.
I managed to do this in Octave. Now, I'm trying to do it in C with fftw and the
sound I obtain is like "shivering". It's recognizable but severely altered.
I've
spent hours on this, it just won't work.
Attached :
- freqfilter.m: my working octave version
- transform.h and transform.c: time -> frequency -> time routines
- transformtest.c: test program
Could anyone tell me what's wrong with this C code?
Thanks in advance
--
Olivier
25 Nov
25 Nov
4:36 p.m.
On Wed, Nov 25, 2009 at 06:52:00PM +0100, Olivier Guilyardi wrote:
I managed to do this in Octave. Now, I'm trying to
do it in C with fftw and the
sound I obtain is like "shivering". It's recognizable but severely altered.
I've
spent hours on this, it just won't work.
Clearly there's something wrong with the input/output/overlap
logic. It's difficult to say exactly what, because you are
making this much more complicated than it should be.
You are using at least 4 variables to control this: begun,
rindex, windex, rlen. They will be either redundant or in
conflict.
Hint: since the basic processing proceeds in blocks of
fftsize/2, both input and output should use that step
as well. Doing the input at twice that size will just
complicate things, and is probably the deeper cause of
whatever error you made.
It should be easy to modify your code into the form
suggested below.
All these buffers are of size fftsize/2:
IP input save buffer
W1, W2 the first and second half of the window
R1, R2 the first and second half of the real
buffer used for the FFT/IFFT
OP output save buffer
IP and OP are per channel, the others can be shared.
init()
{
IP <= 0
OP <= 0
}
process (input, output)
{
// both input, output are of size fftsize/2
1. R1 <= IP * W1
2. R2 <= input * W2
3. IP <= input
4. FFT, callback, IFFT
5. output <= (R1 + OP) * scale
6. OP <= R2
}
No state variables are required at all.
There are some tricks to avoid the copies at
steps 3 and 6, but in general it's not worth
the trouble.
Ciao,
--
FA
Io lo dico sempre: l'Italia è troppo stretta e lunga.
26 Nov
26 Nov
6:30 p.m.
Hi Fons :)
On 11/25/2009 09:32 PM, fons(a)kokkinizita.net wrote:
process (input, output)
{
// both input, output are of size fftsize/2
1. R1 <= IP * W1
2. R2 <= input * W2
3. IP <= input
4. FFT, callback, IFFT
5. output <= (R1 + OP) * scale
6. OP <= R2
}
Now, that works ! I'm attaching the fixed code for the records.
I've tried hard to understand where the problem was, and one thing I'm sure
about is that you fft'ize a given block of input twice, once with the ascending
part of the window applied (in R1), then again the same block of data multiplied
with the descending part of the window (in R2).
I'm still far from properly understanding these matters and I was indeed not
doing this. In my code, a given block of fftsize/2 input frames was either
multiplied with the ascending or the descending half-window, but not by both
subsequently, which obviously resulted in a loss of information.
Now, how I managed to get this running in Octave is a mistery to me ;)
Anyway you method did solve my problem, and I now still have a lot of work,
because I need all of this to run with integers (fixed point fft) and the joys
of overflow..
Thanks!
--
Olivier
7:01 p.m.
On Thu, Nov 26, 2009 at 11:30:27PM +0100, Olivier Guilyardi wrote:
Hi Fons :)
On 11/25/2009 09:32 PM, fons(a)kokkinizita.net wrote:
Thanks, could be useful one day !
process (input, output)
{
// both input, output are of size fftsize/2
1. R1 <= IP * W1
2. R2 <= input * W2
3. IP <= input
4. FFT, callback, IFFT
5. output <= (R1 + OP) * scale
6. OP <= R2
}
Now, that works ! I'm attaching the fixed code for the records. I've tried hard to understand where the problem
was, and one thing I'm sure
about is that you fft'ize a given block of input twice, once with the ascending
part of the window applied (in R1), then again the same block of data multiplied
with the descending part of the window (in R2).
That's correct. But the first time you combine the first half
of the current block with the second half of the previous,
and the second time you combine the second half of the current
block with the first half of the *next*. I think you failed to
do that somehow.
Now, how I managed to get this running in Octave is a
mistery to me ;)
Miracles happen.
Anyway you method did solve my problem, and I now
still have a lot of work,
because I need all of this to run with integers (fixed point fft) and the joys
of overflow..
Integer FFTs can be lots of fun :-)
One trick you may try (if you write the FFT yourself) is to
distribute the 'scale' thing over the internal iterations of
the FFT. Either shift by one bit each iteration in the FFT,
or in the IFFT, or do it each *second* iteration in both.
Which is best depends on the application, and the type of
signal which you process.
Another thing: you are aware that this method is not just
a filter ? It also modulates the signal... How bad it gets
depends on the frequency response you try to obtain.
Ciao,
--
FA
Io lo dico sempre: l'Italia è troppo stretta e lunga.
27 Nov
27 Nov
6:49 a.m.
On 11/26/2009 11:58 PM, fons(a)kokkinizita.net wrote:
I need to read some good book about all of this when I found some time. You know
my profile a bit now I think, I'm basically a long time coder with little
math/signal processing background. What book would you recommend?
I sense that ;) Actually, this code is intended to run on an increasingly
popular Linux distribution: Android. Which currently means arm and no fpu (gcc
-msoft-float).
However, floats are so handy that I'm first going to check how much cpu time the
current code consumes.
I've tried
hard to understand where the problem was, and one thing I'm sure
about is that you fft'ize a given block of input twice, once with the ascending
part of the window applied (in R1), then again the same block of data multiplied
with the descending part of the window (in R2).
That's correct. But the first time you combine the first half
of the current block with the second half of the previous,
and the second time you combine the second half of the current
block with the first half of the *next*. I think you failed to
do that somehow. Anyway you
method did solve my problem, and I now still have a lot of work,
because I need all of this to run with integers (fixed point fft) and the joys
of overflow..
Integer FFTs can be lots of fun :-) One trick you may try (if you write the FFT yourself)
is to
distribute the 'scale' thing over the internal iterations of
the FFT. Either shift by one bit each iteration in the FFT,
or in the IFFT, or do it each *second* iteration in both.
Which is best depends on the application, and the type of
signal which you process.
No I want to use kissfft. I don't have the skills nor the time to implement the
fft myself, although I would quite like to have some sort of control on that
bloody scale thing.
Doing an fft+ifft with kissfft basically scales down the frames by fftsize. So
if fftsize is 1024, all samples are scaled down by 10 bits, which got me crazy
when I first tried with int16 (Q15).
Another thing: you are aware that this method is not
just
a filter ? It also modulates the signal... How bad it gets
depends on the frequency response you try to obtain.
My usual frequency response measure instrument (hear ;) is saying that with the
current code, the output is very similar to the input. So that should be ok.
That said, do you have any recommendation about the fftsize?
Basically, what I want to do is something like spectral noise gating, similar to
what's implemented in Audacity:
http://wiki.audacityteam.org/index.php?title=Noise_Removal
I'm however also interested in so-called spectral subtraction. I tend to think
that although it may be simpler, it won't give such good results as the gating,
which are quite amazing to me, when applied to speech.
--
Olivier
5411
days inactive
5413
days old
linux-audio-dev@lists.linuxaudio.org
4 comments
2 participants
participants (2)
-
fons＠kokkinizita.net -
Olivier Guilyardi