pokey909 wrote:
> Hi!
> 
> I am sampling at 256Hz and I want to filter out all frequencies except
> 1-9Hz. The filter should have the highest attenuation as possible. My
> approach would be to use the FFTW library, transform from time to freq.
> domain and set all bins which are outside of my passband to 0 (like a
> rectangular window).

   ...

You can't hope for a real solution without giving a real spec. "Highest 
possible" is not a real spec. What will you sacrifice to achieve that? 
The most likely factor is time. Show me a good filter that gives results 
in milliseconds, and I'll show you a "better" one that takes seconds. 
Imagine what attenuation I could achieve in a week!

Zeroing all bins outside your desired passband doesn't work. See the 
recent thread "Frequency domain filtering (rectangular window question)" 
for discussion on this. It comes up so often it should be in the FAQ as 
a known trap.

Jerry
-- 
Engineering is the art of making what you want from things you can get.
�����������������������������������������������������������������������

Spend a little time. Read previous (even recent) posts.  I think you
will find it has been covered ad nauseum.

Dirk

pokey909 wrote:
> Hi!
>
> I am sampling at 256Hz and I want to filter out all frequencies except
> 1-9Hz. The filter should have the highest attenuation as possible. My
> approach would be to use the FFTW library, transform from time to freq.
> domain and set all bins which are outside of my passband to 0 (like a
> rectangular window).
> But honestly I dont have a clue how to accomplish that with FFTW. I tried
> a hartley transform, applied my window function and transformed it back,
> but I always get a "half-period-sine-like-shaped" signal. So I did
> obviously something wrong.
>
> I used following code (sketched):
>
> func transform() {
> // Im using the same array here, this is no typo!
> // could this be a problem?
> fftw_create_plan_r2r_1d(sampleNum, smp_array, smp_array, FFTW_DHT,
> FFTW_ESTIMATE);
> execute...
> }
>
> func filter(f1, f2) {
> transform();  // transform into freq. domain
> N=smp_array.size();
> for i=0..N/2 {
>   if((i*N) <= f1 || (i*N) >= f2) // if freq. is outside passband...
>       smp_array[i]=0; smp_array[N-i]=0;  // set power to 0
> }
> transform(); // transform back as DHT is its own inverse (right?)
> }
>
>
> Ok, I bet I totally messed up the filtering routine. From the docs I would
> assume that each index represents th fraction of freq./N. So I thought I
> have to multiply it to recover the real frequencie in Hz. Is that a false
> assumption? If yes, how is it done correctly?
>
> Thanks in advance
> Alex
>
>
>
>
> This message was sent using the Comp.DSP web interface on
> www.DSPRelated.com

Hi!

I am sampling at 256Hz and I want to filter out all frequencies except
1-9Hz. The filter should have the highest attenuation as possible. My
approach would be to use the FFTW library, transform from time to freq.
domain and set all bins which are outside of my passband to 0 (like a
rectangular window).
But honestly I dont have a clue how to accomplish that with FFTW. I tried
a hartley transform, applied my window function and transformed it back,
but I always get a "half-period-sine-like-shaped" signal. So I did
obviously something wrong.

I used following code (sketched):

func transform() {
// Im using the same array here, this is no typo!
// could this be a problem?
fftw_create_plan_r2r_1d(sampleNum, smp_array, smp_array, FFTW_DHT,
FFTW_ESTIMATE);
execute...
}

func filter(f1, f2) {
transform();  // transform into freq. domain
N=smp_array.size();
for i=0..N/2 {
  if((i*N) <= f1 || (i*N) >= f2) // if freq. is outside passband...
      smp_array[i]=0; smp_array[N-i]=0;  // set power to 0
}
transform(); // transform back as DHT is its own inverse (right?)
}


Ok, I bet I totally messed up the filtering routine. From the docs I would
assume that each index represents th fraction of freq./N. So I thought I
have to multiply it to recover the real frequencie in Hz. Is that a false
assumption? If yes, how is it done correctly?

Thanks in advance
Alex



		
This message was sent using the Comp.DSP web interface on
www.DSPRelated.com