本来是想正余弦变换都做的,根据前人研究:正弦变换的结果比余弦更精确。所以本代码只做了正弦变化,如果读者想实现余弦变换,只需要把对应的正弦变换系数改为余弦变换系数即可。在末尾本文会给出250点正余弦变换系数。
代码如下:
!.. 本代码将实现正余弦变换
!.. 测试函数
!.. ∫exp(-w)sin(wt)dw
!.. 正弦变换离散表达式
!.. f(t)=∑{K(w)*wsl} * coef
!.. coef = 1/t*sqrt(pi/2)
Module testSineTransform
implicit none
integer, parameter :: nt = 46
real(kind=8), parameter :: pi = acos(-1.d0)
real(kind=8) :: t(nt), t1 = 5, t2 = 50 !.. 单位:s
real(kind=8) :: wsl(-149:100), w, delta
contains
subroutine cal_SineTransform
implicit none
integer :: i, j, fileid
real(kind=8) :: s1(nt), s2(nt) !.. s1为数值解, s2为解析解
do i = 1, nt
t(i) = t1 + dble(i-1) * ( t2 - t1 ) / dble( nt - 1)
s2(i) = func2( t(i) )!.. 解析表达式
end do
open( newunit = fileid, file = "sin250.txt" )
do i = -149, 100
read( fileid,* ) wsl(i)
end do
close( fileid )
!.. 计算正弦变换
delta = log(10.d0) / 20.d0 !.. 理论推导而来,无需更改
s1 = 0.d0
open( newunit = fileid, file = "result.dat" )
do i = 1, nt
do j = -149, 100
w = exp( dble(j)*delta ) / t(i) !.. 理论推导而来,无需更改
s1(i) = s1(i) + func1( w ) * wsl(j)
end do
s1(i) = s1(i) / t(i) * sqrt( pi/2.d0 ) !.. 1/t*sqrt(pi/2)为正弦变换出来的系数
write( fileid, '(*(g0,3x))' ) t(i), s1(i), s2(i), abs( s1(i)-s2(i) )/s2(i)*100
end do
close( fileid )
end subroutine cal_SineTransform
real(kind=8) function func1(w) !.. 核函数
implicit none
real(kind=8) :: w
func1 = exp(-1.d0*w)
end function func1
real(kind=8) function func2(t) !.. 解析解
implicit none
real(kind=8) :: t
func2 = t / ( t*t + 1.d0 )
end function func2
End module testSineTransform
Program main
use testSineTransform
implicit none
call cal_SineTransform
End program main
下面给出250点正弦变换系数
正弦变换系数:
1.156447054648637e-016
1.455880584491685e-016
1.832845064354326e-016
2.307415227397174e-016
2.904863665331116e-016
3.657006686082886e-016
4.603898648210914e-016
5.795965001557580e-016
7.296687626330448e-016
9.185985474711453e-016
1.156447054648634e-015
1.455880584491681e-015
1.832845064354321e-015
2.307415227397166e-015
2.904863665331103e-015
3.657006686082864e-015
4.603898648210901e-015
5.795965001557501e-015
7.296687626330329e-015
9.185985474711278e-015
1.156447054648612e-014
1.455880584491647e-014
1.832845064354265e-014
2.307415227397080e-014
2.904863665330965e-014
3.657006686082646e-014
4.603898648210535e-014
5.795965001556953e-014
7.296687626329462e-014
9.185985474709900e-014
1.156447054648389e-013
1.455880584491301e-013
1.832845064353717e-013
2.307415227396208e-013
2.904863665329599e-013
3.657006686080477e-013
4.603898648207092e-013
5.795965001551463e-013
7.296687626320796e-013
9.185985474696154e-013
1.156447054646209e-012
1.455880584487839e-012
1.832845064348231e-012
2.307415227387515e-012
2.904863665315808e-012
3.657006686058640e-012
4.603898648172482e-012
5.795965001496609e-012
7.296687626233827e-012
9.185985474558366e-012
1.156447054624371e-011
1.455880584453228e-011
1.832845064293375e-011
2.307415227300576e-011
2.904863665178019e-011
3.657006685840241e-011
4.603898647826370e-011
5.795965000948059e-011
7.296687625364435e-011
9.185985473180428e-011
1.156447054405990e-010
1.455880584107116e-010
1.832845063744824e-010
2.307415226431180e-010
2.904863663800132e-010
3.657006683656435e-010
4.603898644365244e-010
5.795964995462562e-010
7.296687616670500e-010
9.185985459401473e-010
1.156447052222165e-009
1.455880580645994e-009
1.832845058259313e-009
2.307415217737235e-009
2.904863650021138e-009
3.657006661818234e-009
4.603898609754006e-009
5.795964940607466e-009
7.296687529731007e-009
9.185985321611760e-009
1.156447030383951e-008
1.455880546034779e-008
1.832845003404190e-008
2.307415130797806e-008
2.904863512231297e-008
3.657006443436261e-008
4.603898263641547e-008
5.795964392056795e-008
7.296686660335758e-008
9.185983943714872e-008
1.156446812001731e-007
1.455880199922772e-007
1.832844454852834e-007
2.307414261403871e-007
2.904862134332549e-007
3.657004259617825e-007
4.603894802516352e-007
5.795958906554214e-007
7.296677966382738e-007
9.185970164759857e-007
1.156444628179741e-006
1.455876738807346e-006
1.832838969341425e-006
2.307405567480775e-006
2.904848355357452e-006
3.656982421491810e-006
4.603860191324607e-006
5.795904051742752e-006
7.296591027122000e-006
9.185832376013115e-006
1.156422790075531e-005
1.455842127958865e-005
1.832784114711636e-005
2.307318629427335e-005
2.904710567601305e-005
3.656764044807223e-005
4.603514087536846e-005
5.795355521527330e-005
7.295721667410295e-005
9.184454558911375e-005
1.156204422278983e-004
1.455496047235377e-004
1.832235621566062e-004
2.306449358701818e-004
2.903332902825535e-004
3.654580713105382e-004
4.600053900083832e-004
5.789871945981126e-004
7.287031465075219e-004
9.170683245122596e-004
1.154022098295481e-003
1.452037965023791e-003
1.826756086858394e-003
2.297767484717661e-003
2.889577754406488e-003
3.632790471922652e-003
4.565537633632412e-003
5.735207473181395e-003
7.200470157233521e-003
9.033651050110031e-003
1.132334380467910e-002
1.417727733464931e-002
1.772499346509631e-002
2.212023128150156e-002
2.754164273887034e-002
3.419148067351661e-002
4.228846080133485e-002
5.205416439533162e-002
6.368324102630288e-002
7.729769668970302e-002
9.286191525747088e-002
1.100669688846222e-001
1.281325720732841e-001
1.455759302997095e-001
1.598450602566018e-001
1.675994998118554e-001
1.619683745893242e-001
1.351613440916820e-001
8.220710085812771e-002
-5.378512256133189e-003
-1.293051042068928e-001
-2.723112991886841e-001
-3.978468461558793e-001
-4.394993669496010e-001
-3.176425577428793e-001
1.665756359292827e-002
4.630044483056220e-001
7.298213861736705e-001
3.957559489684660e-001
-4.902404478827365e-001
-1.020512348109581e+000
5.251407980261734e-002
1.253864451347815e+000
-1.317577223583377e-001
-1.549527153893942e+000
1.822344193602740e+000
-1.058204213789976e+000
2.629725097214586e-001
1.625578365785221e-001
-2.871269632839427e-001
2.688993168200287e-001
-2.085718537012279e-001
1.489333044554680e-001
-1.023510671993418e-001
6.920616476794303e-002
-4.650126535306088e-002
3.114481903045837e-002
-2.081379273802720e-002
1.390327480927853e-002
-9.298886523349949e-003
6.224588805575879e-003
-4.160927227996026e-003
2.775386859464604e-003
-1.851896921274413e-003
1.239790988028041e-003
-8.315445781273338e-004
5.558587783122442e-004
-3.696533331522879e-004
2.460602463017792e-004
-1.651618293595525e-004
1.114024919622099e-004
-7.389805836530149e-005
4.938931786488272e-005
-3.300905019751596e-005
2.206139792679832e-005
-1.474460111811994e-005
9.854464474452981e-006
-6.586171392444854e-006
4.401827590217889e-006
-2.941934696111318e-006
1.966224160030928e-006
-1.314113957933395e-006
8.782800707768828e-007
-5.869931432254631e-007
3.923132969292148e-007
-2.622002057839472e-007
1.752399127209358e-007
-1.171205297822945e-007
7.827679370240814e-008
-5.231581895777865e-008
3.496495939305279e-008
-2.336861793837973e-008
1.561827366110047e-008
-1.043837820431878e-008
6.976426582137068e-009
-4.662652272535390e-009
3.116255286086561e-009
-2.082730266047119e-009
1.391980105248554e-009
-9.303214367194074e-010
6.217746736151215e-010
-4.155593212088843e-010
2.777365447188800e-010
-1.856235303493902e-010
1.240603574666303e-010
-8.291498532420789e-011
5.541572611672626e-011
-3.703676348776440e-011
2.475329560347511e-011
-1.654371455635934e-011
1.105689099773430e-011
-7.389805845555023e-012
4.938931788889779e-012
-3.300905020390523e-012
2.206139792849932e-012
-1.474460111857239e-012
9.854464474573356e-013
-6.586171392476662e-013
4.401827590226633e-013
-2.941934696113489e-013
余弦变换系数:
250点余弦变换滤波系数(王华军, 2004)
n=-149~-100 n=-99~-50 n=-49~0 n=1~50 n=51~100
3.25931064E-09 1.03068452E-06 3.25929012E-04 4.47184906E-02 -1.61220740E-04
3.65700668E-09 1.15644705E-06 3.65697770E-04 3.54831893E-02 1.08092999E-04
4.10322899E-09 1.29755494E-06 4.10318805E-04 2.04513710E-02 -7.22432985E-05
4.60389865E-09 1.45588058E-06 4.60384082E-04 -2.05568659E-03 4.82833692E-05
5.16565924E-09 1.63352488E-06 5.16557757E-04 -3.36444840E-02 -3.22698961E-05
5.79596500E-09 1.83284506E-06 5.79584963E-04 -7.63113734E-02 2.15673888E-05
6.50317969E-09 2.05648598E-06 6.50301672E-04 -1.26643432E-01 -1.44144331E-05
7.29668762E-09 2.30741523E-06 7.29645743E-04 -1.86274938E-01 9.63379869E-06
8.18701817E-09 2.58896246E-06 8.18669301E-04 -2.46135314E-01 -6.43869076E-06
9.18598547E-09 2.90486366E-06 9.18552617E-04 -2.90298367E-01 4.30325981E-06
1.03068452E-08 3.25931064E-06 1.03061964E-03 -2.98807588E-01 -2.87605752E-06
1.15644705E-08 3.65700668E-06 1.15635541E-03 -2.44850614E-01 1.92219555E-06
1.29755494E-08 4.10322898E-06 1.29742549E-03 -9.95001792E-02 -1.28468770E-06
1.45588058E-08 4.60389864E-06 1.45569774E-03 1.35603930E-01 8.58613212E-07
1.63352488E-08 5.16565923E-06 1.63326661E-03 4.09054093E-01 -5.73848918E-07
1.83284506E-08 5.79596499E-06 1.83248024E-03 5.77439580E-01 3.83528434E-07
2.05648598E-08 6.50317967E-06 2.05597066E-03 4.57925568E-01 -2.56328896E-07
2.30741523E-08 7.29668760E-06 2.30668732E-03 -6.98345034E-02 1.71315859E-07
2.58896247E-08 8.18701814E-06 2.58793429E-03 -7.05728308E-01 -1.14497912E-07
2.90486366E-08 9.18598542E-06 2.90341135E-03 -7.84315028E-01 7.65239831E-08
3.25931064E-08 1.03068452E-05 3.25725924E-03 2.52561836E-01 -5.11443386E-08
3.65700668E-08 1.15644704E-05 3.65410908E-03 1.11765483E+00 3.41820077E-08
4.10322899E-08 1.29755492E-05 4.09913616E-03 7.96337157E-02 -2.28453369E-08
4.60389865E-08 1.45588057E-05 4.59811762E-03 -1.54929662E+00 1.52685420E-08
5.16565924E-08 1.63352486E-05 5.15749377E-03 9.69765096E-01 -1.02046372E-08
5.79596500E-08 1.83284503E-05 5.78443175E-03 5.22594045E-01 6.82020717E-09
6.50317969E-08 2.05648593E-05 6.48688994E-03 -1.38177462E+00 -4.55824398E-09
7.29668762E-08 2.30741515E-05 7.27368023E-03 1.42319844E+00 3.04647464E-09
8.18701817E-08 2.58896236E-05 8.15452380E-03 -1.10540419E+00 -2.03609279E-09
9.18598547E-08 2.90486352E-05 9.14009379E-03 7.66503550E-01 1.36081024E-09
1.03068452E-07 3.25931043E-05 1.02420355E-02 -5.12775503E-01 -9.09489244E-10
1.15644705E-07 3.65700639E-05 1.14729492E-02 3.41609163E-01 6.07851604E-10
1.29755494E-07 4.10322858E-05 1.28463164E-02 -2.28175614E-01 -4.06253922E-10
1.45588058E-07 4.60389807E-05 1.43763378E-02 1.52559552E-01 2.71517338E-10
1.63352488E-07 5.16565843E-05 1.60776460E-02 -1.01971829E-01 -1.81466961E-10
1.83284506E-07 5.79596384E-05 1.79648247E-02 6.81663603E-02 1.21282340E-10
2.05648598E-07 6.50317806E-05 2.00516669E-02 -4.55882998E-02 -8.10583142E-11
2.30741523E-07 7.29668532E-05 2.23500306E-02 3.04810683E-02 5.41748312E-11
2.58896247E-07 8.18701492E-05 2.48681715E-02 -2.03615795E-02 -3.62074189E-11
2.90486366E-07 9.18598088E-05 2.76082657E-02 1.35977270E-02 2.41990082E-11
3.25931064E-07 1.03068387E-04 3.05629294E-02 -9.08987662E-03 -1.61732600E-11
3.65700668E-07 1.15644614E-04 3.37101508E-02 6.08272853E-03 1.08092999E-11
4.10322899E-07 1.29755364E-04 3.70064047E-02 -4.06761271E-03 -7.22432985E-12
4.60389865E-07 1.45587876E-04 4.03767590E-02 2.71487245E-03 4.82833692E-12
5.16565924E-07 1.63352230E-04 4.37019801E-02 -1.81128418E-03 -3.22698961E-12
5.79596500E-07 1.83284142E-04 4.68002258E-02 1.21137650E-03 2.15673888E-12
6.50317969E-07 2.05648083E-04 4.94047004E-02 -8.12064069E-04 -1.44144331E-12
7.29668762E-07 2.30740795E-04 5.11324187E-02 5.43426304E-04 9.63379869E-13
8.18701817E-07 2.58895218E-04 5.14507096E-02 -3.61946442E-04 -6.43869076E-13
9.18598547E-07 2.90484914E-04 4.96316608E-02 2.40820925E-04 4.30325981E-13