| How to perform the 1-D signal filter? [message #58441] |
Fri, 01 February 2008 02:20  |
duxiyu@gmail.com
Messages: 88
Registered: March 2007
|
Member |
|
|
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
|
|
|
|
| Re: How to perform the 1-D signal filter? [message #58529 is a reply to message #58441] |
Mon, 04 February 2008 08:48  |
jdu
Messages: 7
Registered: February 2008
|
Junior Member |
|
|
On 2月4日, 下午5时28分, Wox <nom...@hotmail.com> wrote:
> On Sat, 2 Feb 2008 11:27:47 -0800 (PST), jdu.u...@gmail.com wrote:
>> I found that the the value of the peak in the second figure is lower
>> than those in first and third figure. Why does it happen?
>> And the filtered signal_1 is not equal to the real part of signal_2.
>> Does it mean that the filter in time domain is not exactly equal to it
>> in frequency domain?
>
> Well, I'm not sure. If you execute the code below, you will see that
> the two filters are similar. Maybe this is just the result of the
> intrinsic approximation of time domain filtering?
>
> ; 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)
>
> ; Frequency domain
> nfreq=ntime/2+1
> freq=findgen(nfreq)/(dtime*ntime)
> fsignal=fft(signal)
>
> ; 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)
>
> ; Time domain Filter in Frequency domain
> ntime2=n_elements(timefilter)
> nfreq2=ntime2/2+1
> freq2=findgen(nfreq2)/(dtime*ntime2)
> ftimefilter=fft(timefilter,1); == ntime2*fft(timefilter)
>
> ; Frequency domain filter (instead of time domain filter)
> steep=50.
> freqfilter= 1./(1.+(freq/f_high)^steep)
> fsignal*=freqfilter
>
> ; Plot
> window,0
> plot,freq,abs(fsignal[0:nfreq-1])^2,xtitle='frequency',ytitl e='spectrum',title='Time
> domain filtered'
> window,1
> plot,freq,abs((fft(signal))[0:nfreq-1])^2,xtitle='frequency' ,ytitle='spectrum',title='Freq
> domain filtered'
> window,2
> plot,freq,abs(freqfilter)^2,yrange=[0,1.1],ystyle=1,xtitle=' frequency',ytitle='spectrum',title='Filters'
> oplot,freq2,abs((ftimefilter)[0:nfreq2-1])^2,col=128
In this example, max(abs(fsignal[0:nfreq-1])^2) = 0.249995 and
max(abs((fft(signal))[0:nfreq-1])^2) = 0.211242.
The peak of spectra in window 0 is not equal to it in window 1.
Does this difference between them result from the arithmetical errors
of the filters?
Is there any methods to eliminate it?
> ftimefilter=fft(timefilter,1); == ntime2*fft(timefilter)
Besides, I don't know why you use the inverse transform instead of
the forward transform in IDL to caculate ftimefilter.
Du
|
|
|
|
| Re: How to perform the 1-D signal filter? [message #58538 is a reply to message #58441] |
Mon, 04 February 2008 01:28 |
Wox
Messages: 184
Registered: August 2006
|
Senior Member |
|
|
On Sat, 2 Feb 2008 11:27:47 -0800 (PST), jdu.ustc@gmail.com wrote:
> I found that the the value of the peak in the second figure is lower
> than those in first and third figure. Why does it happen?
> And the filtered signal_1 is not equal to the real part of signal_2.
> Does it mean that the filter in time domain is not exactly equal to it
> in frequency domain?
Well, I'm not sure. If you execute the code below, you will see that
the two filters are similar. Maybe this is just the result of the
intrinsic approximation of time domain filtering?
; 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)
; Frequency domain
nfreq=ntime/2+1
freq=findgen(nfreq)/(dtime*ntime)
fsignal=fft(signal)
; 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)
; Time domain Filter in Frequency domain
ntime2=n_elements(timefilter)
nfreq2=ntime2/2+1
freq2=findgen(nfreq2)/(dtime*ntime2)
ftimefilter=fft(timefilter,1); == ntime2*fft(timefilter)
; Frequency domain filter (instead of time domain filter)
steep=50.
freqfilter= 1./(1.+(freq/f_high)^steep)
fsignal*=freqfilter
; Plot
window,0
plot,freq,abs(fsignal[0:nfreq-1])^2,xtitle='frequency',ytitl e='spectrum',title='Time
domain filtered'
window,1
plot,freq,abs((fft(signal))[0:nfreq-1])^2,xtitle='frequency' ,ytitle='spectrum',title='Freq
domain filtered'
window,2
plot,freq,abs(freqfilter)^2,yrange=[0,1.1],ystyle=1,xtitle=' frequency',ytitle='spectrum',title='Filters'
oplot,freq2,abs((ftimefilter)[0:nfreq2-1])^2,col=128
|
|
|
|
| Re: How to perform the 1-D signal filter? [message #58544 is a reply to message #58441] |
Sat, 02 February 2008 11:27 |
jdu
Messages: 7
Registered: February 2008
|
Junior Member |
|
|
There is another question about the filter example.
The code of example is below:
; 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_1=convol(signal,timefilter)
fsignal_1=fft(signal_1)
; Frequency domain
nfreq=ntime/2+1
freq=findgen(nfreq)/(dtime*ntime)
fsignal=fft(signal)
; Frequency domain filter (instead of time domain filter)
steep=100.d
y=1./(1.+(freq/f_high)^steep)
freqfilter=[y[0:nfreq-2],reverse(y[0:nfreq-2])]
fsignal_2=fsignal*freqfilter
signal_2=fft(fsignal_2,1)
window,/free
plot,freq,abs(fsignal[0:nfreq-1])^2,title='Orignal
Data',xtitle='frequency',ytitle='spectrum'
window,/free
plot,freq,abs(fsignal_1[0:nfreq-1])^2,title='Filter in Time
domain',xtitle='frequency',ytitle='spectrum'
window,/free
plot,freq,abs(fsignal_2[0:nfreq-1])^2,title='Filter in Frequency
domain',xtitle='frequency',ytitle='spectrum'
I found that the the value of the peak in the second figure is lower
than those in first and third figure. Why does it happen?
And the filtered signal_1 is not equal to the real part of signal_2.
Does it mean that the filter in time domain is not exactly equal to it
in frequency domain?
Du
|
|
|
|
| Re: How to perform the 1-D signal filter? [message #58550 is a reply to message #58441] |
Fri, 01 February 2008 15:27 |
Kenneth P. Bowman
Messages: 585
Registered: May 2000
|
Senior Member |
|
|
In article <MPG.220d13322fc9fa3c98a25f@news.frii.com>,
David Fanning <news@dfanning.com> wrote:
> 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
The math is there is you want it, but it is not essential, I suppose, to
using a digital filter. The programming examples are as straightforward
as is possible when using FFTs. ;-)
Ken
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
