Re: How to perform the 1-D signal filter? [message #58425] |
Fri, 01 February 2008 06:58  |
Kenneth P. Bowman
Messages: 585 Registered: May 2000
|
Senior Member |
|
|
In article
<e9b5822c-2240-4937-ad85-a53f057d9ec0@p69g2000hsa.googlegroups.com>,
"duxiyu@gmail.com" <duxiyu@gmail.com> wrote:
> Dear all,
>
> Here I give a signal example and hope somone can show me how to
> perferm the frequency filter on it.
>
> ;creat a signal data with two peaks in frequency domain at 2 and 3 Hz.
> t=findgen(1000)/10.
> data=sin(2*!pi*2*t)+sin(2*!pi*3*t)
>
> freq=findgen(501)/100.
> v=fft(data)
> plot,freq,abs(v[0:500])^2,xtitle='frequency',ytitle='spectru m'
>
>
> I want to filter the signal with the frequency higher than 2.5 Hz. How
> do I do this?
>
> I have read the help files about Digital_Filter and Convol, but I do
> not know how to select the parameters for Signal_Filter.
>
> Du
You could look at the chapter on FFTs and digital filtering in my IDL
book (http://idl.tamu.edu). I wrote the chapter in part so that *I*
could refer to it whenever I need to do a digital filter. :-)
Ken Bowman
|
|
|
Re: How to perform the 1-D signal filter? [message #58427 is a reply to message #58425] |
Fri, 01 February 2008 09:16   |
Wox
Messages: 184 Registered: August 2006
|
Senior Member |
|
|
On Fri, 1 Feb 2008 07:16:26 -0700, David Fanning <news@dfanning.com>
wrote:
> Wox writes:
>
>> Example below filters in time or frequency domain:
>>
>>
>> ; Time domain
>> freq1=2.
>> freq2=3.
>> freq3=4.
>> dtime=0.1
>> ntime=1000
>>
>> time=dtime*findgen(ntime)
>> signal=sin(2*!pi*freq1*time)+sin(2*!pi*freq2*time)+sin(2*!pi *freq3*time)
>>
>> ; Time domain Filter
>> f_low = 0
>> f_high = 2.5
>> timefilter = DIGITAL_FILTER(f_low*2*dtime, f_high*2*dtime, 50.,40)
>> signal=convol(signal,timefilter)
>>
>> ; Frequency domain
>> nfreq=ntime/2+1
>>
>> freq=findgen(nfreq)/(dtime*ntime)
>> fsignal=fft(signal)
>>
>> ; Frequency domain filter (instead of time domain filter)
>> if n_elements(timefilter) eq 0 then begin
>> steep=20.
>> freqfilter= 1./(1.+(freq/f_high)^steep)
>> fsignal*=freqfilter
>> endif
>>
>> plot,freq,abs(fsignal[0:nfreq-1])^2,xtitle='frequency',ytitl e='spectrum'
>
> Wonderful example, but I'm trying to understand this whole
> subject. Do you think you could flush this out with a little
> explanation of what you are doing and why you choose the terms
> you use, etc.? What kind of frequency filter are you constructing
> here? I don't necessarily see it doing any filtering of the signal,
> at least if I pass it the original signal, rather than the signal
> that had already been filtered in the time domain, as written
> in your example.
>
> Cheers,
>
> Confused
Ok, sorry for the confusion, but I was just illustrating that you can
do the same filtering in the frequency domain as in the time domain.
You do one or the other, not both at the same time. Btw, convolution
in one domain becomes multiplication in the other:
filtered = signal "convol" filter
fft(filtered) = fft(signal) x fft(filter)
But I guess you already knew all this.
The filter used is the Kaiser-Bessel filter. At least I think
digital_filter is using this. For the filter I constructed in the
fourier domain, I'm not quit sure whether it is really identical to
the KB filter, but if you plot it, it looks like a nice lowpass filter
to me :-).
I'm just typing this in a hurry... Did I answer your questions?
|
|
|
Re: How to perform the 1-D signal filter? [message #58429 is a reply to message #58425] |
Fri, 01 February 2008 08:57   |
Wox
Messages: 184 Registered: August 2006
|
Senior Member |
|
|
On Fri, 1 Feb 2008 06:54:42 -0800 (PST), "duxiyu@gmail.com"
<duxiyu@gmail.com> wrote:
> In the example, you set A=50 and Nterms=40 and I do not know where
> their values come from and how I get them.
The Nterms=40 determines the width of the filter: width=2.Nterms+1
The A=50 has something to do with the Gibbs phenomena, i.e. the
wiggles you get when you Fourier expand a function with a finite
number of terms (as opposed to infinite). The wiggles are visible at
discontinuities.
I don't know anything about the Kaiser-Bessel filter (I think that's
the one used here). So maybe someone else can explain this.
> Seondly, I want to know the relation between the filter in time and
> frequency domain.
> if the following command is excuted, I think newfsignal is exactly
> equal to fft(newsignal), isn't it?
> "newsignal=convol(signal,timefilter)
> fsignal=fft(signal)
> newfsignal=fsignal*freqfilter"
You are absolutly right. That's why in the code comments I wrote: ";
Frequency domain filter (instead of time domain filter)". This means
you can choose in which domain you are filtering, the time or
frequency(fourier) domain. You could even have a filter in the time
domain, take the fourier transform and use that to apply the filter in
the frequency domain. Keep in mind however that every fft introduces
some errors.
> And I also do not konw why you set steep=20.
If you plot the lowpass filter (in the fourier domain):
steep=20.
freqfilter= 1./(1.+(freq/f_high)^steep)
plot,freq,freqfilter
Do this for different "steep" and you will see that the edge becomes
sharper with higher "steep".
|
|
|
Re: How to perform the 1-D signal filter? [message #58434 is a reply to message #58425] |
Fri, 01 February 2008 06:54   |
duxiyu@gmail.com
Messages: 88 Registered: March 2007
|
Member |
|
|
Thank you very much.
But I still have some questions about the parameters in the procedure
Digital_Filter.
Result = DIGITAL_FILTER( Flow, Fhigh, A, Nterms [, /DOUBLE] )
In the example, you set A=50 and Nterms=40 and I do not know where
their values come from and how I get them.
Seondly, I want to know the relation between the filter in time and
frequency domain.
if the following command is excuted, I think newfsignal is exactly
equal to fft(newsignal), isn't it?
"newsignal=convol(signal,timefilter)
fsignal=fft(signal)
newfsignal=fsignal*freqfilter"
And I also do not konw why you set steep=20.
Du
On Feb 1, 9:27 pm, Wox <nom...@hotmail.com> wrote:
> On Fri, 1 Feb 2008 02:20:52 -0800 (PST), "dux...@gmail.com"
>
>
>
> <dux...@gmail.com> wrote:
>> Dear all,
>
>> Here I give a signal example and hope somone can show me how to
>> perferm the frequency filter on it.
>
>> ;creat a signal data with two peaks in frequency domain at 2 and 3 Hz.
>> t=findgen(1000)/10.
>> data=sin(2*!pi*2*t)+sin(2*!pi*3*t)
>
>> freq=findgen(501)/100.
>> v=fft(data)
>> plot,freq,abs(v[0:500])^2,xtitle='frequency',ytitle='spectru m'
>
>> I want to filter the signal with the frequency higher than 2.5 Hz. How
>> do I do this?
>
>> I have read the help files about Digital_Filter and Convol, but I do
>> not know how to select the parameters for Signal_Filter.
>
>> Du
>
> Example below filters in time or frequency domain:
>
> ; Time domain
> freq1=2.
> freq2=3.
> freq3=4.
> dtime=0.1
> ntime=1000
>
> time=dtime*findgen(ntime)
> signal=sin(2*!pi*freq1*time)+sin(2*!pi*freq2*time)+sin(2*!pi *freq3*time)
>
> ; Time domain Filter
> f_low = 0
> f_high = 2.5
> timefilter = DIGITAL_FILTER(f_low*2*dtime, f_high*2*dtime, 50.,40)
> signal=convol(signal,timefilter)
>
> ; Frequency domain
> nfreq=ntime/2+1
>
> freq=findgen(nfreq)/(dtime*ntime)
> fsignal=fft(signal)
>
> ; Frequency domain filter (instead of time domain filter)
> if n_elements(timefilter) eq 0 then begin
> steep=20.
> freqfilter= 1./(1.+(freq/f_high)^steep)
> fsignal*=freqfilter
> endif
>
> plot,freq,abs(fsignal[0:nfreq-1])^2,xtitle='frequency',ytitl e='spectrum'
|
|
|
Re: How to perform the 1-D signal filter? [message #58436 is a reply to message #58434] |
Fri, 01 February 2008 06:16   |
David Fanning
Messages: 11724 Registered: August 2001
|
Senior Member |
|
|
Wox writes:
> Example below filters in time or frequency domain:
>
>
> ; Time domain
> freq1=2.
> freq2=3.
> freq3=4.
> dtime=0.1
> ntime=1000
>
> time=dtime*findgen(ntime)
> signal=sin(2*!pi*freq1*time)+sin(2*!pi*freq2*time)+sin(2*!pi *freq3*time)
>
> ; Time domain Filter
> f_low = 0
> f_high = 2.5
> timefilter = DIGITAL_FILTER(f_low*2*dtime, f_high*2*dtime, 50.,40)
> signal=convol(signal,timefilter)
>
> ; Frequency domain
> nfreq=ntime/2+1
>
> freq=findgen(nfreq)/(dtime*ntime)
> fsignal=fft(signal)
>
> ; Frequency domain filter (instead of time domain filter)
> if n_elements(timefilter) eq 0 then begin
> steep=20.
> freqfilter= 1./(1.+(freq/f_high)^steep)
> fsignal*=freqfilter
> endif
>
> plot,freq,abs(fsignal[0:nfreq-1])^2,xtitle='frequency',ytitl e='spectrum'
Wonderful example, but I'm trying to understand this whole
subject. Do you think you could flush this out with a little
explanation of what you are doing and why you choose the terms
you use, etc.? What kind of frequency filter are you constructing
here? I don't necessarily see it doing any filtering of the signal,
at least if I pass it the original signal, rather than the signal
that had already been filtered in the time domain, as written
in your example.
Cheers,
Confused
--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.dfanning.com/
Sepore ma de ni thui. ("Perhaps thou speakest truth.")
|
|
|
Re: How to perform the 1-D signal filter? [message #58438 is a reply to message #58436] |
Fri, 01 February 2008 05:27   |
Wox
Messages: 184 Registered: August 2006
|
Senior Member |
|
|
On Fri, 1 Feb 2008 02:20:52 -0800 (PST), "duxiyu@gmail.com"
<duxiyu@gmail.com> wrote:
> Dear all,
>
> Here I give a signal example and hope somone can show me how to
> perferm the frequency filter on it.
>
> ;creat a signal data with two peaks in frequency domain at 2 and 3 Hz.
> t=findgen(1000)/10.
> data=sin(2*!pi*2*t)+sin(2*!pi*3*t)
>
> freq=findgen(501)/100.
> v=fft(data)
> plot,freq,abs(v[0:500])^2,xtitle='frequency',ytitle='spectru m'
>
>
> I want to filter the signal with the frequency higher than 2.5 Hz. How
> do I do this?
>
> I have read the help files about Digital_Filter and Convol, but I do
> not know how to select the parameters for Signal_Filter.
>
> Du
>
>
Example below filters in time or frequency domain:
; Time domain
freq1=2.
freq2=3.
freq3=4.
dtime=0.1
ntime=1000
time=dtime*findgen(ntime)
signal=sin(2*!pi*freq1*time)+sin(2*!pi*freq2*time)+sin(2*!pi *freq3*time)
; Time domain Filter
f_low = 0
f_high = 2.5
timefilter = DIGITAL_FILTER(f_low*2*dtime, f_high*2*dtime, 50.,40)
signal=convol(signal,timefilter)
; Frequency domain
nfreq=ntime/2+1
freq=findgen(nfreq)/(dtime*ntime)
fsignal=fft(signal)
; Frequency domain filter (instead of time domain filter)
if n_elements(timefilter) eq 0 then begin
steep=20.
freqfilter= 1./(1.+(freq/f_high)^steep)
fsignal*=freqfilter
endif
plot,freq,abs(fsignal[0:nfreq-1])^2,xtitle='frequency',ytitl e='spectrum'
|
|
|
Re: How to perform the 1-D signal filter? [message #58548 is a reply to message #58425] |
Sat, 02 February 2008 07:20  |
duxiyu@gmail.com
Messages: 88 Registered: March 2007
|
Member |
|
|
On 2月1日, 下午10时58分, "Kenneth P. Bowman" <k-bow...@null.edu> wrote:
> In article
> < e9b5822c-2240-4937-ad85-a53f057d9...@p69g2000hsa.googlegroup s.com >,
>
>
>
> "dux...@gmail.com" <dux...@gmail.com> wrote:
>> Dear all,
>
>> Here I give a signal example and hope somone can show me how to
>> perferm the frequency filter on it.
>
>> ;creat a signal data with two peaks in frequency domain at 2 and 3 Hz.
>> t=findgen(1000)/10.
>> data=sin(2*!pi*2*t)+sin(2*!pi*3*t)
>
>> freq=findgen(501)/100.
>> v=fft(data)
>> plot,freq,abs(v[0:500])^2,xtitle='frequency',ytitle='spectru m'
>
>> I want to filter the signal with the frequency higher than 2.5 Hz. How
>> do I do this?
>
>> I have read the help files about Digital_Filter and Convol, but I do
>> not know how to select the parameters for Signal_Filter.
>
>> Du
>
> You could look at the chapter on FFTs and digital filtering in my IDL
> book (http://idl.tamu.edu). I wrote the chapter in part so that *I*
> could refer to it whenever I need to do a digital filter. :-)
>
> Ken Bowman
It is a pity that I cannot buy your book in my country.
Du
|
|
|
Re: How to perform the 1-D signal filter? [message #58549 is a reply to message #58427] |
Sat, 02 February 2008 07:14  |
duxiyu@gmail.com
Messages: 88 Registered: March 2007
|
Member |
|
|
On Feb 2, 1:16 am, Wox <nom...@hotmail.com> wrote:
> On Fri, 1 Feb 2008 07:16:26 -0700, David Fanning <n...@dfanning.com>
> wrote:
>
>
>
>> Wox writes:
>
>>> Example below filters in time or frequency domain:
>
>>> ; Time domain
>>> freq1=2.
>>> freq2=3.
>>> freq3=4.
>>> dtime=0.1
>>> ntime=1000
>
>>> time=dtime*findgen(ntime)
>>> signal=sin(2*!pi*freq1*time)+sin(2*!pi*freq2*time)+sin(2*!pi *freq3*time)
>
>>> ; Time domain Filter
>>> f_low = 0
>>> f_high = 2.5
>>> timefilter = DIGITAL_FILTER(f_low*2*dtime, f_high*2*dtime, 50.,40)
>>> signal=convol(signal,timefilter)
>
>>> ; Frequency domain
>>> nfreq=ntime/2+1
>
>>> freq=findgen(nfreq)/(dtime*ntime)
>>> fsignal=fft(signal)
>
>>> ; Frequency domain filter (instead of time domain filter)
>>> if n_elements(timefilter) eq 0 then begin
>>> steep=20.
>>> freqfilter= 1./(1.+(freq/f_high)^steep)
>>> fsignal*=freqfilter
>>> endif
>
>>> plot,freq,abs(fsignal[0:nfreq-1])^2,xtitle='frequency',ytitl e='spectrum'
>
>> Wonderful example, but I'm trying to understand this whole
>> subject. Do you think you could flush this out with a little
>> explanation of what you are doing and why you choose the terms
>> you use, etc.? What kind of frequency filter are you constructing
>> here? I don't necessarily see it doing any filtering of the signal,
>> at least if I pass it the original signal, rather than the signal
>> that had already been filtered in the time domain, as written
>> in your example.
>
>> Cheers,
>
>> Confused
>
> Ok, sorry for the confusion, but I was just illustrating that you can
> do the same filtering in the frequency domain as in the time domain.
> You do one or the other, not both at the same time. Btw, convolution
> in one domain becomes multiplication in the other:
>
> filtered = signal "convol" filter
> fft(filtered) = fft(signal) x fft(filter)
>
> But I guess you already knew all this.
>
> The filter used is the Kaiser-Bessel filter. At least I think
> digital_filter is using this. For the filter I constructed in the
> fourier domain, I'm not quit sure whether it is really identical to
> the KB filter, but if you plot it, it looks like a nice lowpass filter
> to me :-).
>
> I'm just typing this in a hurry... Did I answer your questions?
Thank you again for your detailed explaination.
Now I am clear about it.
Best regards,
Du
|
|
|
Re: How to perform the 1-D signal filter? [message #58551 is a reply to message #58425] |
Fri, 01 February 2008 17:07  |
David Fanning
Messages: 11724 Registered: August 2001
|
Senior Member |
|
|
Kenneth P. Bowman writes:
> You could look at the chapter on FFTs and digital filtering in my IDL
> book (http://idl.tamu.edu). I wrote the chapter in part so that *I*
> could refer to it whenever I need to do a digital filter. :-)
Humm. I found quite a lot of interest in the last three
or four chapters in that book, and had a nice little read
this afternoon. Unlike Chloe, I have no need to impress
anyone anymore, so I'm just going to steal one or two
of those ideas, if you don't mind. :-)
Thanks,
David
--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.dfanning.com/
Sepore ma de ni thui. ("Perhaps thou speakest truth.")
|
|
|
Re: How to perform the 1-D signal filter? [message #58566 is a reply to message #58425] |
Fri, 01 February 2008 10:30  |
David Fanning
Messages: 11724 Registered: August 2001
|
Senior Member |
|
|
Kenneth P. Bowman writes:
> You could look at the chapter on FFTs and digital filtering in my IDL
> book (http://idl.tamu.edu). I wrote the chapter in part so that *I*
> could refer to it whenever I need to do a digital filter. :-)
Ah, right. I remember that chapter. Here you are going
along in a nice, beginning book and all of a sudden you
are surrounded by formidable equations! All I really
remember about it was being dazed and thinking, "Boy,
that was a big ending."
I'll have another look. :-)
Cheers,
David
--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.dfanning.com/
Sepore ma de ni thui. ("Perhaps thou speakest truth.")
|
|
|