Prev: Suggestion for DS-CDMA simulation issue (just have a look)
Next: Practical differences between BCH and Reed-Solomon codes?
From: JochenH on 9 Nov 2009 06:03
I'm not entirely shure if this is the right place to ask, so if my post
does not fit please ignore or suggest a better place.

I want to do a FFT on real data using the FFTW3 library. The code i'm
working on is based on hermitian matrices, therefore i need to use the
r2c transform.
As far as i understand it, the first 2n+1 entries with a c2c transform
on real data should be the same as with a r2c transform.
However, i can't reproduce the c2c output with a r2c transform.
Thanks.


Fortran source code:

program kint

use fftw3

implicit none

complex(kind=8), allocatable, dimension(:) :: fr,fk,fk2

real(kind=8) :: dx,dk,a,pi
complex(kind=8) :: resultk
integer :: nx,i,plan

nx = 36
dx = 1.d0
pi = 3.141592653589793d0
dk = 2.d0*pi/(dble(nx)*dx)
a = 0.2d0

allocate(fr(nx),fk(nx))

do i=1,nx
fr(i) = cmplx(exp(-a * (dx*(dble(i-1)-dble(nx)/2.d0)-2)**2),0.d0)
end do

!**************** c2c *************

call dfftw_plan_dft_1d(plan,nx,fr,fk,FFTW_BACKWARD,FFTW_ESTIMATE)
call dfftw_execute(plan)
call dfftw_destroy_plan(plan)

do i=1,nx
write(*,'(i,f,f)') i, dreal(fk(i)),dimag(fk(i))
end do

fk = fk * (dx / sqrt(2.d0 * pi))

deallocate(fk)

allocate(fk2(nx/2+1))

write(*,'(a)') ''

!**************** c2r *************

call dfftw_plan_dft_r2c_1d(plan,nx,fr,fk2,FFTW_BACKWARD,FFTW_ESTIMATE)
call dfftw_execute(plan)
call dfftw_destroy_plan(plan)

do i=1,nx/2+1
write(*,'(i,f,f)') i, dreal(fk2(i)),dimag(fk2(i))
end do

deallocate(fr,fk2)

end program

Output:

1 3.9633273373814508 0.0000000000000000
2 -3.5851642010512617 -1.3048930541394301
3 2.6071563537030196 2.1876639348137057
4 -1.4066938032409493 -2.4364651379056217
5 0.3742336568932738 2.1223845346278751
6 0.2656516214009662 -1.5065852108032667
7 -0.5031627214892469 0.8715033980940043
8 0.4699105738175311 -0.3943017891765944
9 -0.3256090393846692 0.1185119983440311
10 0.1813788371283939 0.0000000000000000
11 -0.0826751513495624 -0.0300912942046938
12 0.0302950931975201 0.0254206015285206
13 -0.0082365236849697 -0.0142660775001120
14 0.0011041449519306 0.0062619171921073
15 0.0003949462596618 -0.0022398515421467
16 -0.0003770206872807 0.0006530189858746
17 0.0001780852757594 -0.0001494312892079
18 -0.0000658764105576 0.0000239770525833
19 0.0000347119594307 -0.0000000000000002
20 -0.0000658764105575 -0.0000239770525833
21 0.0001780852757595 0.0001494312892079
22 -0.0003770206872806 -0.0006530189858747
23 0.0003949462596617 0.0022398515421467
24 0.0011041449519307 -0.0062619171921072
25 -0.0082365236849697 0.0142660775001120
26 0.0302950931975202 -0.0254206015285206
27 -0.0826751513495624 0.0300912942046940
28 0.1813788371283939 0.0000000000000000
29 -0.3256090393846690 -0.1185119983440312
30 0.4699105738175311 0.3943017891765944
31 -0.5031627214892470 -0.8715033980940043
32 0.2656516214009662 1.5065852108032667
33 0.3742336568932738 -2.1223845346278751
34 -1.4066938032409491 2.4364651379056212
35 2.6071563537030196 -2.1876639348137052
36 -3.5851642010512617 1.3048930541394301

