- c - 在位数组中找到第一个零
- linux - Unix 显示有关匹配两种模式之一的文件的信息
- 正则表达式替换多个文件
- linux - 隐藏来自 xtrace 的命令
我正在学习如何使用 Fortran 处理 FFTW 包。为了生成一个易于验证的示例,我计算了一个二维平面的功率谱,我用两个不同的叠加波填充它。这样,我就可以确切地知道功率谱中的峰值在哪里。
由于 FFTW 文档针对 C 的内容更为广泛,我首先在 C 中实现了该算法,这给了我相当满意的结果: “lambda1”和“lambda2”是已知的预期波长。蓝线是得到的功率谱。
我不知道去哪里寻找可能的错误。代码执行顺利。有人可以帮忙吗?
这是 C 代码,编译时使用:gcc stackexchange.c -o a.out -I/home/myname/.local/include -lfftw3 -lm -g -Wall
(gcc 5.4 .0)
#include <fftw3.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
int main(void)
{
// parameters
int Nx = 200;
int Ny = 100;
int nsamples = 200;
float pi = 3.1415926;
float physical_length_x = 20;
float physical_length_y = 10;
float lambda1 = 0.5;
float lambda2 = 0.7;
float dx = physical_length_x/Nx;
float dy = physical_length_y/Ny;
float dkx = 2*pi/physical_length_x;
float dky = 2*pi/physical_length_y;
// counters/iterators
int ind, i, j, ix, iy, ik;
// power spectra stuff
float *Pk, *distances_k, d, kmax;
float *Pk_field;
// FFTW vars and arrays
fftw_complex *in, *out;
fftw_plan my_plan;
// allocate arrays for input/output
in = (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * Nx * Ny);
out = (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * Nx * Ny);
// Create Plan
int n[2]; n[0] = Nx; n[1] = Ny;
my_plan = fftw_plan_dft(2, n, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
// fill up array with a wave
for (i=0; i<Nx; i++){
for (j=0; j<Ny; j++){
ind = i*Ny + j;
in[ind][0] = cos(2.0*pi/lambda1*i*dx) + sin(2.0*pi/lambda2*j*dy);
}
}
// execute fft
fftw_execute(my_plan); /* repeat as needed */
// Calculate power spectrum: P(kx, ky) = |F[delta(x,y)]|^2
Pk_field = malloc(sizeof(float)*Nx*Ny);
for (i=0; i<Nx; i++){
for (j=0; j<Ny; j++){
ind = i*Ny+j;
Pk_field[ind] = out[ind][0]*out[ind][0] + out[ind][1]*out[ind][1];
}
}
Pk = calloc(nsamples, sizeof(float));
distances_k = malloc(nsamples*sizeof(float));
kmax = sqrt(pow((Nx/2+1)*dkx, 2) + pow((Ny/2+1)*dky, 2));
for(i=0; i<nsamples; i++){
distances_k[i]= 1.0001*i/nsamples*kmax; // add a little more to make sure kmax will fit
}
// histogrammize P(|k|)
for (i=0; i<Nx; i++){
if (i<Nx/2+1)
{ ix = i; }
else
{ ix = -Nx+i; }
for (j=0; j<Ny; j++){
if (j<Ny/2+1)
{ iy = j; }
else
{ iy = -Ny+j; }
ind = i*Ny + j;
d = sqrt(pow(ix*dkx,2)+pow(iy*dky,2));
for(ik=0; ik<nsamples; ik++){
if (d<=distances_k[ik] || ik==nsamples){
break;
}
}
Pk[ik] += Pk_field[ind];
}
}
//-----------------------------------
// write arrays to file.
// can plot them with plot_fftw.py
//-----------------------------------
FILE *filep;
filep = fopen("./fftw_output_2d_complex.txt", "w");
for (i=0; i<nsamples; i++){
fprintf(filep, "%f\t%f\n", distances_k[i], Pk[i]);
}
fclose(filep);
printf("Finished! Written results to ./fftw_output_2d_complex.txt\n");
//----------------------
// deallocate arrays
//----------------------
fftw_destroy_plan(my_plan);
fftw_free(in); fftw_free(out);
return 0;
}
这是 Fortran 代码,编译时使用:gfortran stackexchange.f03 -o a.out -L/home/my_name/.local/lib/-I/home/my_name/.local/include/-lfftw3 -g -Wall
program use_fftw
use,intrinsic :: iso_c_binding
implicit none
include 'fftw3.f03'
integer, parameter :: dp=kind(1.0d0)
integer, parameter :: Nx = 200
integer, parameter :: Ny = 100
integer, parameter :: nsamples = 200
real(dp), parameter :: pi = 3.1415926d0
real(dp), parameter :: physical_length_x = 20.d0
real(dp), parameter :: physical_length_y = 10.d0
real(dp), parameter :: lambda1 = 0.5d0
real(dp), parameter :: lambda2 = 0.7d0
real(dp), parameter :: dx = physical_length_x/real(Nx,dp)
real(dp), parameter :: dy = physical_length_y/real(Ny,dp)
real(dp), parameter :: dkx = 2.d0 * pi / physical_length_x
real(dp), parameter :: dky = 2.d0 * pi / physical_length_y
integer :: i, j, ix, iy, ik
real(dp):: kmax, d
complex(dp), allocatable, dimension(:,:) :: arr_in, arr_out, Pk_field
real(dp), allocatable, dimension(:) :: Pk, distances_k
integer*8 :: my_plan
integer :: n(2) = (/Nx, Ny/)
allocate(arr_in( 1:Nx,1:Ny))
allocate(arr_out(1:Nx,1:Ny))
! call dfftw_plan_dft_2d(my_plan, Nx, Ny, arr_in, arr_out, FFTW_FORWARD, FFTW_ESTIMATE) ! doesn't help either
call dfftw_plan_dft(my_plan, 2, n, arr_in, arr_out, FFTW_FORWARD, FFTW_ESTIMATE)
! fill up wave
do i = 1, Nx
do j = 1, Ny
arr_in(i,j) =cmplx(cos(2.0*pi/lambda1*i*dx)+sin(2.0*pi/lambda2*j*dy) , 0.d0, kind=dp)
enddo
enddo
! execute fft
call dfftw_execute_dft(my_plan, arr_in, arr_out)
allocate(Pk_field(1:Nx, 1:Ny))
allocate(distances_k(1:nsamples))
allocate(Pk(1:nsamples))
Pk=0
distances_k=0
! Get bin distances
kmax = sqrt(((Nx/2+1)*dkx)**2+((Ny/2+1)*dky)**2)*1.001
do i = 1, nsamples
distances_k(i) = i*kmax/nsamples
enddo
! Compute P(k) field, distances field
do i = 1, Nx
do j = 1, Ny
Pk_field(i,j) = arr_out(i,j)*conjg(arr_out(i,j))
if (i<=Nx/2+1) then
ix = i
else
ix = -Nx+i
endif
if (j<=Ny/2+1) then
iy = j
else
iy = -Ny+j
endif
d = sqrt((ix*dkx)**2+(iy*dky)**2)
do ik=1, nsamples
if (d<=distances_k(ik) .or. ik==nsamples) exit
enddo
Pk(ik) = Pk(ik)+real(Pk_field(i,j))
enddo
enddo
! write file
open(unit=666,file='./fftw_output_2d_complex.txt', form='formatted')
do i = 1, nsamples
write(666, '(2E14.5,x)') distances_k(i), Pk(i)
enddo
close(666)
deallocate(arr_in, arr_out, Pk, Pk_field, distances_k)
call dfftw_destroy_plan(my_plan)
write(*,*) "Finished! Written results to fftw_output_2d_complex.txt"
end program use_fftw
为了您的方便,这是我用来绘制结果的简短 python 脚本:
#!/usr/bin/python3
#====================================
# Plots the results of the FFTW
# example programs.
#====================================
import numpy as np
import matplotlib.pyplot as plt
from sys import argv
from time import sleep
errormessage="""
I require an argument: Which output file to plot.
Usage: ./plot_fftw.py <case>
options for case:
1 fftw_1d_complex.txt
2 fftw_2d_complex.txt
3 fftw_1d_real.txt
4 fftw_2d_real.txt
Please select a case: """
#----------------------
# Hardcoded stuff
#----------------------
file_dict={}
file_dict['1'] = ('fftw_output_1d_complex.txt', '1d complex fftw')
file_dict['2'] = ('fftw_output_2d_complex.txt', '2d complex fftw')
file_dict['3'] = ('fftw_output_1d_real.txt', '1d real fftw')
file_dict['4'] = ('fftw_output_2d_real.txt', '2d real fftw')
lambda1=0.5
lambda2=0.7
#------------------------
# Get case from cmdline
#------------------------
case = ''
def enforce_integer():
global case
while True:
case = input(errormessage)
try:
int(case)
break
except ValueError:
print("\n\n!!! Error: Case must be an integer !!!\n\n")
sleep(2)
if len(argv) != 2:
enforce_integer()
else:
try:
int(argv[1])
case = argv[1]
except ValueError:
enforce_integer()
filename,title=file_dict[case]
#-------------------------------
# Read and plot data
#-------------------------------
k, Pk = np.loadtxt(filename, dtype=float, unpack=True)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_title("Power spectrum for "+title)
ax.set_xlabel("k")
ax.set_ylabel("P(k)")
# ax.plot(k, Pk, label='power spectrum')
ax.semilogx(k[k>0], Pk[k>0], label='power spectrum') # ignore negative k
ax.plot([2*np.pi/lambda1]*2, [Pk.min()-1, Pk.max()+1], label='expected lambda1')
ax.plot([2*np.pi/lambda2]*2, [Pk.min()-1, Pk.max()+1], label='expected lambda2')
ax.legend()
plt.show()
最佳答案
在直方图的计算中,循环的第一个索引对应于平均值 i。 e. “零频率”。
在C中,当i=0,j=0,ix=0,iy=0,d=0。那是正确的。对于fortran中数组的同一项,索引i=1,则ix=1,d不再为null。
那是零频率,但由于 C 和 Fortran 之间的索引差异,整个直方图被轻微污染。特别是,修改可能以不同方式影响正频率和负频率,从而触发图中观察到的双峰。
您能否尝试将 Fortran 中 ix 和 iy 的计算修改为:
if (i-1<Nx/2+1) then
ix = i-1
else
ix = -Nx+i-1
endif
if (j-1<Ny/2+1) then
iy = j-1
else
iy = -Ny+j-1
endif
关于c - fortran 2d-FFTW 与 C 2d-FFTW 结果不一致,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/50726456/
关闭。这个问题需要debugging details .它目前不接受答案。 编辑问题以包含 desired behavior, a specific problem or error, and th
我试图用这种形式简单地获取数字 28 integer+space+integer+integer+space+integer我试过这个正则表达式 \\s\\d\\d\\s 但我得到了两个数字11 和
最近一直在学习D语言。我一直对运行时感到困惑。 从我能收集到的关于它的信息中,(这不是很多)我知道它是一种有助于 D 的一些特性的运行时。像垃圾收集一样,它与您自己的程序一起运行。但是既然 D 是编译
想问一下这两个正则表达式有区别吗? \d\d\d 与 \d{3} 我已经在我的本地机器上使用 Java 和 Windows 操作系统对此进行了测试,两者都工作正常并且结果相同。但是,当在 linux
我正在学习 Go,而且我坚持使用 Go 之旅(exercise-stringer.go:https://tour.golang.org/methods/7)。 这是一些代码: type IPAddr
我在Java正则表达式中发现了一段令我困惑的代码: Pattern.compile( "J.*\\d[0-35-9]-\\d\\d-\\d\\d" ); 要编译的字符串是: String string
我在 ruby 代码上偶然发现了这个。我知道\d{4})\/(\d\d)\/(\d\d)\/(.*)/是什么意思,但是\1-\2-\3-\4 是什么意思? 最佳答案 \1-\2-\3-\4 是 b
我一直在努力解决这个问题,这让我很恼火。我了解 D 运行时库。它是什么,它做什么。我也明白你可以在没有它的情况下编译 D 应用程序。就像 XoMB 所做的那样。好吧,XoMB 定义了自己的运行时,但是
我有两个列表列表,子列表代表路径。我想找到所有路径。 List> pathList1 List> pathList2 当然是天真的解决方案: List> result = new ArrayList>
我需要使用 Regex 格式化一个字符串,该字符串包含数字、字母 a-z 和 A-Z,同时还包含破折号和空格。 从用户输入我有02-219 8 53 24 输出应该是022 198 53 24 我正在
目标是达到与this C++ example相同的效果: 避免创建临时文件。我曾尝试将 C++ 示例翻译为 D,但没有成功。我也尝试过不同的方法。 import std.datetime : benc
tl;dr:你好吗perfect forwarding在 D? 该链接有一个很好的解释,但例如,假设我有这个方法: void foo(T)(in int a, out int b, ref int c
有什么方法可以在 D 中使用abstract auto 函数吗? 如果我声明一个类如下: class MyClass { abstract auto foo(); } 我收到以下错误: mai
有没有人为内存中重叠的数组切片实现交集?算法在没有重叠时返回 []。 当 pretty-print (使用重叠缩进)内存中重叠的数组切片时,我想要这个。 最佳答案 如果您确定它们是数组,那么只需取 p
我已经开始学习 D,但我在使用 Andrei Alexandrescu 所著的 The D Programming Language 一书中提供的示例时遇到了一些麻烦。由于 int 和 ulong 类
如何创建一个不可变的类? 我的目标是创建一个实例始终不可变的类。现在我只是用不可变的方法和构造函数创建了一个“可变”类。我将其称为 mData,m 表示可变。然后我创建一个别名 alias immut
不久前我买了《The D Programming Language》。好书,很有教育意义。但是,我在尝试编译书中列出的语言功能时遇到了麻烦:扩展函数。 在这本书中,Andrei 写了任何可以像这样调用
我在 D http://www.digitalmars.com/d/2.0/lazy-evaluation.html 中找到了函数参数的惰性求值示例 我想知道如何在 D 中实现可能的无限数据结构,就像
这个问题在这里已经有了答案: 12 年前关闭。 Possible Duplicate: Could anyone explain these undefined behaviors (i = i++
当前是否可以跨模块扫描/查询/迭代具有某些属性的所有函数(或类)? 例如: source/packageA/something.d: @sillyWalk(10) void doSomething()
我是一名优秀的程序员,十分优秀!