R help - Nov 2009 - fft and filtering puzzle

I am puzzled by a filtering problem using fft(). I don't blame R.

I have a waveform y consisting of the sum of 2 sinewaves having freqs f1 and f2.
I do s = fft() of y.
Remove s's spike at freq=f2
Do inverse fft on s.
The resulting waveform still has a lot of f2 in it! But the filtering
should have removed it all.
What is going on, and how to fix??

Thanks very much for any help.
Bill


Below is code illustrating the problem.

n<-1024
#f1=2, f2=7
y<- 9*cos(2*pi*2.0*(1:n)/n) + 3*cos(2*pi*7.0*(1:n)/n)
plot(y,type='l')

s<-2*Re(fft(y))/n
plot(s[1:20], type='l')

#note freq[i] is s[i+1]; e.g. f=7.0 peak is at s[8]
coef<- 1-(((1:n)>8-1) & ((1:n)<8+1))
s2<- fft(y)*coef   # I have now removed f2

plot(2*Re(s2[1:20])/n, type='l')   #f2 is removed from spectrum

y2<-Re(fft(s2,inverse=TRUE)/n)  #now get back to time domain

plot(y,type='l')  #original waveform f1+f2
lines(y2,type='l',col='blue' )  #filtered waveform
lines(9*cos(2*pi*2.0*(1:n)/n), type='l',col='red' )  #what it
should
look like -- f2 removed

I think I've got it!

I only filtered out the pos freq and not the negative! (At the far end
of s). I will check now.

It would be nice to forget about that bit s[n/2:n]. But I get if I get
rid of it all, then when I back transform into time domain I imagine
I'll get artifacts.

Bill

On Sat, Nov 28, 2009 at 9:53 AM, William Simpson
<william.a.simpson at gmail.com> wrote:
> I am puzzled by a filtering problem using fft(). I don't blame R.
>
> I have a waveform y consisting of the sum of 2 sinewaves having freqs f1
and f2.
> I do s = fft() of y.
> Remove s's spike at freq=f2
> Do inverse fft on s.
> The resulting waveform still has a lot of f2 in it! But the filtering
> should have removed it all.
> What is going on, and how to fix??
>
> Thanks very much for any help.
> Bill
>
>
> Below is code illustrating the problem.
>
> n<-1024
> #f1=2, f2=7
> y<- 9*cos(2*pi*2.0*(1:n)/n) + 3*cos(2*pi*7.0*(1:n)/n)
> plot(y,type='l')
>
> s<-2*Re(fft(y))/n
> plot(s[1:20], type='l')
>
> #note freq[i] is s[i+1]; e.g. f=7.0 peak is at s[8]
> coef<- 1-(((1:n)>8-1) & ((1:n)<8+1))
> s2<- fft(y)*coef ? # I have now removed f2
>
> plot(2*Re(s2[1:20])/n, type='l') ? #f2 is removed from spectrum
>
> y2<-Re(fft(s2,inverse=TRUE)/n) ?#now get back to time domain
>
> plot(y,type='l') ?#original waveform f1+f2
> lines(y2,type='l',col='blue' ) ?#filtered waveform
> lines(9*cos(2*pi*2.0*(1:n)/n), type='l',col='red' ) ?#what
it should
> look like -- f2 removed
>

Yep, that was it. Needed to remove neg freq from opposite end of s[]

coef<- 1-(((1:n)>8-1) & ((1:n)<8+1) |  (((1:n)>(n-6)-1) &
((1:n)<(n-6)+1))     )

Sorry to bother you.

Cheers
Bill