Hi, I am using FFTW to try another test, f(t)=exp(-x), x>=0, the real part of the Fourier transforms should be 1/(1+x^2). I use xmax=1000, and xmin=0, let N=8192. I found the results seems shifted upper, every output value is larger than the analytical solution. My code is like this: program test implicit none include "fftw3.f" integer N parameter(N=8192) integer pi parameter(pi=3.14159265357) double complex in double complex out double precision f dimension out(N) dimension in(N) dimension f(N) integer i real*8 :: x,y,xmax,xmin,dx integer*8 plan OPEN(UNIT=7,FILE='fftin') OPEN(UNIT=8,FILE='fftout') ! write(*,*) 'Input array:' xmax=1000. xmin=0. dx=(xmax-xmin)/N write(6,*) 'XMAX',xmax,'xmin',xmin,'dx',dx,'N',N,'omega',1.0/(2.0*dx)*2*pi do i = 1,N/2+1 x=xmin+dx*(i-1) in(i)=exp(-x)*(-1)**(i+1) if(i.gt.N/2+1) in(i)=0.0 enddo do i=1,N x=xmin+dx*(i-1) write(7,*) i,x,real(in(i)), imag(in(i)) enddo close(7) !call dfftw_plan_r2r_1d (plan, N, in, out, fftw_r2hc, FFTW_ESTIMATE) !call dfftw_plan_dft_r2c_1d(plan,N,in,out,FFTW_ESTIMATE) call dfftw_plan_dft_1d (plan, N, in, out, FFTW_FORWARD, FFTW_ESTIMATE) call dfftw_execute(plan) !write(*,*) 'Output array:' do i = 1,N f(i)=(-N/2+i-1)/(N*dx) write(8,*) i,f(i)*2*pi,real(out(i))*dx,imag(out(i))*dx,1./(1.+(f(i)*2.*pi)**2),-f(i)*2.*pi/(1.+(f(i)*2.*pi)**2) enddo close(8) STOP can any one help me about it? Thanks
another FFTW problem
Started by ●July 3, 2008
Reply by ●July 3, 20082008-07-03
jv_ju wrote:> Hi, I am using FFTW to try another test, f(t)=exp(-x), x>=0, the real part > of the Fourier transforms should be 1/(1+x^2). I use xmax=1000, and xmin=0, > let N=8192. I found the results seems shifted upper, every output value is > larger than the analytical solution. > My code is like this: > program test > > implicit none > > include "fftw3.f" > > integer N > parameter(N=8192) > integer pi > parameter(pi=3.14159265357) > > double complex in > double complex out > double precision f > dimension out(N) > dimension in(N) > dimension f(N) > > integer i > real*8 :: x,y,xmax,xmin,dx > integer*8 plan > > OPEN(UNIT=7,FILE='fftin') > OPEN(UNIT=8,FILE='fftout') > > ! write(*,*) 'Input array:' > > xmax=1000. > xmin=0. > dx=(xmax-xmin)/N > > write(6,*) > 'XMAX',xmax,'xmin',xmin,'dx',dx,'N',N,'omega',1.0/(2.0*dx)*2*pi > > > do i = 1,N/2+1 > x=xmin+dx*(i-1) > in(i)=exp(-x)*(-1)**(i+1) > if(i.gt.N/2+1) in(i)=0.0 > enddo > do i=1,N > x=xmin+dx*(i-1) > write(7,*) i,x,real(in(i)), imag(in(i)) > enddo > close(7) > > > !call dfftw_plan_r2r_1d (plan, N, in, out, fftw_r2hc, FFTW_ESTIMATE) > !call dfftw_plan_dft_r2c_1d(plan,N,in,out,FFTW_ESTIMATE) > call dfftw_plan_dft_1d (plan, N, in, out, FFTW_FORWARD, FFTW_ESTIMATE) > > call dfftw_execute(plan) > > !write(*,*) 'Output array:' > > do i = 1,N > f(i)=(-N/2+i-1)/(N*dx) > write(8,*) > i,f(i)*2*pi,real(out(i))*dx,imag(out(i))*dx,1./(1.+(f(i)*2.*pi)**2),-f(i)*2.*pi/(1.+(f(i)*2.*pi)**2) > enddo > close(8) > STOP > > can any one help me about it?An FFT doesn't compute a continuous transform of a continuous function. Because what is being transformed is necessarily discrete in time, the result of the transform is necessarily periodic and discrete in frequency. 1/(1+x^2) meets neither of those constraints. Your tool can't produce it. Jerry -- Engineering is the art of making what you want from things you can get. �����������������������������������������������������������������������
Reply by ●July 10, 20082008-07-10
> > >An FFT doesn't compute a continuous transform of a continuous function. >Because what is being transformed is necessarily discrete in time, the >result of the transform is necessarily periodic and discrete in >frequency. 1/(1+x^2) meets neither of those constraints. Your tool can't>produce it. > >Jerry >-- >Engineering is the art of making what you want from things you can get. >����������������������������������������������������������������������� >How about you consider it as a periodic function with the period of 1000?
Reply by ●July 10, 20082008-07-10
> > >An FFT doesn't compute a continuous transform of a continuous function. >Because what is being transformed is necessarily discrete in time, the >result of the transform is necessarily periodic and discrete in >frequency. 1/(1+x^2) meets neither of those constraints. Your tool can't>produce it. > >Jerry >-- >Engineering is the art of making what you want from things you can get. >����������������������������������������������������������������������� >How about you consider it as a periodic function with the period of 1000?
Reply by ●July 12, 20082008-07-12
jv_ju wrote:>> >> An FFT doesn't compute a continuous transform of a continuous function. >> Because what is being transformed is necessarily discrete in time, the >> result of the transform is necessarily periodic and discrete in >> frequency. 1/(1+x^2) meets neither of those constraints. Your tool can't>> produce it....> How about you consider it as a periodic function with the period of 1000?Will wishing make it so? Jerry -- Engineering is the art of making what you want from things you can get. �����������������������������������������������������������������������