FFT：fortran 与 python
FFT: fortran vs. python
我有计算离散信号（具有两个不同频率的双正弦信号）的 FFT 的 Fortran 代码，摘自：
y = 0.5*np.sin(2 * np.pi * ff1 * t) + 0.1*np.sin(2 * np.pi * ff2 * t)
当我用 Fortran 代码计算 FFT 并将其与用 python 计算的结果进行比较时，我可以看到：
1. 两个图表中选秀权的分散是由于四舍五入？我能以某种方式消除或减少它吗？
python中使用的代码是：
import numpy as np
import matplotlib.pyplot as plt
from scipy import fft
Fs = 2048 # sampling rate = number of lines in the input file
Ts = 1.0/Fs # sampling interval
data = np.loadtxt('input.dat')
t = data[:,0]
y = data[:,1]
plt.subplot(2,1,1)
plt.plot(t,y,'ro')
plt.xlabel('time')
plt.ylabel('amplitude')
plt.subplot(2,1,2)
n = len(y) # length of the signal
k = np.arange(n)
T = n/Fs # equal 1
frq = k/T # two sides frequency range
freq = frq[range(n/2)] # one side frequency range
Y = np.fft.fft(y)/n # fft computing and normalization
Y = Y[range(n/2)]
plt.plot(freq, abs(Y), 'r-')
plt.axis([0, 20, 0, .4])
plt.xlabel('freq (Hz)')
plt.ylabel('|Y(freq)|')
plt.show()
2. 我的 fortran 程序中的 fft 幅度不等于 python 中计算的幅度，是 |Y(freq)|在 python 中计算等于：
ABS (AR(I)**2+ AI(I)**2) / n_tot
其中 AR 和 AI 是信号的实部和虚部，n_tot 是总点数。
使用的 Fortran 代码在这里：
PROGRAM fft
IMPLICIT NONE
INTEGER, PARAMETER :: N=2048 ! tot_num of points
INTEGER, PARAMETER :: M=11 !! this is the exp in subroutine: N1 = 2**M
INTEGER :: I,J
REAL(8) :: PI,F1,T
REAL(8), DIMENSION (N) :: AR,AI,O,time
!
PI = 4.D0*DATAN(1.D0) ; F1 = 1.d0/SQRT(real(N))
open(unit=6,file="input.dat")
do i = 1,n
read(6,*) time(i),ar(i)
end do
!
DO I = 1, N
AI(I) = 0.D0
END DO
CALL FFT (AR,AI,N,M)
!
OPEN(unit=6,file="output.dat")
DO I = 1, 20
O(I) = I-1
AR(I) = (F1*AR(I))**2+(F1*AI(I))**2 !! this is for the |y(freq)|
AR(I) = dabs(AR(I)) !! absolute value of y(freq)
WRITE(6,"(3F16.10)") O(I),AR(I)
END DO
CLOSE(6)
END PROGRAM fft
!
SUBROUTINE FFT(AR,AI,N,M)
!
! An example of the fast Fourier transform subroutine with N = 2**M.
! AR and AI are the real and imaginary part of data in the input and
! corresponding Fourier coefficients in the output.
! Copyright (c) Tao Pang 1997.
!
IMPLICIT NONE
INTEGER, INTENT (IN) :: N,M
INTEGER :: N1,N2,I,J,K,L,L1,L2
REAL(8) :: PI,A1,A2,Q,U,V
REAL(8), INTENT (INOUT), DIMENSION (N) :: AR,AI
!
PI = 4.D0*ATAN(1.D0)
N2 = N/2
!
N1 = 2**M
IF(N1.NE.N) STOP 'Indices do not match'
!
! Rearrange the data to the bit reversed order
!
L = 1
DO K = 1, N-1
IF (K.LT.L) THEN
A1 = AR(L)
A2 = AI(L)
AR(L) = AR(K)
AR(K) = A1
AI(L) = AI(K)
AI(K) = A2
END IF
J = N2
DO WHILE (J.LT.L)
L = L-J
J = J/2
END DO
L = L+J
END DO
!
! Perform additions at all levels with reordered data
!
L2 = 1
DO L = 1, M
Q = 0.D0
L1 = L2
L2 = 2*L1
DO K = 1, L1
U = DCOS(Q)
V = -DSIN(Q)
Q = Q + PI/L1
DO J = K, N, L2
I = J + L1
A1 = AR(I)*U-AI(I)*V
A2 = AR(I)*V+AI(I)*U
AR(I) = AR(J)-A1
AR(J) = AR(J)+A1
AI(I) = AI(J)-A2
AI(J) = AI(J)+A2
END DO
END DO
END DO
END SUBROUTINE FFT
在 Python 脚本中实际上是 |Y(freq)|与实际信号幅度的一半相关 y 这在 Fortran 程序中是不正确的。
在您的 Python 代码的以下几行中
Y = np.fft.fft(y)/n
plot(freq, abs(Y), 'r-')
Y是FFT得到的复数数组，所以abs(Y)是由|Y[i]|和|z|组成的数组= sqrt( Re{z}^2 + Im{z}^2 ).但是，在您的 Fortran 代码的以下几行中
AR(I) = (F1*AR(I))**2+(F1*AI(I))**2 !! this is for the |y(freq)|
AR(I) = dabs(AR(I)) !! absolute value of y(freq)
（其中 F1 = 1/sqrt(N)），实部和虚部平方求和后缺少 sqrt()。因此，将这些行替换为
AR(I) = sqrt( AR(I)**2 + AI(I)**2 ) / N
给出了预期的结果（假设 FFT() 与 np.fft.fft() 具有相同的归一化）。例如，AR(6)和AR(9)变为0.25和0.05，这与Python得到的结果一致。
下面是一些代码比较（只是为了好玩！），它们给出了所有相同的结果。
Python:
from numpy import pi, sin
N = 2048
t = np.linspace( 0.0, 1.0, N+1 )[:-1]
y = 0.5 * sin(2 * pi * 5 * t) + 0.1 * sin(2 * pi * 8 * t)
z = np.fft.fft( y )
for i in range( 10 ) + range( N-9, N ):
print ("%4d" + "%16.10f" * 3) % \
( i, z[i].real, z[i].imag, abs( z[i] ) / N )
Fortran:
time = [( dble(i-1) / dble(N), i=1,N )]
AR = 0.5d0 * sin(2 * pi * 5 * time) + 0.1d0 * sin(2 * pi * 8 * time)
AI = 0.0d0
call FFT (AR,AI,N,M)
do i = 1, N
if ( i > 10 .and. i < N-8 ) cycle
print "(i4, 3f16.10)", &
i-1, AR(i), AI(i), sqrt( AR(i)**2 + AI(i)**2 ) / N
enddo
朱莉娅：
N = 2048
t = linspace( 0.0, 1.0, N+1 )[1:N]
y = 0.5 * sin(2 * pi * 5 * t) + 0.1 * sin(2 * pi * 8 * t)
z = fft( y )
for i in [ 1:10 ; N-8:N ]
@printf( "%4d%16.10f%16.10f%16.10f\n",
i-1, real( z[i] ), imag( z[i] ), abs( z[i] ) / N )
end
我有计算离散信号（具有两个不同频率的双正弦信号）的 FFT 的 Fortran 代码，摘自：
y = 0.5*np.sin(2 * np.pi * ff1 * t) + 0.1*np.sin(2 * np.pi * ff2 * t)
当我用 Fortran 代码计算 FFT 并将其与用 python 计算的结果进行比较时，我可以看到：
1. 两个图表中选秀权的分散是由于四舍五入？我能以某种方式消除或减少它吗？
python中使用的代码是：
import numpy as np
import matplotlib.pyplot as plt
from scipy import fft
Fs = 2048 # sampling rate = number of lines in the input file
Ts = 1.0/Fs # sampling interval
data = np.loadtxt('input.dat')
t = data[:,0]
y = data[:,1]
plt.subplot(2,1,1)
plt.plot(t,y,'ro')
plt.xlabel('time')
plt.ylabel('amplitude')
plt.subplot(2,1,2)
n = len(y) # length of the signal
k = np.arange(n)
T = n/Fs # equal 1
frq = k/T # two sides frequency range
freq = frq[range(n/2)] # one side frequency range
Y = np.fft.fft(y)/n # fft computing and normalization
Y = Y[range(n/2)]
plt.plot(freq, abs(Y), 'r-')
plt.axis([0, 20, 0, .4])
plt.xlabel('freq (Hz)')
plt.ylabel('|Y(freq)|')
plt.show()
2. 我的 fortran 程序中的 fft 幅度不等于 python 中计算的幅度，是 |Y(freq)|在 python 中计算等于：
ABS (AR(I)**2+ AI(I)**2) / n_tot
其中 AR 和 AI 是信号的实部和虚部，n_tot 是总点数。
使用的 Fortran 代码在这里：
PROGRAM fft
IMPLICIT NONE
INTEGER, PARAMETER :: N=2048 ! tot_num of points
INTEGER, PARAMETER :: M=11 !! this is the exp in subroutine: N1 = 2**M
INTEGER :: I,J
REAL(8) :: PI,F1,T
REAL(8), DIMENSION (N) :: AR,AI,O,time
!
PI = 4.D0*DATAN(1.D0) ; F1 = 1.d0/SQRT(real(N))
open(unit=6,file="input.dat")
do i = 1,n
read(6,*) time(i),ar(i)
end do
!
DO I = 1, N
AI(I) = 0.D0
END DO
CALL FFT (AR,AI,N,M)
!
OPEN(unit=6,file="output.dat")
DO I = 1, 20
O(I) = I-1
AR(I) = (F1*AR(I))**2+(F1*AI(I))**2 !! this is for the |y(freq)|
AR(I) = dabs(AR(I)) !! absolute value of y(freq)
WRITE(6,"(3F16.10)") O(I),AR(I)
END DO
CLOSE(6)
END PROGRAM fft
!
SUBROUTINE FFT(AR,AI,N,M)
!
! An example of the fast Fourier transform subroutine with N = 2**M.
! AR and AI are the real and imaginary part of data in the input and
! corresponding Fourier coefficients in the output.
! Copyright (c) Tao Pang 1997.
!
IMPLICIT NONE
INTEGER, INTENT (IN) :: N,M
INTEGER :: N1,N2,I,J,K,L,L1,L2
REAL(8) :: PI,A1,A2,Q,U,V
REAL(8), INTENT (INOUT), DIMENSION (N) :: AR,AI
!
PI = 4.D0*ATAN(1.D0)
N2 = N/2
!
N1 = 2**M
IF(N1.NE.N) STOP 'Indices do not match'
!
! Rearrange the data to the bit reversed order
!
L = 1
DO K = 1, N-1
IF (K.LT.L) THEN
A1 = AR(L)
A2 = AI(L)
AR(L) = AR(K)
AR(K) = A1
AI(L) = AI(K)
AI(K) = A2
END IF
J = N2
DO WHILE (J.LT.L)
L = L-J
J = J/2
END DO
L = L+J
END DO
!
! Perform additions at all levels with reordered data
!
L2 = 1
DO L = 1, M
Q = 0.D0
L1 = L2
L2 = 2*L1
DO K = 1, L1
U = DCOS(Q)
V = -DSIN(Q)
Q = Q + PI/L1
DO J = K, N, L2
I = J + L1
A1 = AR(I)*U-AI(I)*V
A2 = AR(I)*V+AI(I)*U
AR(I) = AR(J)-A1
AR(J) = AR(J)+A1
AI(I) = AI(J)-A2
AI(J) = AI(J)+A2
END DO
END DO
END DO
END SUBROUTINE FFT
在 Python 脚本中实际上是 |Y(freq)|与实际信号幅度的一半相关 y 这在 Fortran 程序中是不正确的。
在您的 Python 代码的以下几行中
Y = np.fft.fft(y)/n
plot(freq, abs(Y), 'r-')
Y是FFT得到的复数数组，所以abs(Y)是由|Y[i]|和|z|组成的数组= sqrt( Re{z}^2 + Im{z}^2 ).但是，在您的 Fortran 代码的以下几行中
AR(I) = (F1*AR(I))**2+(F1*AI(I))**2 !! this is for the |y(freq)|
AR(I) = dabs(AR(I)) !! absolute value of y(freq)
（其中 F1 = 1/sqrt(N)），实部和虚部平方求和后缺少 sqrt()。因此，将这些行替换为
AR(I) = sqrt( AR(I)**2 + AI(I)**2 ) / N
给出了预期的结果（假设 FFT() 与 np.fft.fft() 具有相同的归一化）。例如，AR(6)和AR(9)变为0.25和0.05，这与Python得到的结果一致。
下面是一些代码比较（只是为了好玩！），它们给出了所有相同的结果。
Python:
from numpy import pi, sin
N = 2048
t = np.linspace( 0.0, 1.0, N+1 )[:-1]
y = 0.5 * sin(2 * pi * 5 * t) + 0.1 * sin(2 * pi * 8 * t)
z = np.fft.fft( y )
for i in range( 10 ) + range( N-9, N ):
print ("%4d" + "%16.10f" * 3) % \
( i, z[i].real, z[i].imag, abs( z[i] ) / N )
Fortran:
time = [( dble(i-1) / dble(N), i=1,N )]
AR = 0.5d0 * sin(2 * pi * 5 * time) + 0.1d0 * sin(2 * pi * 8 * time)
AI = 0.0d0
call FFT (AR,AI,N,M)
do i = 1, N
if ( i > 10 .and. i < N-8 ) cycle
print "(i4, 3f16.10)", &
i-1, AR(i), AI(i), sqrt( AR(i)**2 + AI(i)**2 ) / N
enddo
朱莉娅：
N = 2048
t = linspace( 0.0, 1.0, N+1 )[1:N]
y = 0.5 * sin(2 * pi * 5 * t) + 0.1 * sin(2 * pi * 8 * t)
z = fft( y )
for i in [ 1:10 ; N-8:N ]
@printf( "%4d%16.10f%16.10f%16.10f\n",
i-1, real( z[i] ), imag( z[i] ), abs( z[i] ) / N )
end