2015-02-04 2 views
1

불연속 푸리에 변환을 계산하는 서브 루틴이 포함 된 Fortran 90의 FFTPACK5.1을 사용하는 데 문제가 있습니다. 나는 그것을 설치하고 루틴을 사용할 수 있지만 모든 것이 주파수가있는 간단한 사인파인지 확인하는 경우 나는 A가 아닌 (주파수 공간의 스펙트럼에서) 0이 아닌 계수를 얻지 만, 2A. 스펙트럼에 변화가 있고 나는 왜 이해하지 못합니다. 나는 정확히 주파수 축 스텝을 계산한다는 것을 거의 확신한다 (그러나 나는 의문이다) :fftpack5.1을 사용한 푸리에 변환 관련 문제

N은 원래의 사인파의 포인트 수이고 Fech 나의 샘플 주파수는 주파수 축 스텝을 df로 계산한다. (i) = Fech (i-1)/N이다.

저는 rfft1f 루틴을 사용하고 있습니다. 따라서 누군가가 경험을 쌓고 내 문제를 알고 있다면 여기에서 잘못된 점을 이해하는 데 정말 도움이 될 것입니다. 여기

내 코드입니다 :

! n: number of samples in the discret signal 
integer (kind = 4), parameter :: n = 4096 
real, parameter :: deuxpi=6.283185307 
!frequence is the frequence of the signal 
!fech is the frequence of sampling 
real :: frequence,fech 
integer :: kk 
! r is the signal i want to process 
! t is the built time and f is the built frequency 
real (kind = 4) r(n),t(n),f(n) 


!Arrays routines need to work (internal recipe): 
real (kind = 4), allocatable, dimension (:) :: work 
real (kind = 4), allocatable, dimension (:) :: wsave 

!size of arrays wsave and work for internal recipe 
lensav = n + int (log (real (n, kind = 4))/log (2.0E+00)) + 4 
lenwrk = n 
allocate (work(1:lenwrk)) 
allocate (wsave(1:lensav)) 


! initializes rttft1f, wsave array 
call rfft1i (n, wsave, lensav, ier) 


frequence=0.5 
fech=20 

! I built here the signal 
do kk=1,n 
t (kk) = (kk-1)/(fech) 
f (kk) = fech*(kk-1)/n 
r (kk) = sin(deuxpi * t(kk) * frequence ) 
end do 


!and I call the rfft1f to build the Discrete Fourier Transform: 
call rfft1f (n, inc, r, lenr, wsave, lensav, work, lenwrk, ier) 

!I get back r which contains now the fourier coefficients and plot it 

I 0.5 (CF 코드)의 주파수에서 디랙이 간단한 사인파로 기대하고있어 대신 난에, 1에서 디랙를 얻을 수 주파수 도메인. 실수 값 시퀀스로 변환 이산 푸리에 변환을 계산하는 루틴의 전형으로

Spectrum

+1

코드, 예상되는 동작 및 실제 동작을 보여줍니다. –

답변

0

은, 그 결과 복잡한 가치 스펙트럼이 절반 만 반환됩니다 : 게다가, 스펙트럼 내가 무엇을 얻을 여기에 ... 이상한 보인다 스펙트럼. 값을 원래의 N 요소 배열에 맞추려면 첫 번째 값의 실수 부분 (항상 실제 값) 만 반환됩니다. 마찬가지로 n의 짝수 값의 경우 n/2-th 복합 값 (항상 실제 값)의 실수 부분이 반환됩니다. 다른 모든 복소수 값의 경우 실수 부와 허수 부는 인터리브됩니다.

r(1)  -> 0 
r(2), r(3) -> Delta 
r(4), r(5) -> 2*Delta 
... 
r(n)  -> (n/2)*Delta 

홀수 n 경우 :

이 너무도 n 경우, 해당 주파수에 의해 주어진다

r(1)  -> 0 
r(2), r(3) -> Delta 
r(4), r(5) -> 2*Delta 
... 
r(n-1),r(n) -> ((n-1)/2) * Delta 

델타 fech/n 같이 주어진다.

복잡한 값을 플롯하려면 실제 (인덱스 1,2,4,6, ...) & 허수 (인덱스 3,5,7, ...) 부분을 다음과 같이 플롯해야합니다. 두 개의 개별 그래프 또는 진폭 & 단계 (두 개의 별도 그래프로).

+0

답변을 주셔서 감사합니다. 원래 n 개의 배열 크기가 -n/2에서 n/2로 확장되어 있기 때문에 대칭을 알기 때문에 긍정적 인 부분을 알기 때문에 긍정적 인 부분을 얻었습니다. 2로 나눕니 까? – Moonpalacio

+0

'n'포인트 FFT는 'n'시간 도메인 값을 취하여 'n'복소 주파수 도메인 값으로 변환합니다. 빈도 단계는 'fech/n'입니다. 스펙트럼의 절반을 매핑 할 때'fech/2'를 다루는'n/2' 포인트가 있으므로'fech/2/(n/2) = fech/n'도 다시 나타납니다. 2로 나누기위한 주된 이유는 복잡한 값을 부동 소수점 숫자의 배열 (실수/허수 분수 교대)에 저장하기 때문입니다. – SleuthEye