XPEHOMETP
Silver Member | Редактировать | Профиль | Сообщение | Цитировать | Сообщить модератору
Andrew10
В том-то самый смак, чтобы понимать, как эта штука работает! А с пониманием по поводу FFTW - проблематично. Как и с пониманием любых программ, написанных по методу псевдо-динамического распределения памяти. Когда объявляют некий рабочий массив, в качестве такой братской могилы, и успешно валят туда и вещественные данные, и целочисленные. Да там вообще волосы дыбом встанут, если в этом FFTW попытаться хоть чуток разобраться.
По поводу программы по разложению Фурье. Выкладывал на одном из форумов. Это несколько доработанный вариант по сравнению с тем, по которому я реально работал. Но я его не компилировал. Просто отказался потом от этого в пользу вариантов, найденых в сети.
Код: subroutine foft(dat,nft,ftcoef,free,start)
c---------------------------------------
c dat - входной массив данных удвоенной точности
c nft - число точек в массиве dat
c ftcoef - массив коэффициентов разложения Фурье размерности (nft/2 * 2) - на выход
c Один столбец у него для косинусов, другой для синусов
c Тут для простоты примем, что nft - четное число, тогда длина ftcoef будет ровно nft/2
c free - свободный член разложения Фурье (умные дяди программеры стандартно запихивают его
c в выходной массив, но мы так делать не будем, чтобы все было как в формуле) - на выход
c
implicit none
integer nft,j,n,m,pos
real*8 datt(nft),trig(nft,2),ftcoef(nft/2,2),free,x,twopi
parameter (twopi = 6.28318530717959)
c cчитаем свободный член как среднее арифметическое:
c----------------------------------------------------------
x = 0.0
do 5 j = 1, nft
5 x = x + dat(j)
free = x/nft
c Инициируем trig(nft,2) - рабочий массив косинусов-синусов:
c----------------------------------------------------------
trig(1,1) = 1.0d
trig(1,2) = 0.0d
do 10 j = 2,nft
trig(j,1) = cos(twopi*(j-1)/nft)
10 trig(j,2) = sin(twopi*(j-1)/nft)
do 30 m = 1,2
do 30 j = 1, nft/2-1
x = 0.0
do 20 n = 1, nft
pos = mod((n-1)*j,nft)+1
20 x = x + trig(pos,m)*dat(n)
ftcoef(j,m) = 2*x/nft
30 continue
ftcoef(nft/2,2) = 0.0d
x = 0.0d
do 40 n = 1, nft
pos = mod((n-1)*nft/2,nft)+1
40 x = x + trig(pos,1)*dat(n)
ftcoef(nft/2,1) = 2*x/nft
c Boт поскольку коэффициент ftcoef(nft/2,2) всегда и везде будет 0, ушлые программеры на
c его место запихивают свободный член разложения. Но мы так не сделаем.
return
end
|
|
. 