1 0.2136039381609404 0.0000000000000000
2 0.1900433461077655 0.0893645283276201
3 0.1295805768614734 0.1524641483244400
4 0.0551876363483824 0.1777596430087896
5 -0.0123561349742845 0.1704519990727029
6 -0.0631106751037243 0.1435240234118874
7 -0.0966909073778188 0.1084505545109052
8 -0.1168298917580710 0.0716723675914216
9 -0.1273355034507535 0.0354863778605393
10 -0.1305808314668788 0.0000000000000000
11 -0.1273355034507535 -0.0354863778605393
12 -0.1168298917580710 -0.0716723675914216
13 -0.0966909073778188 -0.1084505545109052
14 -0.0631106751037243 -0.1435240234118873
15 -0.0123561349742845 -0.1704519990727029
16 0.0551876363483824 -0.1777596430087895
17 0.1295805768614734 -0.1524641483244400
18 0.1900433461077654 -0.0893645283276201
19 0.2136039381609404 0.0000000000000000

From: dbd on 9 Nov 2009 11:33
On Nov 9, 3:03 am, JochenH <jochenh...(a)t-online.de> wrote:
> I'm not entirely shure if this is the right place to ask, so if my post
> does not fit please ignore or suggest a better place.
>
> I want to do a FFT on real data using the FFTW3 library. The code i'm
> working on is based on hermitian matrices, therefore i need to use the
> r2c transform.
> As far as i understand it, the first 2n+1 entries with a c2c transform
> on real data should be the same as with a r2c transform.
> However, i can't reproduce the c2c output with a r2c transform.
> Thanks.
>
> Fortran source code:
>
> program kint
>
>    use fftw3
>
>    implicit none
>
>    complex(kind=8), allocatable, dimension(:) :: fr,fk,fk2
>
>    real(kind=8) :: dx,dk,a,pi
>    complex(kind=8) :: resultk
>    integer :: nx,i,plan
>
>    nx = 36
>    dx = 1.d0
>    pi  = 3.141592653589793d0
>    dk = 2.d0*pi/(dble(nx)*dx)
>    a = 0.2d0
>
>    allocate(fr(nx),fk(nx))
>
>    do i=1,nx
>      fr(i) = cmplx(exp(-a * (dx*(dble(i-1)-dble(nx)/2.d0)-2)**2),0.d0)
>    end do
>
> !**************** c2c *************
>
>    call dfftw_plan_dft_1d(plan,nx,fr,fk,FFTW_BACKWARD,FFTW_ESTIMATE)
>    call dfftw_execute(plan)
>    call dfftw_destroy_plan(plan)
>
>    do i=1,nx
>      write(*,'(i,f,f)') i, dreal(fk(i)),dimag(fk(i))
>    end do
>
>    fk = fk * (dx / sqrt(2.d0 * pi))
>
>    deallocate(fk)
>
>    allocate(fk2(nx/2+1))
>
> write(*,'(a)') ''
>
> !**************** c2r *************
>
>    call dfftw_plan_dft_r2c_1d(plan,nx,fr,fk2,FFTW_BACKWARD,FFTW_ESTIMATE)
>    call dfftw_execute(plan)
>    call dfftw_destroy_plan(plan)
>
>    do i=1,nx/2+1
>      write(*,'(i,f,f)') i, dreal(fk2(i)),dimag(fk2(i))
>    end do
>
>    deallocate(fr,fk2)
>
> end program
>
> Output:
>
>             1       3.9633273373814508       0.0000000000000000
>             2      -3.5851642010512617      -1.3048930541394301
>             3       2.6071563537030196       2.1876639348137057
>             4      -1.4066938032409493      -2.4364651379056217
>             5       0.3742336568932738       2.1223845346278751
>             6       0.2656516214009662      -1.5065852108032667
>             7      -0.5031627214892469       0.8715033980940043
>             8       0.4699105738175311      -0.3943017891765944
>             9      -0.3256090393846692       0.1185119983440311
>            10       0.1813788371283939       0.0000000000000000
>            11      -0.0826751513495624      -0.0300912942046938
>            12       0.0302950931975201       0.0254206015285206
>            13      -0.0082365236849697      -0.0142660775001120
>            14       0.0011041449519306       0.0062619171921073
>            15       0.0003949462596618      -0.0022398515421467
>            16      -0.0003770206872807       0.0006530189858746
>            17       0.0001780852757594      -0.0001494312892079
>            18      -0.0000658764105576       0.0000239770525833
>            19       0.0000347119594307      -0.0000000000000002
>            20      -0.0000658764105575      -0.0000239770525833
>            21       0.0001780852757595       0.0001494312892079
>            22      -0.0003770206872806      -0.0006530189858747
>            23       0.0003949462596617       0.0022398515421467
>            24       0.0011041449519307      -0.0062619171921072
>            25      -0.0082365236849697       0.0142660775001120
>            26       0.0302950931975202      -0.0254206015285206
>            27      -0.0826751513495624       0.0300912942046940
>            28       0.1813788371283939       0.0000000000000000
>            29      -0.3256090393846690      -0.1185119983440312
>            30       0.4699105738175311       0.3943017891765944
>            31      -0.5031627214892470      -0.8715033980940043
>            32       0.2656516214009662       1.5065852108032667
>            33       0.3742336568932738      -2.1223845346278751
>            34      -1.4066938032409491       2.4364651379056212
>            35       2.6071563537030196      -2.1876639348137052
>            36      -3.5851642010512617       1.3048930541394301
>
>             1       0.2136039381609404       0.0000000000000000
>             2       0.1900433461077655       0.0893645283276201
>             3       0.1295805768614734       0.1524641483244400
>             4       0.0551876363483824       0.1777596430087896
>             5      -0.0123561349742845       0.1704519990727029
>             6      -0.0631106751037243       0.1435240234118874
>             7      -0.0966909073778188       0.1084505545109052
>             8      -0.1168298917580710       0.0716723675914216
>             9      -0.1273355034507535       0.0354863778605393
>            10      -0.1305808314668788       0.0000000000000000
>            11      -0.1273355034507535      -0.0354863778605393
>            12      -0.1168298917580710      -0.0716723675914216
>            13      -0.0966909073778188      -0.1084505545109052
>            14      -0.0631106751037243      -0.1435240234118873
>            15      -0.0123561349742845      -0.1704519990727029
>            16       0.0551876363483824      -0.1777596430087895
>            17       0.1295805768614734      -0.1524641483244400
>            18       0.1900433461077654      -0.0893645283276201
>            19       0.2136039381609404       0.0000000000000000

What did you expect to be (or see to be) different and why?

Dale B. Dalrymple
From: JochenH on 9 Nov 2009 15:11
>
> What did you expect to be (or see to be) different and why?
>
I expected the first 17 numbers in arrays fk and the numbers in array
fk2 to be equal.
As far as i could understand, the difference between r2c and c2c fft is
simply that
the second half of the data in r2c transforms is omitted because it can
easily be calculated
from the first half by using the hermitian symmetry, that is by taking
the complex conjugate and reversing the order. However, that should not
affect the first half, which therefore should be the same for r2c and
c2c transforms. Am i confusing different things?

Jochen
From: dbd on 9 Nov 2009 16:43
On Nov 9, 12:11 pm, JochenH <jochenh...(a)t-online.de> wrote:
> > What did you expect to be (or see to be) different and why?
>
> I expected the first 17 numbers in arrays fk and the numbers in array
> fk2 to be equal.
> As far as i could understand, the difference between r2c and c2c fft is
> simply that
> the second half of the data in r2c transforms is omitted because it can
> easily be calculated
> from the first half by using the hermitian symmetry, that is by taking
> the complex conjugate and reversing the order. However, that should not
> affect the first half, which therefore should be the same for r2c and
> c2c transforms. Am i confusing different things?
>
> Jochen

Why are both transforms FFTW_BACKWORD?

Why don't you check for NULL plans?

Dale B. Dalrymple
 | 
Pages: 1
Prev: Suggestion for DS-CDMA simulation issue (just have a look)
Next: Practical differences between BCH and Reed-Solomon codes?